Bi-Directional Evolutionary Topology Optimization with Adaptive Evolutionary Ratio for Nonlinear Structures

Abstract

Current topology optimization methods for nonlinear continuum structures often suffer from low computational efficiency and limited applicability to complex nonlinear problems. To address these issues, this paper proposes an improved bi-directional evolutionary structural optimization (BESO) method tailored for maximizing stiffness in nonlinear structures. The optimization program is developed in Python and can be combined with Abaqus software to facilitate finite element analysis (FEA). To accelerate the speed of optimization, a novel adaptive evolutionary ratio (ER) strategy based on the BESO method is introduced, with four distinct adaptive ER functions proposed. The Newton-Raphson method is utilized for iteratively solving nonlinear equilibrium equations, and the sensitivity information for updating design variables is derived using the adjoint method. Additionally, this study extends topology optimization to account for both material nonlinearity and geometric nonlinearity, analyzing the effects of various nonlinearities. A series of comparative studies are conducted using benchmark cases to validate the effectiveness of the proposed method. The results show that the BESO method with adaptive ER significantly improves the optimization efficiency. Compared to the BESO method with a fixed ER, the convergence speed of the four adaptive ER BESO methods is increased by 37.3%, 26.7%, 12% and 18.7%, respectively. Given that Abaqus is a powerful FEA platform, this method has the potential to be extended to large-scale engineering structures and to address more complex optimization problems. This research proposes an improved BESO method with novel adaptive ER, which significantly accelerates the optimization process and enables its application to topology optimization of nonlinear structures.

Introduction

Topology optimization is a mathematical method used to optimize material distribution within a design domain based on specific load cases, constraints, and performance metrics. In Ref. [1], topology optimization has made significant advancements in addressing linear problems and has achieved considerable success in structural design. However, topology optimization presents substantial challenges due to the large number of design variables and complex constraints, which often make it more difficult to find optimal solutions compared to other optimization methods. Despite these difficulties, the bi-directional evolutionary structural optimization (BESO) method has emerged as a robust tool in both academic research and engineering applications, enabling efficient exploration of optimal structural configurations through iterative addition and removal of elements [2].

Practical engineering scenarios often involve nonlinear phenomena. To achieve more realistic designs, researchers continuously explore the extension of topology optimization to address various nonlinear aspects in structural design, such as material nonlinear, geometrically nonlinear. Xu et al. [3, 4] investigated the BESO method for material and geometric nonlinearity under stress constraints, while considering both strength and stiffness. Banh et al. [5] proposed a topology optimization method that simultaneously considers geometric nonlinearity and design-dependent pressure loads. Wang et al. [6] proposed a topology optimization method considering both geometrical and load nonlinearities. Shobeiri [7] examined the optimal topology of structures subjected to dynamic loading, considering material nonlinearity, geometric nonlinearity, and contact nonlinearity using the BESO method. Chen et al. [8] conducted research on nonlinear topology optimization for flexoelectric soft dielectrics undergoing large deformations. Da Silva et al. [9] demonstrated that the maximum output displacement method, based on nonlinear analysis, provides better solutions under large displacements while simultaneously satisfying stress and manufacturing requirements.

The primary challenge in nonlinear structural topology optimization lies in finding the optimal topology structure or, in other words, developing a robust methodology. In nonlinear structural topology optimization, particularly for structures experiencing large deformations, numerical instability often arises due to excessive deformation in low-stiffness regions. Numerical instability is typically characterized by three key features: excessive distortion, low stiffness regions, and non-positive definite tangent stiffness matrix. Researchers have developed various techniques to address these challenges, such as the additive hyperelasticity technique [10, 11], hyperelastic constitutive material model [12], self-adaptive material interpolation scheme [13], linear material interpolation scheme [14], stress relaxation strategy [15], element removal strategies and the redundant degrees of freedoms removal technique [16, 17, 18, 19, 20–21]. Additionally, some studies have demonstrated that energy interpolation method can effectively address excessive distortion and numerical instability in low-stiffness regions in topology optimization [22, 23–24]. Furthermore, the BESO method with discrete formulation has shown strong numerical stability.

Another major challenge in nonlinear structural topology optimization is the extensive computation time. The computation time for nonlinear topology optimization is significantly influenced by the size of the structure, especially in three-dimensional (3D) optimization problems, where computational efficiency becomes crucial. To address this issue, researchers have focused on optimizing algorithms to improve the efficiency of topology optimization. Sigmund [25] published the first open-source topology optimization code for the SIMP method. Building upon this, Andreassen et al. [26] performed algorithmic optimizations to further enhance computational efficiency and compactness. Subsequently, Liu et al. [27] introduced the compact Top3 d MATLAB code for tackling 3D topology optimization problems, and Ferrari et al. [28] proposed a new code that further improved optimization efficiency. Zhang et al. [29] proposed a reduced-order modeling method to significantly improve the computational efficiency of topology optimization. Despite these advancements, research focused on enhancing the optimization efficiency of the BESO method remains limited. Therefore, to improve optimization speed, this study proposes an adaptive evolutionary ratio strategy based on the BESO method, along with four corresponding adaptive ER functions.

Several open-source codes developed for nonlinear topology optimization have played pivotal roles in advancing the field. Zhu et al. [30] developed an 89-line geometrically nonlinear topology optimization code based on the Solid Isotropic Material with Penalization (SIMP) method in FreeFEM. Han et al. [31] presented a geometrically nonlinear topology optimization program, based on the BESO method, aimed at minimizing the compliance of statically loaded structures. Subsequently, Zhao et al. [32] proposed two comprehensive MATLAB code sets for 3D geometrically nonlinear topology optimization. Wei et al. [33] incorporated stress constraints in the optimization of porous composite materials, demonstrating that the maximum von Mises stress in the L-bracket structure was reduced by 29% when compared to the case without stress.

However, most topology optimization algorithms rely on MATLAB programs rather than commercial FEA software. This limitation hinders the seamless integration and collaborative development of topology optimization programs with finite element simulation tools. In this paper, Abaqus is utilized for FEA, and the Python program is employed for topology optimization. Building upon Xie’s research on the BESO method [2], this study introduces an adaptive evolutionary ratio strategy to enhance the optimization efficiency of the BESO method. The method is then extended to the design of structures considering both material and geometric nonlinearities. Given that Abaqus is a powerful FEA platform, the proposed topology optimization approach has the potential to be applied to large-scale engineering structures and to address more complex optimization problems.

