Sage Journals: Discover world-class research

Abstract

In this paper, a topology optimization method combining Bi-directional Evolutionary Structural Optimization (BESO) and Fast Non-dominated Sorting Genetic Algorithm II (NSGA-II) is proposed, which is called BESO-NSGA-II. To effectively fuse the two algorithms, a semi-random double-penalized crossover and mutation operator is developed. The proposed algorithm extends the traditional single-objective optimization to a double objective optimization problem, and investigates the influence of the stress norm parameters on the optimization results. To verify the effectiveness of the new method, several examples of 2D structural stiffness and stress optimization are given. At the same time, the optimization design of each model was repeated six times to demonstrate the ability of the new method to converge to a solution. Compared with BESO, BESO-NSGA-II can effectively realize the double objective optimization of stiffness and stress, and eliminate the stress concentration phenomenon on the premise of ensuring that the structure has greater stiffness.

Keywords

Bi-directional Evolutionary Structural Optimization (BESO)fast non-dominated sorting genetic algorithm II (NSGA-II)double-penalized operator stiffness maximization stress minimization double objective optimization

Introduction

Topology optimization, first proposed by Bendsøe and Kikuchi,¹ is an innovative structural optimization method that can quickly find the optimal material layout under specific conditions. After development, topology optimization has been more and more used in various fields. The most commonly used topology optimization methods include Solid Isotropic Material Interpolation Penalty (SIMP),² homogenization method,¹ level set method^3–5 and the Evolutionary Structural Optimization (ESO) method.⁶ For the common checkerboard phenomenon, gray area problem and local optimal solution phenomenon in topology optimization,⁷ various scholars have proposed many effective methods, including sensitivity filtering,⁸ perimeter constraint,⁹ and continuous method¹⁰ etc.

Early topology optimization design problems usually take stiffness as the optimization objective and volume as the constraint variable, however these designs tend to generate large stress concentrations, especially using the BESO algorithm.¹¹ Because the boundary of the optimization result is not smooth, the structure will appear stress concentration, which will cause serious damage to the structure in engineering problems. With the deepening of the research on topology optimization problems, many scholars have discussed the use of stress as the objective function or constraint, and summed up three major challenges of the stress problem: “singularity” problem, the local nature of stress and highly nonlinear stress behavior. In response to these three challenges, various scholars have proposed effective solutions.^12–16 Xia et al.¹⁷ used the BESO algorithm to optimize the design of the stress, using the P-norm stress aggregation scheme to replace the overall stress level, and adjusting the sensitivity number through the filtering scheme, which effectively avoided the three major problems of stress. Khodamorad Nabaki et al.¹⁸ employed a modified P-norm to concentrate the local stresses in a global function and took into account the effects of artificial material interpolation schemes in the sensitivity values. Deng et al.¹⁹ also used the P-norm to achieve an approximate measurement of global stress, and released a code for stress sensitivity analysis based on 3D model. Although the global stress measurement cannot obtain accurate stress values, it can effectively reduce the stress concentration phenomenon in the structure, so the P-norm global stress measurement is widely used in stress optimization problems.

For the optimization of stress and stiffness, Han²⁰ introduced Lagrange multipliers to comprehensively consider the volume constraint and von Mises stress constraint to achieve the design of maximum stiffness while reducing stress concentration. Nicola Ferro et al.²¹ combined the SIMP method with anisotropic adaptive meshes. Under the mass minimization model, taking L-bracket as an example, the influence of stress constraints and stiffness constraints on the optimization results was discussed, and the importance of considering both stress and stiffness constraints on the optimization results is verified. Erik Holmberg et al.²² proposed three optimization models for stiffness, stress and mass by using the SIMP method. Through comparison, it is proved that the simultaneous optimization of stiffness and stress can effectively reduce the stress concentration phenomenon and ensure the stiffness of the structure. However, the above studies all consider stress as a constraint. The multi-objective optimization problem of topology optimization needs to be considered when co-optimizing with stiffness and stress as objectives. In topology optimization problems, multi-objective optimization usually performs weighted calculation, and the multi-objective function values are considered globally in a single dimension through normalization processing.²³ However, determining the weight of each objective function has become a new problem, which leads to a complicated computational process. Furthermore, for classical topology optimization algorithms, the SIMP method performs a single-point optimization search based on sensitivity. The optimization result will show a local optimum to a certain extent. The BESO method may fail in addition and removal of materials in some aspects.²⁴ To accurately perform multi-dimensional optimization on multi-objective function problems, other optimization algorithms are needed for collaborative optimization.