The remainder of this paper is organized as follows: Section 2 introduces the improved BESO method. Section 3 discusses the validation of the improved BESO algorithm. Section 4 presents nonlinear topology optimization through 3D numerical cases. Finally, Section 5 draws conclusions.

The Improved BESO Method

Nonlinear Analysis

Material Nonlinearity

For most materials, there is a distinct yield stress, denoted as σ_s. If the stress increases with the increase of deformation after reaching σ_s, the material enters a strain-hardening state. Considering a typical multilinear elastic-plastic material model, the stress-strain curve of the material is shown in Figure 1.

[See PDF for image]

Figure 1

The real stress-strain relation of steel material

When considering material nonlinearity, the energy distribution in the elements differs from that in linear elasticity. In material nonlinearity, the energy of each element typically consists of two components: the elastic strain energy $E_e^e$ and the plastic strain energy $E_e^p$ [34], as shown in Figure 2.

[See PDF for image]

Figure 2

Strain energy components: elastic strain energy $E_e^e$ and plastic strain energy $E_e^p$

In Abaqus, the elastic and plastic strain energies of elements can be obtained from the field output. The"SENER"keyword corresponds to the elastic strain energy density, relative to the current volume, while"PENER"represents the energy dissipated by rate-independent and rate-dependent plasticity per unit volume. Specifically, for the unit volume, the "SENER" value represents elastic strain energy, whereas the "PENER" value represents plastic strain energy.

Geometric Nonlinearity

In this study, it is assumed that the structures undergo large deformation with small strain. To model this nonlinear behaviour, the total Lagrangian (TL) formulation is utilized, in which all static and kinematic variables are referred to the initial undeformed configuration of the structure and the integrals are calculated with respect to that configuration. Because of the transformations, the second Piola–Kirchhoff stress tensor, has to be introduced with the Green–Lagrange strain tensor. Considering TL formulation for a general body subjected to physical force p₀ and surface load q₀ on the surface A and displacement field δu, the virtual work equation within the element can be expressed as

1

$$\underset{e_0}{\int }\delta {\mathit{E}}^{\text{T}}\mathit{S}\text{d}V=\underset{e_0}{\int }\delta {\mathit{u}}^{\text{T}}{\mathit{p}}_0\text{d}V+\underset{A_{e0}}{\int }\delta {\mathit{u}}^{\text{T}}{\mathit{q}}_0\text{d}A,$$where the integral domains e₀ and A_e0 are the volume of the element in its initial configuration and the corresponding surface area, respectively.

By integrating the equilibrium equations of all the elements, the equilibrium equation of the whole finite element system can be obtained. If the integral over the whole area of the object V₀ and over the entire boundary A₀ means summation, then the equilibrium equation of the system can be simplified as

2

$$\underset{V_0}{\int }{\mathit{B}}_0^{\text{T}}\mathit{S}\text{d}V=\underset{V_0}{\int }{\mathit{N}}^{\text{T}}{\mathit{p}}_0\text{d}V+\underset{A_0}{\int }{\mathit{N}}^{\text{T}}{\mathit{q}}_0\text{d}A.$$

The incremental Newton–Raphson method is utilized to find the equilibrium solution at every evolutionary iteration. While the response of a linear structure can be determined by simply solving a set of linear equations, the equilibrium of a nonlinear structure needs to be found by an iterative procedure. The residual force, R, is defined as the discrepancy between the internal force vector and the external force vector as

3

$${\mathit{R}\mathbf{=}\mathit{F}\mathbf{-}\mathit{F}}^{\text{int}}=0.$$

The internal force vector can be expressed as

4

$${\mathit{F}}^{\text{int}}=\sum _{e=1}^M\underset{e}{\int }{\mathit{C}}^{\text{eT}}\mathit{B}\sigma \text{d}v=\sum _{e=1}^M{\mathit{C}}^{\text{eT}}{\mathit{F}}_e^{\text{int}},$$where ${\mathit{F}}^{\text{int}}$ is the nodal force vector of an element and ${\mathit{C}}^{\text{eT}}$ is a matrix which transforms the nodal force vector of an element to the global nodal force vector.

Quite often, the equilibrium Eq. (3) is solved incrementally and iteratively using the Newton-Raphson method, which requires the determination of the tangential stiffness matrix in each step as

5

$${\mathit{K}}_{\text{t}}=-\frac{\partial \mathit{R}}{\partial \mathit{U}}.$$

Abaqus uses the Newton-Raphson method to solve nonlinear equilibrium equations, the solution process of which is illustrated in Figure 3. The applied load (F) is divided into smaller increments. For each increment, the displacement is computed using the tangent stiffness matrix (${\mathit{K}}_{\text{t}}$), and the residual force is calculated. The unbalanced force, which is the difference between the applied and the external forces, is determined. The iterative process continues with recalculating the tangent stiffness matrix and updating the displacement and unbalanced force (Figure 3).

[See PDF for image]

Figure 3

The process of solving nonlinear equilibrium equations

Optimization Formulation

The objective of design optimization is to maximize structural stiffness, which is equivalent to minimizing the total external work W depicted as the shaded area in Figure 4. In practice, W is approximated through numerical integration using the trapezoidal rule.

[See PDF for image]

Figure 4

Force-displacement curve in nonlinear finite element analysis

Thus, the optimization problem with a volume constraint can be formulated as

6

$$\begin{array}{c}\text{Minimize:}f(\mathit{x})=W=\underset{n\to \infty }{lim}\left[\frac{1}{2}\sum _{i=1}^n\mathrm{\Delta }{\mathit{U}}^{\text{T}}({\mathit{F}}_i+{\mathit{F}}_{i-1})\right],\\ \text{Subject to:}\left\{\begin{array}{c}{\mathit{R}\mathbf{=}\mathit{F}\mathbf{-}\mathit{F}}^{\text{int}}=0,\\ V(\mathit{x})=\sum _{e=1}^M{v}_e{x}_e=V_{}^{\ast },\\ \mathit{x}=\{x_e\},x_e=1\text{or}{x}_{\text{min}},e=1,\cdots ,M,\\ \end{array}\right)\\ \end{array}$$where n is the total number of displacement increments, $\mathrm{\Delta }\mathit{U}$ is the increment of the nodal displacement vector and ${\mathit{F}}_i$ is the external nodal force vector at the i-th load increment. x is the binary design variables where 1 denotes solid elements and x_min is set to 0.001 to represent void elements. M is the total number of discrete elements. V is the total volume of the structure, v_e is the volume of the element, and V^* is the imposed volume constraint.

Optimization Algorithm

Sensitivity Analysis

In order to perform the topology optimization, the sensitivity of the objective function f with respect to the topology design variables x needs to be provided

7

$$\frac{\partial f}{\partial x_e}=\underset{n\to \infty }{lim}\left[\begin{array}{c}\frac{1}{2}\sum _{i=1}^n({\mathit{U}}_i^{\text{T}}-{\mathit{U}}_{i-1}^{\text{T}})\left(\frac{\partial {\mathit{F}}_i}{\partial x_e}+\frac{\partial {\mathit{F}}_{i-1}}{\partial x_e}\right)+\\ \frac{1}{2}\sum _{i=1}^n\left(\frac{\partial {\mathit{U}}_i^{\text{T}}}{\partial x_e}-\frac{\partial {\mathit{U}}_{i-1}^{\text{T}}}{\partial x_e}\right)({\mathit{F}}_i+{\mathit{F}}_{i-1})\\ \end{array}\right].$$

It is noted that the second term is zero because the variation of displacement equals zero, $\frac{\partial {\mathit{U}}_i^{\text{T}}}{\partial x_e}=0$, $\frac{\partial {\mathit{U}}_{i-1}^{\text{T}}}{\partial x_e}=0$. The sensitivity can be described as

8

$$\frac{\partial f}{\partial x_e}=\underset{n\to \infty }{lim}\left[\frac{1}{2}\sum _{i=1}^n({\mathit{U}}_i^{\text{T}}-{\mathit{U}}_{i-1}^{\text{T}})\left(\frac{\partial {\mathit{F}}_i}{\partial x_e}+\frac{\partial {\mathit{F}}_{i-1}}{\partial x_e}\right)\right].$$

The derivation of the sensitivity requires using the adjoint method. An adjoint equation is introduced by adding a series of vectors of Lagrangian multipliers ${\lambda }_i$ into the objective function as

9

$$f(\mathit{x})=\underset{n\to \infty }{lim}\frac{1}{2}\sum _{i=1}^n\left[\left({\mathit{U}}_i^{\text{T}}-{\mathit{U}}_{i-1}^{\text{T}}\right)\left({\mathit{F}}_i+{\mathit{F}}_{i-1}\right)+{\lambda }_i^{\text{T}}\left({\mathit{R}}_i+{\mathit{R}}_{i-1}\right)\right],$$where ${\mathit{R}}_i$ and ${\mathit{R}}_{i-1}$ are the residual forces in load increment steps i and i − 1. From Eq. (3), the sensitivity of the modified objective function becomes

10

$$\begin{array}{c}\frac{\partial f}{\partial x_e}=\underset{n\to \infty }{lim}\frac{1}{2}\sum _{i=1}^n\left[\begin{array}{c}({\mathit{U}}_i^{\text{T}}-{\mathit{U}}_{i-1}^{\text{T}})\left(\frac{\partial {\mathit{F}}_i}{\partial x_e}+\frac{\partial {\mathit{F}}_{i-1}}{\partial x_e}\right)+\\ {\lambda }_i^{\text{T}}\left(\frac{\partial {\mathit{R}}_i}{\partial {\mathit{U}}_i}\frac{\partial {\mathit{U}}_i}{\partial x_e}+\frac{\partial {\mathit{R}}_{i-1}}{\partial {\mathit{U}}_{i-1}}\frac{\partial {\mathit{U}}_{i-1}}{\partial x_e}+\frac{\partial \left({\mathit{R}}_i+{\mathit{R}}_{i-1}\right)}{\partial x_e}\right)\\ \end{array}\right]\\ \begin{array}{cc}& =\end{array}\underset{n\to \infty }{lim}\frac{1}{2}\sum _{i=1}^n\left[\begin{array}{c}({\mathit{U}}_i^{\text{T}}-{\mathit{U}}_{i-1}^{\text{T}})\left(\frac{\partial {\mathit{F}}_i}{\partial x_e}+\frac{\partial {\mathit{F}}_{i-1}}{\partial x_e}\right)+\\ {\lambda }_i^{\text{T}}\left(\frac{\partial {\mathit{F}}_i}{\partial x_e}-\frac{\partial {\mathit{F}}_i^{\text{int}}}{\partial x_e}+\frac{\partial {\mathit{F}}_{i-1}}{\partial x_e}-\frac{\partial {\mathit{F}}_{i-1}^{\text{int}}}{\partial x_e}\right)\\ \end{array}\right].\\ \end{array}$$

To eliminate the unknowns $\frac{\partial {\mathit{F}}_i}{\partial x_e}+\frac{\partial {\mathit{F}}_{i-1}}{\partial x_e}$, let the Lagrangian multipliers ${\lambda }_i$ be

11

$${\lambda }_i=-\left({\mathit{U}}_i^{\text{T}}-{\mathit{U}}_{i-1}^{\text{T}}\right).$$

By substituting ${\lambda }_i$ into Eq. (10), the sensitivity of the objective function is expressed by

12

$$\frac{\partial f}{\partial x_e}=\underset{n\to \infty }{lim}\frac{1}{2}\sum _{i=1}^n({\mathit{U}}_i^{\text{T}}-{\mathit{U}}_{i-1}^{\text{T}})\left(\frac{\partial {\mathit{F}}_i^{\text{int}}}{\partial x_e}+\frac{\partial {\mathit{F}}_{i-1}^{\text{int}}}{\partial x_e}\right).$$

By introducing the material interpolation scheme, the internal force vector can be expressed by

13

$${\mathit{F}}^{\text{int}}=\sum _{e=1}^M{x}_e^p{\mathit{C}}^{\text{eT}}{\mathit{F}}_{e0}^{\text{int}},$$where ${\mathit{F}}_{e0}^{\text{int}}$ is the internal force vector for solid element. Substituting Eq. (13) into Eq. (12), the variation of the objective function can be expressed as