In recent years, Genetic Algorithms (GA)²⁵ have also been developed and gradually applied to topology optimization due to their powerful global search capabilities.²⁶ Aguilar Madeira et al.²⁷ and Madeira et al.²⁸ made great contributions to the application of GA to topology optimization, and they developed a computational model for handling linear elastic structural topologies where multiple objective functions are required optimization. However, this method cannot remove the lowest stress element, nor can it realize the addition or deletion of elements in subsequent iterations, and there may be errors in the deletion of elements. The combination of GA with traditional optimization algorithms has also been extensively studied. The combination of GA and ESO²⁴ greatly improves the performance of the ESO algorithm, but material addition cannot be achieved during optimization. Zuo et al.²⁹ proposed topology optimization methods for GA and BESO. It implements material addition and removal in an iterative process, but remains stuck in a single-objective optimization problem. NSGA-II is developed on the basis of GA³⁰ and improved on the basis of the original GA. The NSGA-II algorithm mainly includes operators such as fast non-dominated sorting and crowding calculation, competitive selection and elite strategy to achieve fast optimization. For multi-objective global optimization problems, NSGA-II has greater potential and advantages than traditional genetic algorithms.

Aiming at the multi-dimensional optimization problem of von Mises stress and stiffness, in order to improve the global search ability of the optimization algorithm, a joint algorithm combining NSGA-II and BESO is proposed, called BESO-NSGA-II. This algorithm aims to reduce the maximum stress while maximizing the stiffness of the structure under the premise of satisfying the volume constraints. P-norm global stress aggregation scheme was used as an approximation to the maximum stress. A double penalty operator is proposed to ensure a “clean” topology and improve the stability and convergence of the algorithm. Finally, the effectiveness of this algorithm in double objective optimization is proved by several examples.

The arrangement of this paper is as follows. Section 2 presents design models for stress minimization and stiffness maximization. Section 3 presents the derivation and filtering scheme of the sensitivity. Section 4 introduces the specific implementation of the BESO-NSGA-II algorithm. Section 5 verifies the effectiveness of the proposed method through a series of 2D examples. The last section gives the conclusion.

Problem statements

The Classical topology optimization problem pursues the maximum stiffness under a certain volume fraction. In the convex optimization problem, the maximum value problem is usually transformed into a minimum value problem. At this time, the maximum stiffness problem can be equivalently replaced by the minimum mean compliance problem can be written as:

{\begin{matrix} find : {x_{i}} (i = 1, 2, . . ., N) \\ \min : C = \frac{1}{2} U^{T} KU \\ s . t : {\begin{matrix} KU = F \\ V^{*} - \sum_{i = 1}^{N} V_{i} x_{i} = 0, x_{i} = 1 or x_{\min} \end{matrix} \end{matrix}

(1)

where C is the mean compliance, x_i is the design variable of the model, V^* is the volume constraint, V_i is the volume of the ith element. There are two constraints in the formula, the first constraint is equilibrium constraint of the structure, and the second constraint is volume constraint. K is the global stiffness matrix. F and U represent the external force vector and displacement vector that the structure is subjected to, and N represents the number of discrete elements. x_i = 1 when the ith element is filled with solid element, otherwise x_i = x_min = 1.0 ×10⁻⁶.

For the optimization results with the average compliance as the objective function, there is often a large stress concentration phenomenon, such as the reentrant corner of the L-bracket. So another solution is to seek to minimize the stress of the structure. The mathematical model that minimizes the P-norm global stress measure can be written as:

{\begin{matrix} find : {x_{i}} (i = 1, 2, . . ., N) \\ \min : σ_{PN} \\ s . t : {\begin{matrix} KU = F \\ V^{*} - \sum_{i = 1}^{N} V_{i} x_{i} = 0, x_{i} = 1 or x_{\min} \end{matrix} \end{matrix}

(2)

For the double objective optimization problem of stiffness and stress, the advantage of the BESO-NSGA-II algorithm is that the objective functions in equations (1) and (2) can be simply combined. There is also no need to normalize the target values during the sensitivity value ranking stage. Only by the fast non-dominated sorting operator in the BESO-NSGA-II algorithm, the sensitivity of the cells can be distinguished.

For the stress-based topology optimization problem, the effective stress vector model is expressed as follows³¹:

σ_{i} = D_{i} B_{i} u_{i}

(3)

where D_i represents the stiffness matrix of the solid material, u_i represents the nodal displacement vector of the ith element, and B_i represents the strain matrix of the ith element. According to the expression, the stress of the element is independent of the design variables. For two-dimensional plane stress problem, the von Mises stress can be defined as:

σ_{vm, i} = {(σ_{ix}^{2} + σ_{iy}^{2} - σ_{ix} σ_{iy} + 3 τ_{ixy}^{2})}^{\frac{1}{2}}

(4)

For the second biggest challenge of stress: the local nature of stress, here the P-norm global stress is used to measure the equivalent maximum stress:

σ_{PN} = {(\sum_{i = 1}^{N} σ_{vm, i}^{pl})}^{\frac{1}{pl}}

(5)

where σ_vm,i represents the von Mises stress at the centroid of the ith element. N represents the number of total elements, and pl is the stress norm parameter. When pl = 1, represents the mean stress, and when pl approaches infinity, is approximately equal to the maximum stress. However, the study found that when the value of pl is 3–5, the model is suitable for most topology optimization models.¹⁸

Sensitivity numerical calculation and filter scheme

Assuming that the Poisson’s ratio is independent of the design variables, the stiffness matrix K can be expressed as the element stiffness matrix and the design variables as³²:

K = \sum_{i} x_{i}^{p} K_{i}

(6)

where K_i represents the element stiffness matrix of the solid element. In the application process, in order to avoid the singularity of the stiffness matrix, an extremely low stiffness is assigned to the cavity element. Update design variables via gradient-based methods. Usually, the adjoint method is used to solve the gradient of the objective function to a single element.³³ For the minimum compliance solid-cavity optimization problem, the element sensitivity value can be obtained by deriving the objective function.

α_{i} = \frac{\partial C}{\partial x_{i}} = - \frac{1}{2} {px}_{i}^{p - 1} u_{i}^{T} K_{i} u_{i}

(7)

For the sensitivity calculation of the P-norm global stress, according to the chain rule:

\frac{\partial σ_{PN} (x)}{\partial x_{i}} = \frac{\partial σ_{PN} (x)}{\partial σ_{vm, i}} \frac{\partial σ_{vm, i} (x)}{\partial x_{i}} = \frac{\partial σ_{PN} (x)}{\partial σ_{vm, i}} (\frac{\partial σ_{vm, i} (x)}{\partial σ_{i}}) \frac{\partial σ_{i} {(x)}^{T}}{\partial x_{i}}

(8)

Taking the derivation of equation (5), we can get:

\frac{\partial σ_{PN} (x)}{\partial σ_{vm, i}} = {(\sum_{i = 1}^{nele} σ_{vm, i}^{pl})}^{(\frac{1}{pl} - 1)} {(σ_{vm, i}^{pl} (x))}^{pl - 1}

(9)

For the two-dimensional plane model, based on the definition of von Mises stress in equation (4), its derivative with respect to $σ_{i}$ can be expressed as:

\frac{\partial σ_{vm, i} (x)}{\partial σ_{ix}} = \frac{1}{2 σ_{vm, i} (x)} (2 σ_{ix} (x) - σ_{iy} (x))

(10)

\frac{\partial σ_{vm, i} (x)}{\partial σ_{iy}} = \frac{1}{2 σ_{vm, i} (x)} (2 σ_{iy} (x) - σ_{ix} (x))

(11)

\frac{\partial σ_{vm, i} (x)}{\partial τ_{ixy}} = \frac{3}{σ_{vm, i} (x)} τ_{ixy} (x)

(12)

In equation (3), the stress vector can be derived from the design variable x_i to obtain:

\frac{\partial σ_{i} (x)}{\partial x_{i}} = \frac{\partial D (x)}{\partial x_{i}} Bu + DB \frac{\partial u (x)}{\partial x_{i}}

(13)

\frac{\partial F}{\partial x_{i}} = \frac{\partial K (x)}{\partial x_{i}} u + K \frac{\partial u (x)}{\partial x_{i}}

(14)

Assuming that the specified load vector is independent of the design variables, equation (12) is equal to 0 and can be written as:

\frac{\partial u (x)}{\partial x_{i}} = - K^{- 1} \frac{\partial K (x)}{\partial x_{i}} u

(15)

Substituting equation (15) into equation (8), we can get:

\begin{matrix} \frac{\partial σ_{PN} (x)}{\partial x_{i}} = \frac{\partial σ_{PN} (x)}{\partial σ_{vm, i}} \frac{\partial σ_{vm, i} (x)}{\partial x_{i}} \\ = \frac{\partial σ_{PN} (x)}{\partial σ_{vm, i}} {(\frac{\partial σ_{vm, i} (x)}{\partial σ_{i}})}^{T} \times (\frac{\partial D (x)}{\partial x_{i}} Bu - DB K^{- 1} \frac{\partial K (x)}{\partial x_{i}} u) \end{matrix}

(16)

The accompanying variable is introduced, which can be defined as:

λ^{T} = \frac{\partial σ_{PN} (x)}{\partial σ_{vm, i}} {(\frac{\partial σ_{vm, i} (x)}{\partial σ_{i}})}^{T} DB K^{- 1}

(17)

Substituting equation (17) into equation (16), we can get:

\frac{\partial σ_{PN} (x)}{\partial x_{i}} = \frac{\partial σ_{PN} (x)}{\partial σ_{vm, i}} {(\frac{\partial σ_{vm, i} (x)}{\partial σ_{i}})}^{T} \frac{\partial D (x)}{\partial x_{i}} Bu - λ^{T} (\frac{\partial K (x)}{\partial x_{i}} u)

(18)

Then the sensitivity can be expressed as:

α_{i} = - [\frac{\partial σ_{PN} (x)}{\partial σ_{vm, i}} {(\frac{\partial σ_{vm, i} (x)}{\partial σ_{i}})}^{T} \frac{\partial D (x)}{\partial x_{i}} Bu - λ^{T} (\frac{\partial K (x)}{\partial x_{i}} u)]

(19)

It should be noted that the sensitivity number only exists in solid elements, and the sensitivity for hollow elements is directly set to 0. Sensitivity is corrected by a sensitivity filtering strategy. This not only effectively reduces the checkerboard phenomenon, but also enables the addition and deletion of elements by deriving the sensitivity of empty elements by filtering the sensitivity values of physical elements within the radius. The filtering strategy adopted in this paper is divided into two steps to complete. The first step is to equalize the sensitivity values of all elements associated with the node to the element sensitivity values of the node,¹¹ as shown in the following equation:

α_{i}' = \frac{\sum_{i = 1}^{H} w (r_{ik}) α_{k}^{n}}{\sum_{i = 1}^{H} w (r_{ik})}

(20)

In the formula, H is the total number of nodes in the design domain, w(r_ik) is the weight factor, r_ik is the distance from the center of node i to the center of node k, and the weight factor can be expressed as:

w (r_{ik}) = {\begin{matrix} r_{min} - r_{ik} for r_{ik} < r_{min} \\ 0 for r_{ik} \geq r_{min} \end{matrix}

(21)

The second step uses a convergence strategy that considers historical iteration information to correct the sensitivity value. The specific operation is to calculate the average value of the sensitivity of the current iteration and the previous iteration as the update value of the sensitivity.¹¹

\bar{α_{i}} = \frac{1}{2} (α_{i, g}' + α_{i, g - 1}')

(22)

In the formula, $α_{i, g}'$ and $α_{i, g - 1}'$ are the sensitivity values of the gth and (g − 1)th iteration steps of the ith element, respectively. Finally, the effect of each element on the structural properties is measured according to this sensitivity value.

During the optimization process, the material redistribution and the update calculation of the objective function are continuously carried out. When the optimization reaches the specified target mass fraction and the following convergence conditions are met, the optimization process is terminated. The convergence error is mainly calculated from the objective function value.¹¹

error = \frac{| \sum_{i = 1}^{N} (M_{k - i + 1} - M_{k - N - i + 1}) |}{\sum_{i = 1}^{N} M_{k - i + 1}} \leq τ

(23)

In the formula, M represents the value of the objective function, k is the current number of iterations, and N is 5, indicating that the optimization results of five consecutive iterations have a certain stability. $τ$ is the allowable convergence error. In this paper, it is set to 0.1%.

Implementation of BESO-NSGA-II

The main steps of traditional NSGA-II optimization are initialization population, objective function calculation, fast non-dominated sorting and crowding degree calculation, competition selection, crossover mutation, and elite selection strategy. Solve the optimization problem by iteratively. For the BESO-NSGA-II, the main difference from the traditional NSGA-II lies in the crossover and mutation operator and the elite selection strategy. The details are as follows.

Pre-processing process

The initial optimization model is pre-processed, the corresponding mechanical model is established according to the actual application conditions, and the load and boundary conditions are set. The finite element discretization process is performed on the structure, and the number of finite element meshes is used as the initial population number of the BESO-NSGA-II algorithm. Consider each element as an independent individual for analysis. In order to ensure that the number of individuals in the entire population is always equal to the number of elements in the design domain, the number of populations is set to remain unchanged during the iterative process. The gene of each individual is represented by a binary string of length 6. The initial string is defined in this paper as an array of all 1 s, indicating that all elements are solid-phase materials, and an array of all 0 s means that the element to be deleted in this iteration.

Fast non-dominated sorting and crowding degree calculation

Perform finite element analysis on the structure, and calculate the sensitivity value of each element according to equations (7) and (19). The stiffness and stress sensitivity values are matched for each individual and used as the basis for sorting. Traverse all individuals in the population, if all the sensitivity values of individual i are not worse than that of individual j, and at least one sensitivity value of individual i is better than that of individual j, then individual i is said to be in the dominant position. Otherwise, individual i is marked as a non-dominated individual. Similarly, the comparison of sensitivity values is continued for all non-dominant individuals until the entire population is ranked according to dominance in the non-dominant position. Individuals with high stress and stiffness sensitivities are ranked higher in the ranking and will be retained in subsequent iterations. Individuals with low-level sensitivity values are ranked lower, and have a high probability of being eliminated, which means deleting elements. For individuals of the same level, the crowding degree is calculated, and the expression is:

cd = {\begin{matrix} \infty α^{max} = α^{min} \\ \frac{α_{i + 1} - α_{i - 1}}{α^{max} - α^{min}} α^{max} \neq α^{min} \end{matrix}

(24)

In the formula, α^max and α^min represent the maximum and minimum values of sensitivity in this iteration cycle. α_i_+ 1 represents the sensitivity value calculated by the (i + 1)th element, and α_i_-1 represents the sensitivity value calculated by the (i − 1)th element.

Cross and mutation operator

The crossover mutation operator plays an extremely important role in NSGA-II, which enriches the population diversity. The crossover mutation operator proposed in this paper adopts a semi-random method, which avoids the checkerboard phenomenon and non-convergence under the fully random operator to a certain extent. The cross-mutation operator only has the transformation between 0 and 1.

The double penalty operator is mainly used in the crossover mutation operator. First, the population is divided into high-sensitivity regions and low-sensitivity regions according to the sensitivity thresholds calculated for each generation. Two individuals are randomly selected as the parents in the design domain, and two crossover positions are randomly selected to perform the crossover operation on the parents. In order to speed up the convergence, the algorithm gradually penalizes the crossover rate to a certain maximum value as the number of iterations increases during the optimization process. The mutation operator is to mutate the binary code of each individual according to a certain mutation rate. Here, a double penalty operator is proposed. In the iterative loop of the optimization algorithm, with the increase of the number of iterations, the crossover rate and the mutation rate are gradually penalized to a certain maximum value, and the maximum value is limited to 1, such as equations (25) and (26), where p_c is crossover rate, p_c, _min is minimum crossover rate, p_c, _max is maximum crossover rate, p_m mutation rate, p_m, _min is minimum mutation rate, p_m, _max is maximum mutation rate, prg is a metric that controls the rate of crossover variation. ph is the penalty factor to control the crossover mutation rate, and the value is 3 in this paper.

p_{c} = p_{c, min} + (p_{c, max} - p_{c, min}) pr g^{ph}

(25)

p_{m} = p_{m, min} + (p_{m, max} - p_{m, min}) pr g^{ph}

(26)

In each iteration step, for the high-sensitivity region, as shown in equation (27), the mutation rate (probability of binary code mutation to 1) is gradually penalized to a certain minimum value as the sensitivity value increases from high to low. Individuals with higher sensitivity values have a greater probability of having more 1 s, which are easier to retain in the optimization process. For low-sensitivity regions, as shown in equation (28), as the sensitivity increases from high to low, the mutation rate (The probability that the binary code is mutated to 0.) is penalized to a certain maximum value. At this time, individuals with lower sensitivity values have a greater probability of having more 0 s and are more likely to be deleted in the optimization process. It is worth noting that individuals of the same rank have the same rate of variation. Through the double penalty mechanism, the fast convergence of the optimization process and a clearer and more stable topology are ensured.

p_{high, m} = p_{high, m} - \frac{p_{m, max} - p_{m, min}}{r_{high, max} - r_{high, min} + 1}

(27)

p_{low, m} = p_{low, m} + \frac{p_{m, max} - p_{m, min}}{r_{low, max} - r_{low, min} + 1}

(28)

In the formula, p_high,m is the mutation rate of high-sensitivity individuals, p_low,m is the mutation rate of low-sensitivity individuals, r_high,max is the maximum individual level in the high-sensitivity region, and r_high,min is the minimum individual level in the high-sensitivity region. In the same way, r_low,max and r_low,min are the maximum and minimum values of individual levels in the low-sensitivity area. The initial value of p_high,m is p_m,max, and the initial value of p_low,m is p_m.

Selection mechanism

For the traditional NSGA-II algorithm, through the elite selection strategy, the elite individuals in the parent and child are selected to form a new population. In this process, the elite individuals in the offspring will directly replace the inferior individuals in the parent. For BESO-NSGA-II, the layout of solid materials and void materials in the structure needs to be optimally designed, which means that the poor individuals in the parent generation also need to be preserved. Therefore, a selection mechanism for this optimization algorithm is proposed. Because the number of 0 s in an individual reflects the quality of the individual. Therefore, in the iterative process, the individuals whose binary codes are all 0 are assigned to hollow materials, and the rest of the individuals are saved for the subsequent optimization iterative process, which can effectively prevent the accidental deletion of elements to a certain extent. Figure 1 shows the flow chart of the algorithm implementation:

Figure 1.

Algorithm implementation flowchart.

Numerical examples

In this section, we will demonstrate the effectiveness of our method with several two-dimensional examples. For all cases, it is assumed that the Young’s modulus of the solid material is 1.0 MPa, the Poisson’s ratio is 0.3, the volume evolution ratio ER = 0.02, and the filter radius is about one to three times of the size of one element.¹¹ For all cases, filter radius r_min = 2.5. The convergence criterion is calculated according to equation (23).

Cantilever beams

This example considers the cantilever beam model shown in Figure 2. The right middle node is subjected to a downward load of 5 N. In order to eliminate the influence of local stress concentration, the load is evenly distributed on five nodes. Set the volume fraction to 30%.

Figure 2.

Cantilever beam problem.

In order to verify the effectiveness of the double penalty operator proposed in this paper, the algorithm with the double penalty operator added and the algorithm without the double penalty operator are compared, as shown in Figure 3. The results show that the optimization results of the two algorithms are almost the same. However, the algorithm adding the double penalty operator can converge quickly, and the convergence result can be obtained in 160 generations. The algorithm without the double penalty operator needs to go through 366 iterations to get the optimal solution. Therefore, the double penalty operator is introduced in the examples in this paper to ensure the fast convergence of the results.

Figure 3.

Algorithm comparison results.

Due to the randomness of NSGA-II results, there is no guarantee that the same results will be obtained every time, but fluctuate up and down within the optimal solution range. In order to verify the stability of the proposed algorithm, six repeated experiments were carried out on the model under the same parameters. Figure 4 shows the results of six iterations of the optimization when pl = 3. Similar results can be obtained for the material layout of the structure and the distribution of stress. The two objective functions shown in Table 1 can take close values, and there is no obvious fluctuation around the average value, which verifies the stability of the proposed algorithm. Figure 5 shows the evolution of topology optimization for two objective functions and volumes. It is proved that the algorithm can obtain convergent solutions under both stiffness and stress objectives.

Figure 4.

Six optimization results under the stress norm parameter pl = 3.

Table 1.

Comparison of six repeated optimization results.

	(a)	(b)	(c)	(d)	(e)	(f)
Mean compliance	1230.25	1229.94	1228.58	1232.78	1231.86	1231.08
Maximum stress	2.15	2.13	2.16	2.18	2.14	2.12

Figure 5.

Topology optimization evolution of cantilever beams: (a) evolution curve of volume versus mean compliance and (b) evolution curve of volume and P-norm stress.

In order to verify the effectiveness of the algorithm proposed in this paper in the double objective optimization problem, the optimization results under the same parameters are compared with the results of the classical algorithm BESO optimization. Table 2 shows the comparison results of the two algorithms. For the classical BESO optimization algorithm, when the value of the P-norm stress parameter pl is 3–5, a relatively stable solution will be obtained, and if the value of pl is greater than 6, serious oscillations will occur and the design will fail.¹⁷ Through testing, this algorithm can extend pl to 8 and ensure the stability of the structure. It can be obtained from the analysis that with the increase of pl value, the mean compliance shows an upward trend, and the maximum stress value of the structure gradually decreases. Comparing the objective function values at pl = 3 and 8, the mean compliance increased by 2.88%, while the stress value decreased by 23.26%. It shows that the algorithm can effectively reduce the maximum stress value of the structure under the condition that the influence on the structural stiffness is small.

Table 2.

Comparison of BESO-NSGA-II and traditional BESO results.

BESO-NSGA-II considering stress and compliance
pl = 3	pl = 4	pl = 6	pl = 8


BESO considering stress and compliance
BESO considering stress	pl = 3	pl = 4	pl = 5


BESO considering compliance

Table 2 shows the optimal design results for structural stiffness and stress using the classical BESO optimization algorithm. Table 3 is a comparison of the optimization target values corresponding to the two algorithms. Compared with the stress design using the classical BESO, the maximum stress of the optimization result of the BESO-NSGA-II is 1.65, which is better than the 1.79 calculated by the BESO, and the mean compliance value of 1266.78 is also smaller than that of the BESO of 1286.47. Compared with the stiffness design using the classical BESO, the mean compliance value of the BESO-NSGA-II is higher, the increase ratio is only 1.75%, which means that the structural stiffness is reduced by 1.75%, but the stress is reduced by 37.02%. Therefore, compared with the BESO, this algorithm has certain advantages in double objective optimization. It is worth noting that this algorithm pursues the equilibrium state of the optimal solution achieved by the two objectives, so there is no need to guarantee that the target values calculated by BESO-NSGA-II are all better than the target values of the stress design and stiffness design performed by BESO.

Table 3.

Optimization objective function value comparison.

BESO-NSGA-II considering stress and compliance
	pl = 3	pl = 4		pl = 6	pl = 8
Mean compliance	1230.25	1232.71		1264.98	1266.78
Maximum stress	2.15	2.04		1.76	1.65
BESO considering stress and compliance
BESO considering stress				BESO considering compliance
	pl = 3	pl = 4	pl = 5
Mean compliance	1245.02	1249.35	1286.47	mean compliance	1244.64
Maximum stress	2.21	1.90	1.79	maximum stress	2.62

L-bracket design

As a classic example of stress optimization problem, L-bracket has been studied by many scholars.^17,18 This example considers the optimal design of an L-shaped beam with a length of 80 mm. As shown in Figure 6, the top of the L-shaped beam is fixed, and the right end is subjected to a vertical downward size of 3 N. In order to avoid stress concentration, the load is distributed on three nodes. The volume fraction is 25%.

Figure 6.

Design domain and support conditions for L-bracket.

As shown in Figure 7, in order to verify the stability of the algorithm, the same model was optimized six times under the condition of pl = 3. It can be observed that the layout of materials and the distribution of stress are highly similar in the six optimization results. Table 4 shows the objective function values of the six optimization results. It can be noticed that all the values are relatively close, and there is no obvious deviation near the average value, which verifies the stability of the proposed algorithm. Figure 8 shows the corresponding optimization iteration curve. It can be seen that all objectives can converge to a stable solution.

Figure 7.

Six optimization results under the stress norm parameter pl = 3.

Table 4.

Comparison of six repeated optimization results.

	(a)	(b)	(c)	(d)	(e)	(f)
Mean compliance	1230.25	1229.94	1228.58	1232.78	1231.86	1231.08
Maximum stress	2.15	2.13	2.16	2.18	2.14	2.12

Figure 8.

Topology optimization evolution of L-bracket: (a) evolution curve of volume versus mean compliance and (b) evolution curve of volume and P-norm stress.

Table 5 compares the results between BESO-NSGA-II and traditional BESO. It can be observed that under the L-bracket model, the pl value can also be extended to 8. With the increase of the pl value, the stress concentration of the L-bracket is gradually reduced, the material at the reentrant corner is also gradually removed, and the stress concentration at the tip disappears. At the same time, Table 5 also shows the optimization of the L-bracket model using the BESO method. It is worth noting that the scope of use of pl at this time is consistent with that described in Nabaki et al.¹⁸

Table 5.

Comparison of BESO-NSGA-II and traditional BESO results.

BESO-NSGA-II considering stress and compliance
pl = 3	pl = 4	pl = 6	pl = 8


BESO considering stress and compliance
BESO considering stress	pl = 3	pl = 4	pl = 5


BESO considering compliance

As shown in Table 6, compared with the stress design using BESO, the optimization result of BESO-NSGA-II achieves the minimum stress when pl = 8, and its value is 1.63, which is 12.37% lower than that calculated by BESO, and the corresponding mean compliance value decreased by 8.58%. Compared with the stiffness design using BESO, the mean compliance value of BESO-NSGA-II increased by 15.32%, but decreased by 34.27% for stress. Khodamorad¹⁸ also conducted an optimization study on the stress and stiffness of the L-bracket, and the compliance increased by 17%–25% under the condition that the stress was reduced by 25%–37%. For literature,^34,35 optimization is performed with stress as a constraint, constraining the structural stress value below a specified value. There is no guarantee that the stress level of the structure will always be at the optimal solution. Therefore, the currently proposed optimization algorithms can provide better solutions in double objective optimization.

Table 6.

Optimization objective function value comparison.

BESO-NSGA-II considering stress and compliance
	pl = 3	pl = 4		pl = 6	pl = 8
Mean compliance	947.13	1018.72		1056.27	1106.96
Maximum stress	2.38	2.36		1.93	1.63
BESO considering stress and compliance
BESO considering stress				BESO considering compliance
	pl = 3	pl = 4	pl = 5
Mean compliance	1016.65	1096.40	1210.86	mean compliance	937.40
Maximum stress	2.29	2.17	1.86	maximum stress	2.48

Conclusion

In this paper, a topology optimization method combining BESO and NSGA-II is proposed, which realizes the double objective optimization design of stress minimization and stiffness maximization. Due to the discrete nature of the proposed algorithm, the “singularity” problem in the stress optimization problem is avoided. The locality of stress and the high degree of nonlinearity are addressed through P-norm global stress measurements, filtering schemes, and calculation of sensitivity numbers. This method theoretically avoids three major challenges in stress optimization problems. Due to the influence of the algorithm of NSGA-II, the algorithm has a certain randomness. After six repeated designs for each 2D example, highly similar results can be obtained, which proves the stability of the algorithm. The implementation of the double penalty operator enables the algorithm to obtain a convergent and clean topology. Compared with the single objective optimization problem, the proposed algorithm can effectively achieve a balance between the stiffness optimization and the stress optimization, and can not only improve the stiffness of the structure but also reduce the maximum stress. This algorithm provides theoretical support for multi-objective optimization problems and is expected to be extended to three-dimensional structural optimization problems. In addition, this algorithm is also expected to be extended to the optimal design of two-scale lattice structures to reduce stress concentration on the basis of improving structural stiffness.

Footnotes

Handling Editor: Chenhui Liang

Authors’ contributions

Shuo Feng: Revised the article critically for professional knowledge content. Xing Wu: Proposed a new joint algorithm, drafted this manuscript. Jun Qi: Participated in the experiment of this manuscript. Qidong Han: Participated in the experiment of this manuscript. Zhi He: Revised the article critically for professional knowledge content.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China [12102065], the Fundamental Research Funds for the Central Universities [300102258113].

ORCID iD

Shuo Feng

Code availability

Custom code

References

Bendsøe

Kikuchi

Generating optimal topologies in structural design using a homogenization method. Comput Methods Appl Mech Eng 1988; 71: 197–224.

Bendsøe

MP.

Optimal shape design as a material distribution problem. Struct Optim 1989; 1: 193–202.

Allaire

Jouve

Toader

AM.

Structural optimization using sensitivity analysis and a level-set method. J Comput Phys 2004; 194: 363–393.

Sethian

Wiegmann

Structural boundary design via level set and immersed interface methods. J Comput Phys 2000; 163: 489–528.

Wang

Guo

A level set method for structural topology optimization. Comput Methods Appl Mech Eng 2003; 192: 227–246.

Huang

Xie

YM.

A further review of ESO type methods for topology optimization. Struct Multidiscipl Optim 2010; 41: 671–683.

Zhou

Shyy

Thomas

HL.

Checkerboard and minimum member size control in topology optimization. Struct Multidiscipl Optim 2001; 21: 152–158.

Sigmund

Design of material structures using topology optimization. Lyngby: Technical University of Denmark, 1994.

Ambrosio

Buttazzo

An optimal design problem with perimeter penalization. Calc Var Partial Differ Equ 1993; 1: 55–69.

10.

Borrvall

Topology optimization of elastic continua using restriction. Arch Comput Methods Eng 2001; 8: 351–385.

11.

Huang

Xie

YM.

Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method. Finite Elem Anal Des 2007; 43: 1039–1049.

12.

Bruggi

On an alternative approach to stress constraints relaxation in topology optimization. Struct Multidiscipl Optim 2008; 36: 125–141.

13.

Bruggi

Duysinx

Topology optimization for minimum weight with compliance and stress constraints. Struct Multidiscipl Optim 2012; 46: 369–384.

14.

Duysinx

Sigmund

New developments in handling stress constraints in optimal material distribution. In: 7th AIAA/USAF/NASA/ISSMO symposium on multidisciplinary analysis and optimization, 1998, pp. 4906 . St. Louis: American Institute of Aeronautics and Astronautics, Inc.

15.

Norato

Bruns

, et al. Stress-based topology optimization for continua. Struct Multidiscipl Optim 2010; 41: 605–620.

16.

Yang

Chen

CJ.

Stress-based topology optimization. Struct Optim 1996; 12: 98–105.

17.

Xia

Zhang

Xia

, et al. Stress-based topology optimization using bi-directional evolutionary structural optimization method. Comput Methods Appl Mech Eng 2018; 333: 356–370.

18.

Nabaki

Shen

Huang

Stress minimization of structures based on bidirectional evolutionary procedure. J Struct Eng 2019; 145: 04018256.

19.

Deng

Vulimiri

AC.

An efficient 146-line 3D sensitivity analysis code of stress-based topology optimization written in MATLAB. Eng Optim 2022; 23: 1733–1757.

20.

Han

Duan

, et al. Stress-based topology optimization of continuum structures for the elastic contact problems with friction. Struct Multidiscipl Optim 2022; 65: 54–658.

21.

Ferro

Micheletti

Perotto

Compliance–stress constrained mass minimization for topology optimization on anisotropic meshes. SN Appl Sci 2020; 2: 1–11.

22.

Holmberg

Torstenfelt

Klarbring

Stress constrained topology optimization. Struct Multidiscipl Optim 2013; 48: 33–47.

23.

Teimouri

Asgari

Multi-objective BESO topology optimization for stiffness and frequency of continuum structures. Struct Eng Mech 2019; 72: 181–190.

24.

Liu

, et al. Genetic evolutionary structural optimization. J Constr Steel Res 2008; 64: 305–311.

25.

Vaissier

Pernot

Chougrani

, et al. Genetic-algorithm based framework for lattice support structure optimization in additive manufacturing. Comput Des 2019; 110: 11–23.

26.

Goldberg

Samtani

MP.

Engineering optimization via genetic algorithm. In: Electronic computation, proceedings of the 9th conference on electronic computation, Birmingham, UK, 23–26 February 1986, pp. 471–482. Reston: American Society of Civil Engineers.

27.

Aguilar Madeira

Rodrigues

Pina

Multi-objective optimization of structures topology by genetic algorithms. Adv Eng Softw 2005; 36: 21–28.

28.

Madeira

JFA

Pina

Rodrigues

HC.

GA topology optimization using random keys for tree encoding of structures. Struct Multidiscipl Optim 2010; 40: 227–240.

29.

Zuo

Xie

Huang

Combining genetic algorithms with BESO for topology optimization. Struct Multidiscipl Optim 2009; 38: 511–523.

30.

Deb

Agrawal

Pratap

, et al. A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II. In: Parallel problem solving from nature. Berlin, Heidelberg: Springer, 2000, pp.849–858.

31.

Duysinx

Bendsøe

MP.

Topology optimization of continuum structures with local stress constraints. Int J Numer Methods Eng 1998; 43: 1453–1478.

32.

Huang

Xie

YM.

Evolutionary topology optimization of continuum structures with an additional displacement constraint. Struct Multidiscipl Optim 2010; 40: 409–416.

33.

Bendsoe

Sigmund

Topology optimization: theory, methods, and applications. New York: Springer Science & Business Media, 2003.

34.

Fan

Xia

Lai

, et al. Evolutionary topology optimization of continuum structures with stress constraints. Struct Multidiscipl Optim 2019; 59: 647–658.

35.

Amir

Efficient stress-constrained topology optimization using inexact design sensitivities. Int J Numer Methods Eng 2021; 122: 3241–3272.

Novel multi-objective topology optimization method for stiffness and stress of continuum structures

Abstract

Keywords

Introduction

Problem statements

Sensitivity numerical calculation and filter scheme

Implementation of BESO-NSGA-II

Pre-processing process

Fast non-dominated sorting and crowding degree calculation

Cross and mutation operator

Selection mechanism

Numerical examples

Cantilever beams

L-bracket design

Conclusion

Footnotes

Authors’ contributions

Declaration of conflicting interests

Funding

ORCID iD

Code availability

References