14

$$\frac{\partial f}{\partial x_e}=\underset{n\to \infty }{lim}\frac{1}{2}p{x}_e^{p-1}\sum _{i=1}^n({\mathit{U}}_i^{\text{T}}-{\mathit{U}}_{i-1}^{\text{T}})\left({\mathit{C}}^{\text{eT}}{\mathit{F}}_{e,i}^{\text{int}}+{\mathit{C}}^{\text{eT}}{\mathit{F}}_{e,i-1}^{\text{int}}\right),$$where ${\mathit{F}}_{e,i}^{\text{int}}$ denotes the internal force vector for solid element in the ith increment. In order to minimize the external work, the elemental sensitivity number is defined, representing the relative ranking of the elemental sensitivity, as

15

$$\begin{array}{c}{\alpha }_e=\frac{1}{p}\frac{\partial f}{\partial x_e}=x_e^{p-1}\underset{n\to \infty }{lim}\frac{1}{2}\sum _{i=1}^n({\mathit{U}}_i^{\text{T}}-{\mathit{U}}_{i-1}^{\text{T}})\left({\mathit{C}}^{\text{eT}}{\mathit{F}}_{e,i}^{\text{int}}+{\mathit{C}}^{\text{eT}}{\mathit{F}}_{e,i-1}^{\text{int}}\right)\\ \begin{array}{cc}& =x_e^{p-1}\underset{n\to \infty }{lim}\sum _{i=1}^n(E_i^e-E_{i-1}^e)\end{array}=x_e^{p-1}{E}_n^e,\\ \end{array}$$where $E_n^e$ is the final elemental elastic and plastic strain energy, which can be extracted using the Abaqus keywords"SENER"and"PENER".

In the BESO method only two discrete design variables $x_{min}$ and $1$ are employed. Therefore, the sensitivity can be simplified as

16

$${\alpha }_e=\left\{\begin{array}{c}{E}_n^e,\begin{array}{ccc}& & x_e=1,\end{array}\\ x_{min}^{p-1}{E}_n^e,\begin{array}{cc}& x_e=x_{min},\end{array}\\ \end{array}\right)$$

Typically, the original sensitivity is processed to produce mesh independent optimization solutions, avoiding checkerboard patterns. For sensitivity filtering, this study adopts the filtering scheme [35]. A simplified element sensitivity filtering scheme is used, described by the following equation:

17

$${\widehat{\alpha }}_e=\frac{{\sum }_{j}w\left(r_{\mathit{ej}}\right){\alpha }_j}{{\sum }_{j}w\left(r_{\mathit{ej}}\right)}=\sum _j\left(\frac{w\left(r_{\mathit{ej}}\right)}{{\sum }_{j}w\left(r_{\mathit{ej}}\right)}{\alpha }_j\right)=\sum _j{\eta }_j{\alpha }_j,$$

18

$$w\left(r_{\mathit{ej}}\right)=max\left(0,r_{min}-r_{\mathit{ej}}\right),$$where $r_{\mathit{ej}}$ is the distance between the centers of element e and j; $w$ is the weight function for averaging the original sensitivity; $r_{min}$ is the filtering radius. The weight factor ${\eta }_j$ is independent of sensitivity values and can be precomputed.

To enhance the stability of the optimization iterative process and achieve a convergent solution, the sensitivity is further averaged with its historical information. This is accomplished by averaging the current iteration's sensitivity with the sensitivity from the previous iteration.

19

$$\stackrel{\sim }{{\alpha }_e}=\frac{{\widehat{\alpha }}_e^k+{\widehat{\alpha }}_e^{k-1}}{2},$$where k is the current iteration number.

Adaptive Evolutionary Ratio Strategy

The evolutionary ratio (ER) is a crucial parameter in the BESO method, representing the ratio of the number of elements updated in each iteration to the total number of elements in the design domain. An appropriate ER can accelerate the optimization process; however, if it is too large, unstable convergence may occur, while if it is too small, optimization speed will be reduced, thus affecting efficiency. The traditional BESO method uses a fixed ER (commonly set to 0.02), which ensures convergence stability but results in relatively slow optimization speed. To reduce the number of iterations and enhance optimization efficiency, a newly designed adaptive ER strategy has been proposed.

For topology optimization problems with volume constraints, the optimization algorithm must quickly eliminate unimportant elements during the early iterations and ensure convergence stability in the later iterations. In the early iterations, many low-sensitivity elements are present in the design domain, which makes it preferable to choose a larger ER to quickly eliminate unimportant elements. In the later iterations, a smaller ER should be selected to minimize element additions and deletions, ensuring that changes in the objective function are small, thereby increasing convergence stability and meeting convergence criteria.

Considering the adaptive volume fraction, the evolutionary ratio varies with the volume fraction (v_h). By utilizing the hyperbolic tangent function and incorporating the characteristics required for the ER during optimization iterations, the following adaptive ER has been designed, as shown in Eq. (20). A more intuitive plot of the ER versus volume fraction is shown in Figure 5.

20

$$ER=\frac{E{R}_{\text{max}}-E{R}_{\text{min}}}{2}\text{tanh}(k(v_h-0.5))+\frac{E{R}_{\text{max}}+E{R}_{\text{min}}}{2},$$where, $E{R}_{\text{max}}$ and $E{R}_{\text{min}}$ represent the maximum evolutionary ratio and minimum evolutionary ratio. $k$ is the control factor, adjusting $k$ can alter the steepness of the curve, $v_h$ represents the volume fraction in the current optimization iteration step.

[See PDF for image]

Figure 5

The adaptive evolutionary ratio function curve

Design Variable Update

BESO typically initiates with a topology filled with solid elements and iteratively reduces the structural volume by transitioning element states. During each iteration, the next target volume $V_{k+1}$ is determined based on the current volume $V_k$ and the evolutionary ratio

21

$$V_{k+1}=V_k\left(1\pm ER\right).$$

The update of elements is based on the optimality criterion of BESO, where the sensitivity of retained solid elements consistently exceeds that of void elements. The update scheme relies on an intermediate threshold, which is determined using a simple binary search method based on the required volume fraction in the current iteration step. Elements with sensitivities exceeding this threshold are converted to solid elements (if they were initially void), while elements with sensitivities below the threshold are switched to void elements (if they were initially solid), ultimately achieving the target volume.

The termination of iterations must satisfy two conditions: the volume constraint and the convergence criterion, which is defined as

22

$$change=\frac{\left|\sum _{i=1}^N(W_{k-i+1}-W_{k-N-i+1})\right|}{\sum _{i=1}^N{W}_{k-i+1}}\le \tau ,$$where k is the current iteration step, τ is a allowable convergence tolerance, set to 0.001 in this study, and N is an integer number, in this study N is set to 5, which adequately captures the average change in the objective function over the past 10 iterations.

Validation of the Improved BESO Algorithm

Topology optimization is implemented using the commercial software Abaqus and a Python program. The nonlinear optimization process, considering the adaptive evolutionary ratio, is as follows:

1. Establish the finite element model in Abaqus, define the design domain, material properties, loads, boundary conditions, finite element mesh, and generate the CAE file;

2. Set the initial parameters of BESO: volume fraction V^*, adaptive evolutionary ratio function ER, and filter radius r_min;

3. Estimate the number of iterations required and set the maximum number of iterations M;

4. Loop:

   1. Perform nonlinear finite element analysis iteratively;

   2. Perform sensitivity calculation according to Eq. (16);

   3. Filter and average sensitivity according to Eqs. (17) and (18);

   4. Update the topological design variables;

   5. Calculate the volume fraction V_K+1 for the current iteration step according to Eq. (21);

   6. Calculate the adaptive ER for the current iteration step.

5. End after satisfying the volume constraints and convergence criteria.

The optimization flowchart is shown in Figure 6.

[See PDF for image]

Figure 6

Flowchart of the adaptive evolutionary ratio nonlinear optimization problem

Mathematical Model for Adaptive Evolutionary Ratio

This study employs a novel design featuring an adaptive ER based on changes in the volume fraction, in contrast to the traditional BESO method with a fixed ER. This section provides a detailed explanation of the design process for the adaptive ER associated with volume fraction and validates the BESO algorithm incorporating this adaptive ER through numerical examples.

Taking into full account the characteristics of the topology optimization process and considering the desired properties of the ER during optimization iterations, a more reasonable relationship between the ER and the volume fraction is proposed. Suitable mathematical functions are utilized to design four adaptive ER functions, denoted as ER-v_h functions, as shown in Eq. (23). Their function curves are shown in Figure 7.

23

$$\begin{array}{c}E{R}_1=\frac{E{R}_{\text{max}}-E{R}_{\text{min}}}{2}\text{tanh}(k_1(v_h-0.5))+\frac{E{R}_{\text{max}}+E{R}_{\text{min}}}{2},\\ E{R}_2=(E{R}_{\text{max}}-E{R}_{\text{min}})\frac{1}{1+exp\left(-k_2(v_h-0.5)\right)}+E{R}_{\text{min}},\\ E{R}_3=(E{R}_{\text{max}}-E{R}_{\text{min}})\frac{log\left(1+exp\left(2k_3(v_h-0.5)\right)\right)}{log\left(1+exp(k_3)\right)}+E{R}_{\text{min}},\\ E{R}_4=E{R}_{\text{min}}+(E{R}_{\text{max}}-E{R}_{\text{min}})v_h,\\ \end{array}$$where $E{R}_i$(i=1,2,3,4) denote the evolutionary ratio; $E{R}_{\text{max}}$,$E{R}_{\text{min}}$ are the maximum and minimum evolutionary ratios, respectively. In this study, $E{R}_{\text{max}}$ is set to 0.05, and $E{R}_{\text{min}}$ is set to 0.02. v_h denotes the volume fraction; $k_i$ (i=1,2,3) are control factors, adjusting the values of $k_i$ can alter the steepness of the curves, generally, a larger k_i results in a steeper curve. In this study, the values of k_i(i=1,2,3) are set to 25, 30, and 6, respectively.

[See PDF for image]

Figure 7

The four adaptive evolutionary ratio curves

Validation of the Effectiveness of the Adaptive Evolutionary Ratio BESO Algorithm

To validate the effectiveness of the four newly proposed adaptive ER functions in the BESO algorithm, these adaptive ER optimization algorithms were applied to the optimization problem of a double clamped beam.

Initially, a finite element model of the double clamped beam was established, and the finite element analysis was performed using the commercial software Abaqus. The dimensions of the clamped beam were 80 × 10 × 4 mm, with material properties defined by Young's modulus E of 210 GPa and Poisson's ratio ν of 0.3. The C3D8R elements with the mesh size of 1 mm were employed, resulting in 3200 mesh elements. Boundary conditions were applied to fix all degrees of freedom at left and right ends, and a concentrated downward force F = 300 N was applied at the center of the lower surface of the cantilever beam, as illustrated in Figure 8.

[See PDF for image]

Figure 8

Design domain, boundary condition and load of double clamped beam case

Subsequently, the optimization process was conducted. The four adaptive ER optimization algorithms were individually applied to the double clamped beam example to test their optimization performance. For the same double-clamped beam optimization example, under identical optimization conditions, with only the ER function parameter varying, the number of iterations was chosen as the evaluation criterion for the optimization algorithm's efficiency. A lower number of iterations indicates higher efficiency. With a volume fraction constraint of 0.3, and a filtering radius of 2, the double clamped beam was optimized. The detail results, including the optimization results, evolution history, total number of iterations, and the optimization efficiency improvement, are shown in Table 1.

Table 1. The optimization result of the four adaptive evolutionary ratio BESO methods

ER-vh
Optimization results
Evolution history
Iteration’s number	47	55	66	61
Optimization efficiency improvement (%)	37.3	26.7	12	18.7

Figure 9 clearly shows the improvement in convergence speed with the BESO method using adaptive ER. The number of optimization iterations shows that the BESO method with adaptive ER significantly reduces the total number of iterations. While the original BESO method with fixed ER requires 75 iterations, the newly proposed four adaptive ER BESO methods require only 47, 55, 66, and 61 iterations, respectively, representing decreases of 37.3%, 26.7%, 12%, and 18.7% in the number of iterations. This substantial enhancement in optimization efficiency is particularly noteworthy, with ER₁ yielding the most significant improvement in efficiency.

[See PDF for image]

Figure 9

Iteration's number using BESO with adaptive ER and fixed ER for double clamped beam

The optimized results obtained from the four different adaptive ER BESO algorithms exhibit very similar topological configurations, with only minor differences. These results can be considered nearly identical. This aligns with the expected outcome, as different ERs theoretically only alter the number of optimization iterations, influencing optimization efficiency, without affecting the optimization results.

To further illustrate the improvement in optimization efficiency with the adaptive ER BESO algorithm, a comparison was made between the structural evolution histories of the adaptive ER BESO algorithm and the original fixed ER BESO algorithm. Compliance and topology structures were extracted from both algorithms at equivalent iteration steps, as depicted in Figure 10.

[See PDF for image]

Figure 10

Comparison of the convergence process between the adaptive ER algorithm and the original fixed ER algorithm

It can be observed that the optimization process with the adaptive ER algorithm is completed in only 47 iterations, while the original fixed ER optimization algorithm requires 75 iterations. The number of optimization iterations is decreased by 37.3%, demonstrating a significant improvement in efficiency. The final optimization results are nearly identical, with objective function values of 11.039 and 12.056, indicating little errors and validating the accuracy of the adaptive ER optimization algorithm. A comparison of the topology structure at the 15 th and 30 th iterations for the two different ER algorithms, it is evident that the adaptive ER optimization algorithm evolves more rapidly. At the same iteration steps, the topology structure is more similar to the final optimization result. The adaptive ER algorithm reaches the optimal topology more quickly. This is attributed to the larger ER in the early stages of optimization, allowing the algorithm to quickly remove unimportant structural elements.

It can be concluded that the four adaptive ER optimization algorithms exhibit a significant improvement in convergence speed, resulting in a substantial enhancement in optimization efficiency.

To further validate the effectiveness of the BESO method with adaptive ER, the optimization results were compared with those from the well-known Solid Isotropic Material with Penalization (SIMP) method [28] and the Moving Morphable Components (MMC) method [36].

Figure 11 presents the widely recognized cantilever beam with characteristic dimensionless dimensions. The design domain has dimensions of 80 × 40, with the degrees of freedom constrained at the left end and a vertical load F = 1 applied at the central point of the right-hand end. The volume fraction constraint is 0.5. The design incorporates a filtering radius of r_min = 4 and a penalty factor of p = 3.

[See PDF for image]

Figure 11

Design domain, boundary condition and load of cantilever beam case.

Figure 12 shows the design results of different topology optimization methods. Figure 11a shows the results using the BESO method with adaptive ER, while Figure 11b–d shows the optimization results using the BESO method with fixed ER, the Solid Isotropic Material with Penalization (SIMP) method and the Moving Morphable Components (MMC) method, respectively. C represents the objective function, which is compliance. The optimization results across these methods are largely consistent, with the compliance (C) error being extremely small, thereby validating the effectiveness of the BESO method with adaptive ER.

[See PDF for image]

Figure 12

Design results of different topology optimization methods: a BESO method with adaptive ER, b BESO method with fixed ER, c Solid Isotropic Material with Penalization (SIMP) method, d Moving Morphable Components (MMC) method

Nonlinear Topology Optimization

In this section, two 3D cases are employed to conduct materially and geometrically nonlinear topology optimization analysis using the BESO method with adaptive ER. The ER₁ algorithm demonstrates the fastest optimization speed. Therefore, in subsequent nonlinear optimization studies, the optimization algorithm corresponding to the ER₁ adaptive ER will be adopted.

Double Clamped Beam Case

The finite element model is illustrated in Figure 8. It is worth noting that for material nonlinearity, an elastic-plastic material is defined in this study, as shown in Figure 1. For geometric nonlinearity, a static analysis step with geometric nonlinearity option must be defined in Abaqus.

The improved BESO method is employed to investigate scenarios involving material nonlinearity, geometric nonlinearity, and combined nonlinearity (simultaneous consideration of material nonlinearity and geometric nonlinearity). To facilitate result comparison, optimization is also conducted under linear elastic conditions. Various concentrated loads are considered, and the optimization results are summarized in Table 2.

Table 2. Optimization results of double clamped beam under linear and various nonlinearity

Load (N)	Linear		Material nonlinearity	Geometric nonlinearity	Combined nonlinearity
300
500
750

From the topology optimization results in Table 2, it can be observed that the results of nonlinear topology optimization differ from linear optimization results, especially when the external load is large, indicating a noticeable influence of nonlinearity on the optimization results.

Initially, individual analyses were conducted to explore material nonlinearity, geometric nonlinearity, and combined nonlinearity behaviors. Under linear conditions, optimization results remained consistent across varying external loads, in line with theoretical expectations. However, nonlinear optimization scenarios exhibited evident variations in results as external loads increased, suggesting that the degree of nonlinearity influences the optimal topology of the structure. Notably, significant changes were observed in materially nonlinear optimization results at an external load of 750 N, characterized by structural inversion at both ends. Geometric nonlinearity resulted in more pronounced alterations in topology structures. At 500 N, the structure removed the middle truss section and increased connections to both ends, while 750 N, the clamped beam topology transformed into a V-shape. Combined nonlinearity produced topology structures combining material and geometric nonlinear effects, evident in a V-shaped structure under higher loads.

Under the same external loading conditions, comparative analyses of optimization results under various nonlinearities were conducted. At relatively low loads (F=300 N), the optimization results for linear, material nonlinearity, geometric nonlinearity, and combined nonlinearity cases were similar, forming symmetric topology structures resembling triangular trusses. This similarity is attributed to the minimal nonlinear effects at low loads. As the loads increased to F=500 N, the nonlinear influences became more pronounced, particularly in the case of geometric nonlinearity, leading to structural evolution toward a V-shape. Similar changes were observed in the combined nonlinear cases. At F=750 N, optimization results for all nonlinear conditions varied, with material nonlinearity causing topology inversions at both ends, while geometric nonlinearity and combined nonlinearity resulted in V-shaped structures.

Additionally, the evolution history of the clamped beam under F=500 N and combined nonlinearity conditions was presented in Figure 13. It can be observed that convergence was achieved after 46 iterations. Notably, when employing the BESO method with adaptive ER, the total number of iterations was remarkably low, significantly enhancing optimization efficiency. The objective function (W) history curve reveals that it gradually increases throughout the iterations, reaching its maximum value at the 37 th iteration, followed by slight oscillations before stabilizing until convergence. The volume fraction history curve displays a nonlinear pattern, distinct from previous studies, with a rapid decline in the early iterations and a slowed descent in the later stages. This behavior is attributed to the use of the proposed adaptive ER, which demonstrates faster optimization speed in the initial stages and improved convergence stability in the later stages of optimization.

[See PDF for image]

Figure 13

The evolution history of the clamped beam under F = 500 N and combined nonlinear conditions

L-Shaped Beam Case

The L-shaped beam case is shown in Figure 14 with dimensions in millimeters. The material properties include Young's modulus E=210 GPa, and Poisson's ratio ν=0.3, utilizing the multilinear elastic-plastic model, as shown in Figure 1. The C3D8R elements with mesh size of 1 mm were used. All degrees of freedom at the top end of the structure are fixed. A load F=100 N is applied at the top of the right face, distributed over 9 nodes to avoid stress concentration. The static analysis step with geometric nonlinearity option must be defined in Abaqus.

[See PDF for image]

Figure 14

Design domain, boundary condition and load of L-shaped beam case

A volume fraction of 0.3 is set, with a filtering radius of 2. To address numerical instability during the convergence of nonlinear optimization, the adaptive evolutionary ratio ER₁ (Eq.(23)) function is adjusted, where ER_max =0.05, ER_min=0.01, ensuring a smaller evolutionary ratio for stable optimization results.

Both the adaptive ER BESO method and the fixed ER BESO method were applied to perform combined nonlinear topology optimization on the L-shaped beam, resulting in optimization results and evolution histories. The comparison of results using BESO with adaptive ER and fixed ER for the L-shaped beam are shown in Figure 15.

[See PDF for image]

Figure 15

Comparison of results using BESO with adaptive ER and fixed ER for L-shaped beam

It can be observed that the optimization process with the adaptive ER algorithm is completed in only 69 iterations, while the fixed ER algorithm requires 145 iterations. The convergence speed is increased by 52.4%, indicating a significant improvement in efficiency. The final optimization results are almost identical, with objective function values of 118.7 and 117.4, with very small errors, which verifies the accuracy of the adaptive ER algorithm. Comparing the topology structures at the 20 th and 50 th iterations for the two different ER algorithms, it is evident that the adaptive ER algorithm evolves more rapidly. At the same iteration steps, the topology structure is closer to the final optimization result. The adaptive ER algorithm can reach the optimal topology more quickly. This is attributed to the larger ER in the early stages of optimization, which allows the algorithm to remove unimportant structural elements more efficiently. In the late stages of optimization, a smaller ER is maintained to ensure convergence stability.

Conclusions

This study extends volume-constrained maximum stiffness design based on the BESO method to structures exhibiting material and geometric nonlinearities. To address the slow convergence issue inherent in traditional BESO methods, an innovative strategy is proposed, which introduces an adaptive evolutionary ratio to enhance the efficiency of the BESO method. The optimization results from several benchmark cases validate the effectiveness and practicality of the proposed method. Since Abaqus is a powerful FEA platform, the proposed topology optimization approach has the potential to be applied to large-scale engineering structures and to address more complex optimization problems. The main conclusions are summarized as follows.

1. The four proposed adaptive evolutionary ratio BESO methods exhibit significantly improved convergence speeds compared to the original fixed evolutionary ratio BESO method, with enhancements of 37.3%, 26.7%, 12%, and 18.7%, respectively. Notably, excellent evolution speed is demonstrated in the early iterations.

2. The BESO method with adaptive ER shows efficient optimization speed in the early iterations and maintains convergence stability in the later iterations.

3. Optimization results differ under material nonlinearity, geometric nonlinearity, and combined nonlinearity. Combined nonlinearity can be understood as the superposition of material and geometric nonlinear effects.

4. Varying degrees of nonlinearity lead to different optimal topological structures. Therefore, in engineering practice, the design of the most suitable topology structure should consider the degree of nonlinearity.

Acknowledgements

Not applicable.

Author Contributions

Wenhua Zhang was in charge of the modeling, data analysis and wrote the manuscript; Linli Tian provided fundamental ideas and all support conditions of this paper, reviewed and edited the manuscript. All authors read and approved the final manuscript.

Funding

Supported by National Natural Science Foundation of China (Grant No. 52105271).

Data availability

The data supporting this study are available from the corresponding author upon reasonable request.

Declarations

Competing Interests

The authors declare no competing financial interests.

References

[1] Bendsøe, MP; Kikuchi, N. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering; 1988; 71, 2 pp. 197-224.969217

[2] Zuo, ZH; Xie, YM. A simple and compact Python code for complex 3D topology optimization. Advances in Engineering Software; 2015; 85, pp. 1-11.

[3] Xu, B; Han, Y; Zhao, L. Bi-directional evolutionary topology optimization of geometrically nonlinear continuum structures with stress constraints. Applied Mathematical Modelling; 2020; 80, pp. 771-791.4046538

[4] Han, Y; Xu, B; Wang, Q et al. Topology optimization of material nonlinear continuum structures under stress constraints. Computer Methods in Applied Mechanics and Engineering; 2021; 378, 4225684 113731.

[5] T T Banh, D Lee. Efficient topology optimization for geometrically nonlinear multi-material systems under design-dependent pressure loading. Engineering with Computers, 2024: 1–35.

[6] Wang, B; Bai, J; Lu, S et al. Structural topology optimization considering geometrical and load nonlinearities. Computers & Structures; 2023; 289, 107190.

[7] Shobeiri, V. Bidirectional evolutionary structural optimization for nonlinear structures under dynamic loads. International Journal for Numerical Methods in Engineering; 2019; 121, 5 pp. 888-903.4072522

[8] Chen, X; Yao, S; Yvonnet, J. Nonlinear topology optimization of flexoelectric soft dielectrics at large deformation. Computer Methods in Applied Mechanics and Engineering; 2024; 427, 4739762 117005.

[9] da Silva, GA; Beck, AT; Sigmund, O. Topology optimization of compliant mechanisms considering stress constraints, manufacturing uncertainty and geometric nonlinearity. Computer Methods in Applied Mechanics and Engineering; 2020; 365, 4078502 112972.

[10] Ye, H; Yuan, B; Li, J et al. Geometrically nonlinear topology optimization of continuum structures based on an independent continuous mapping method. Acta Mechanica Solida Sinica; 2021; 34, pp. 658-672.

[11] Chen, Z; Wen, G; Wang, H et al. Multi-resolution nonlinear topology optimization with enhanced computational efficiency and convergence. Acta Mechanica Sinica; 2022; 38, 24422470 421299.

[12] Scherz, L; Kriegesmann, B; Pedersen, CB. A condition number-based numerical stabilization method for geometrically nonlinear topology optimization. International Journal for Numerical Methods in Engineering; 2024; 125, 23 e7574.

[13] Liang, J; Zhang, X; Zhu, B et al. Topology optimization of geometrically nonlinear structures based on a self-adaptive material interpolation scheme. Machines; 2023; 11, 12 1047.

[14] Xu, T; Huang, X; Lin, X et al. Topology optimization for maximizing buckling strength using a linear material model. Computer Methods in Applied Mechanics and Engineering; 2023; 417, 4644275 116437.

[15] Zhang, W; Jiu, L; Meng, L. Buckling-constrained topology optimization using feature-driven optimization method. Structural and Multidisciplinary Optimization; 2022; 65, 1 37.4365191

[16] Zhang, W; Yang, W; Zhou, J et al. Structural Topology Optimization Through Explicit Boundary Evolution. Journal of Applied Mechanics; 2017; 84, 1 011011.

[17] Xue, R; Liu, C; Zhang, W et al. Explicit structural topology optimization under finite deformation via Moving Morphable Void (MMV) approach. Computer Methods in Applied Mechanics and Engineering; 2019; 344, pp. 798-818.3876233

[18] Guo, Y; Du, Z; Liu, C et al. Explicit topology optimization of three-dimensional geometrically nonlinear structures. Acta Mechanica Sinica; 2023; 39, 124642261 423084.

[19] Du, Z; Cui, T; Liu, C et al. An efficient and easy-to-extend Matlab code of the Moving Morphable Component (MMC) method for three-dimensional topology optimization. Structural and Multidisciplinary Optimization; 2022; 65, 5 pp. 1-29.

[20] Behrou, R; Lotfi, R; Carstensen, JV et al. Revisiting element removal for density-based structural topology optimization with reintroduction by Heaviside projection. Computer Methods in Applied Mechanics and Engineering; 2021; 380, 4238944 113799.

[21] Giele, R; Groen, J; Aage, N et al. On approaches for avoiding low-stiffness regions in variable thickness sheet and homogenization-based topology optimization. Structural and Multidisciplinary Optimization; 2021; 64, 1 pp. 39-52.4274755

[22] Zhang, G; Khandelwal, K; Guo, T. Finite strain topology optimization with nonlinear stability constraints. Computer Methods in Applied Mechanics and Engineering; 2023; 413, 4593965 116119.

[23] Zhang, XS; Chi, H; Zhao, Z. Topology optimization of hyperelastic structures with anisotropic fiber reinforcement under large deformations. Computer Methods in Applied Mechanics and Engineering; 2021; 378, 4222255 113496.

[24] Kumar, P; Schmidleithner, C; Larsen, N et al. Topology optimization and 3D printing of large deformation compliant mechanisms for straining biological tissues. Structural and Multidisciplinary Optimization; 2021; 63, pp. 1351-1366.4233685

[25] Sigmund, O. A 99 line topology optimization code written in Matlab. Structural and Multidisciplinary Optimization; 2001; 21, 2 pp. 120-127.

[26] Andreassen, E; Clausen, A; Schevenels, M et al. Efficient topology optimization in MATLAB using 88 lines of code. Structural and Multidisciplinary Optimization; 2010; 43, 1 pp. 1-16.

[27] Liu, K; Tovar, A. An efficient 3D topology optimization code written in Matlab. Structural and Multidisciplinary Optimization; 2014; 50, 6 pp. 1175-1196.3312425

[28] Ferrari, F; Sigmund, O. A new generation 99 line Matlab code for compliance topology optimization and its extension to 3D. Structural and Multidisciplinary Optimization; 2020; 62, 4 pp. 2211-2228.4156382

[29] Zhang, L; Zhang, Y; Van Keulen, F. Topology optimization of geometrically nonlinear structures using reduced-order modeling. Computer Methods in Applied Mechanics and Engineering; 2023; 416, 4632647 116371.

[30] Zhu, B; Zhang, X; Li, H et al. An 89-line code for geometrically nonlinear topology optimization written in FreeFEM. Structural and Multidisciplinary Optimization; 2020; 63, 2 pp. 1015-1027.4214535

[31] Han, Y; Xu, B; Liu, Y. An efficient 137-line MATLAB code for geometrically nonlinear topology optimization using bi-directional evolutionary structural optimization method. Structural and Multidisciplinary Optimization; 2021; 63, 5 pp. 2571-2588.4244972

[32] Zhao, Y; Guo, G; Zuo, W. MATLAB implementations for 3D geometrically nonlinear topology optimization: 230-line code for SIMP method and 280-line code for MMB method. Structural and Multidisciplinary Optimization; 2023; 66, 7 pp. 1-20.4604411

[33] Wei, G; Chen, Y; Li, Q et al. Multiscale topology optimisation for porous composite structures with stress-constraint and clustered microstructures. Computer Methods in Applied Mechanics and Engineering; 2023; 416, 4627709 116329.

[34] Jia, J; Da, D; Hu, J et al. Crashworthiness design of periodic cellular structures using topology optimization. Composite Structures; 2021; 271, 114164.

[35] Su, R; Tangaramvong, S; Song, C. Automatic image-based SBFE-BESO approach for topology structural optimization. International Journal of Mechanical Sciences; 2024; 263, 108773.

[36] Li, Z; Xu, H; Zhang, S. Moving morphable component (MMC) topology optimization with different void structure scaling factors. Plos one; 2024; 19, 1 e0296337.