A parameterized level set method for structural topology optimization based on reaction diffusion equation and fuzzy PID control algorithm

1.   Introduction

Structural topology optimization can achieve the pursuit of saving materials or improving the performance of structures with the given design domain and the relevant boundary conditions. Hence it has already been a powerful tool for the design of engineering structure. As an important branch of structural topology optimization, continuum topology optimization has attracted more and more attention. With the development of computational mechanics, there has been significant progress made in the theory and methodology of continuum topology optimization over the past 30 years. During this period, much pioneering research has been published, advancing related theories and techniques, and spurring the use of practical topology optimization methods. Accordingly, a variety of methods and their accompanying technical schemes have been proposed to deal with continuum topology optimization and obtain effective solutions. The methods for continuum topology optimization can be commonly classified into the following categories: the homogenization method (Bendsøe and Kikuchi ^[29]; Bendsøe ^[31]), the SIMP method (Bendsøe and Sigmund ^[30]), the evolutionary structural optimization (ESO) method (Xie and Steven ^[55]), the phase field method (Bourdin and Chambolle ^[1]; Marino et al. ^[27]), the level set method (Allaire et al. ^[7]; Wang et al. ^[35]), the feature driven optimization approach, alternatively named as the moving morphable voids approach (Guo et al. ^[53]; Zhou et al. ^[56]; Zhang et al. ^[51]; Zhang et al. ^[52]; Zhu et al. ^[20]), and so on.

The level set method was first proposed by Osher and Sethian ^[45] to simulate moving boundaries. Initially, the level set method was mainly used to simulate the evolution of interface in multiphase flow (Sussman et al. ^[32]; Sethian and Smereka ^[16]) and image processing technology (Malladi et al. ^[43]; Osher and Paragios ^[46]). Sethian and Wiegmann ^[15] used the immersion interface method to calculate the structural stress, and then used the level set method to change the shape of the structure. Osher and Santosa ^[44] employed the level set method to solve the structural topology optimization problem of the vibration system.

Virtual time $t$ is introduced into the traditional level set method to make it more effective for topology optimization. By calculating the derivative of level set function relative to virtual time $t$ (van Dijk et al. ^[38]), the Hamilton-Jacobi partial differential equation (Peng et al. ^[3]) (PDE) related to virtual time $t$ is established. The upwind method can be employed to solve PDE, thereby updating the level set equation and realizing the structural topology optimization. However, when the upwind method is used to solve PDE for topology optimization, the computational efficiency is not high because the time step is limited by the Courant-Friedrichs-Lewy (CFL) condition. Moreover, the holes in the structure cannot be generated, so the evolution of the structure boundary depends on the initial conditions.

To improve the traditional level set method, many researches have been carried out. Xia et al. ^[42] proposed a semi-Lagrangian method to update the level set function to solve the problem that the time step is limited by the CFL condition when using the upwind method, and compared the method with the upwind method to verify the effectiveness of the method. Zhou and Wang ^[37] used the line search algorithm to determine the time step adaptively, which improves the oscillation problem of the objective function and volume constrained optimization curve when using the semi-Lagrangian method. Eschenauer et al. ^[8] proposed a "bubble" method, which inserts a "bubble" into the structure, and uses the positioning criteria and shape optimization method to determine the optimal position and shape of the "bubble". This "bubble" method for structural topology optimization can generate holes in the structure. Sokolowski and Zochowski ^[19] used topological derivative to calculate the change of objective function when adding or deleting materials in the structure, and gave the topological derivative form of plane elastic system. Burger et al. ^[24] combined the topological derivative with the level set method to calculate the sensitivity of the objective function when the structure topology changes.

To further improve the level set method, Allaire et al. ^[6] used the adjoint method to calculate the shape derivative, and took the shape derivative as the normal velocity of the free boundary in the optimization process. At the same time, in order to control that the level set function near the boundary is neither too steep nor too flat, the operation of reinitializing the level set function periodically was introduced to ensure that the structure can converge to the best shape. Luo et al. ^[57] applied the optimization criterion method (Zhou and Rozvany ^[36]) to the parametric level set method to replace the steepest descent method to update the design variables. Allaire et al. ^[5] combined the topological derivative with shape sensitivity, thereby establishing the algorithm model of alternating descent of shape sensitivity and topological derivative, which deals with the boundary evolution problem and meanwhile can introduce holes in the structure. Wei and Wang ^[40] combined the level set method with the phase field method, and implemented the piecewise constant level set method by using the numerical scheme of the additive operator splitting in the phase field method, which describes the evolution of the structural boundary. Yamada et al. ^[49] proposed a new method to construct the reaction diffusion equation (RDE) by using the virtual energy term in the phase field method, and gave the derivation process of RDE. Choi et al. ^[18] took the design sensitivity of topology optimization as the reaction term of RDE, which can generate new holes in the structure without using the topological derivative. Otomori et al. ^[28] provided a MATLAB implementation of level set method based on RDE, and gave the derivation process of topological derivative.

The method of topology description function (TDF) was proposed by de Ruiter and van Keulen ^[26], which utilizes a set of basis function to describe the geometric structure. Through adjusting the shape, parameters and distribution of the basis function, the structure can be described reasonably by the method of TDF. Thereafter, the TDF method has evolved into the parameterized level set method (Cecil et al. ^[48]; Wang and Wang ^[47]; Wei et al. ^[41]; Cui et al. ^[33]) (PLSM). The PLSM retains the advantage of implicit representation of evolution surface by level set function, thereby transforming topology optimization problem into parameter optimization problem, and meanwhile does not need to solve complex Hamilton-Jacobi PDE directly or carry on the re-initialization operation. Thus good topology optimization effect can be achieved by using the PLSM. Moreover, radial basis function (RBF) and level set equation are used to represent the shape of the structure implicitly in the PLSM. Therefore, the selection of the basis function determines the shape type and detail of the structure boundary described by the level set equation.

According to the supporting size, there are two types of basic functions: compactly supported radial basis function (CS-RBF) and global supported radial basis function (GS-RBF). Wang and Wang ^[47] used the multi-quadric (MQ) function to establish the implicit RBF model interpolating the level set function, transformed Hamilton-Jacobi PDE into ordinary differential equation (ODE), and solved it by finite element method. Luo et al. ^[58] proposed an element-free Galerkin level set method, where CS-RBF is used to construct the element-free shape function and the structural boundary is updated to achieve structural optimization by solving the discrete level set function in time. Ho et al. ^[12] proposed a PLSM based on mobile node by using inverse multi-quadric (IMQ) function, which takes the position of RBF node as the design variable and uses the steepest descent method to push the node to a new position. Wei et al. ^[41] provided a MATLAB code implementation of level set method based on MQ function, and extended it for different boundary conditions and different types of RBF. Cui et al. ^[33] proposed a PLSM combined with CS-RBF and the method of moving asymptotes (Svanberg ^[22]) (MMA) to update the design variables, and introduced shape sensitivity constraint factor in the iteration process of MMA algorithm, which makes the iteration step controllable and speeds up the optimization process.

Three-dimensional (3D) structural topology optimization (Yang et al. ^[54]; Sigmund and Clausen ^[39]) can determine the geometric configuration of engineering structure and improve structural performance for engineering application. Some methods for two-dimensional (2D) structural topology optimization have thus far been successfully extended to 3D problem, which broadens the application scope of structural topology optimization. Du and Olhoff ^[17] studied the topology optimization of 3D structure under design related loads, and proposed an algorithm for generating loading surface of 3D structure based on iso-parametric surface theory. Abolbashari and Keshavarzmanesh ^[25] studied the influence of removal rate, evolution rate and element size on the optimization of 2D and 3D continuum structures by the ESO method. Zhang et al. ^[14] proposed a boundary search algorithm, which can effectively identify the load surface and deal with 3D topology optimization problem under design related loads. Muñoz et al. ^[2] used Cartesian mesh finite element method instead of standard finite element method in topology optimization of 3D structure, and combined h-adaptive mesh optimization method in optimization algorithm to improve the accuracy of boundary definition. Liu and Tovar ^[21] extended the SIMP method to the 3D field, and gave 169 lines of MATLAB code including sensitivity analysis, density filtering and finite element analysis. On this basis, Ferrari and Sigmund ^[4] improved the above-mentioned 169 lines of MATLAB code through more effective assembly of stiffness matrix to accelerate the realization of filter function, and obtained more concise and effective 3D MATLAB code. Zegard and Paulino ^[50] proposed a 3D topology optimization algorithm based on ground structure, which can well deal with the structural optimization problems of concave and internal holes, and gave the corresponding MATLAB code. Wang et al. ^[34] combined the B-spline function and the level set function, divided the whole design domain into a group of 3D subdomain cells, and used B-spline function to update design variables recursively, so as to realize the 3D topology optimization of cellular structure with spatial changes. Recently, Li et al. ^[9,10,11] have developed the RDE method and extended its applications in large-scale 3D topology optimization of multi-physics problems, i.e., mean compliance problem, fluid-structure interaction problem, and thermal fluid problem, by using adaptive mesh refinement or body-fitted mesh adaption.

When the traditional level set method is used for structural topology optimization, the upwind method is often used to solve Hamilton-Jacobi PDE, and the conventional level set function is updated according to the solution of the PDE. In the optimization iterations, the time step is required to meet the CFL condition. Thus the size of the time step is limited, which reduces the computational efficiency. Moreover, the re-initialization operation needs to be performed in each of the optimization iterations to control the gradient of the level set function near the boundary, so the stability of numerical calculation is affected. However, when RDE is used to replace Hamilton-Jacobi PDE update the level set function in structural topology optimization, it allows not only shape but also topological changes, and allows the geometrical complexity of the optimal configuration to be specified. Therefore, the RDE-based level set method has significant advantages over the conventional level set method for structural topology optimization. In this work, a PLSM based on the RDE and the fuzzy adaptive PID control algorithm is proposed to address structural topology optimization of minimizing compliance with volume constraints. The RDE is used as the evolution equation of level set function, and the regularization parameters in RDE are used to adjust and control the geometry complexity of the optimized structure. The topological derivative is not only used to calculate the sensitivity of the objective function when the hole is inserted at any position of the structure, but also used as the reaction term of the RDE to drive the evolution of the level set function. In the optimization iteration process, CS-RBF is used to interpolate the level set function to ensure the smoothness of the approximate surface. Furthermore, the fuzzy PID control algorithm is employed to address the volume constraint, so as to makes the convergence process of structure volume more stable. The proposed method is also applied to 3D structural topology optimization, so as to confirm that this method can enhance the optimization efficiency and the stability of the optimization process on the premise of ensuring the accuracy of the optimization results. The proposed method retains the merits of the PLSM, the RDE and the fuzzy adaptive PID control algorithm. Therefore, it can obtain the reasonable results with relatively high speed and good stability for structural topology optimization

The rest of this paper is organized as follows. In Section 2, the level set method and the corresponding solution process which combines the compactly supported radial basis function and fuzzy PID control algorithm are introduced. In Section 3, the mathematical model of the topology optimization problem which takes the minimum compliance as the objective function and the volume fraction as the constraint condition is established, and the role of the initial PID control parameters is studied. In Section 4, 2D and 3D structural optimization examples are used to verify the feasibility and effectiveness of the proposed method. Finally, the conclusions are drawn in Section 5.

2.   The proposed method

2.1.   Level set method based on RDE

The level set function uses its zero-value curve to describe the structure boundary. By solving the value of level set function at each time, the evolution equation of the structure boundary is obtained, and the shape and topology of the structure in the material domain are clearly described. By taking a 2D structure as an example, the level set function is used to represent the material distribution of the structure as shown in Figure 1.

The relationship between level set function $\phi$ and material distribution in design domain is expressed by the mathematical formula as follows.

where, $D$ represents the design domain of the structure, $\Omega$ represents the solid part of the structure, $\Gamma$ represents the boundary of the structure, and $x$ represents element node after the finite element mesh is divided. $\phi \left( x \right)\;{\text{ > }}\;{\text{0}}$ represents the material domain, $\phi \left( x \right) = 0$ represents the structural boundary, and $\phi \left( x \right)\;{\text{ < 0}}$ represents the void.

A topology optimization problem with minimum compliance as objective function and volume fraction as constraint condition is considered. As shown in Figure 2, structural topology optimization problem generally involves prescribed design domain and boundary conditions. Design domain $D$ usually includes material domain $\Omega$ and void region. With respect to boundary conditions, displacement constraints are usually imposed at the boundary ${\Gamma _u}$, external loads $F$ are often applied at the boundary ${\Gamma _t}$, and body force ${\bf{b}}$ is applied on the overall material domain. Under the volume constraint, the mathematical expression of the objective function is written as follows.

where, $J\left( \Omega \right)$ is the objective function, $j\left( x \right)$ is the function in the structural design domain, $G\left( \Omega \right)$ is the volume constraint function, ${V_{max}}$ is the set volume constraint value, and $x$ represents element node after the structural design domain is divided into finite element meshes.

To introduce the virtual energy term of phase field model into the objective function and regularize the structural topology optimization problem, it needs to be assumed that the level set function $\phi$ and the phase field variables in the phase field method have the same distribution properties. Based on this assumption, the upper and lower bound constraints are imposed on the level set function $\phi$, which stipulates that the value of level set function must be equal to 1 or −1 in the area far enough from the structure boundary, so that the smoothing effect only works on the points near the structure boundary. Therefore, the definition of level set function expressed in Eq (1) is modified, and the revised definition of level set function is shown in Eq (3).

By introducing the virtual energy term in the phase field method into the objective function, the topology optimization problem described by Eq (2) is transformed into the topology optimization problem shown in Eq (4).

where, $\int_D {\frac{1}{2}\tau {{\left| {\nabla \phi } \right|}^2}d\Omega }$ is the added virtual energy term, $\tau$ is the regularization parameter, which is always greater than zero, ${\chi _\Omega }$ is characteristic function, which can transform the structural topology optimization problem into the material allocation problem in the fixed design domain, so as to determine the optimal distribution of the structure in the design domain.

The Lagrange method is used to deal with the topology optimization problem with the given volume constraints, and the constrained optimization problem of Eq (4) can be transformed into an unconstrained optimization problem, which is expressed as follows.

where, $\lambda$ denotes the Lagrange multiplier, which is always positive and used to deal with the volume constraints, and ${\bar J_R}$ is the Lagrange function. The optimal topology of the structure can be obtained by solving Eq (5).

However, it is almost impossible to find the optimal solution of Eq (5) directly. To obtain the final topology optimization results, Equation (5) is solved by using the time evolution equation. The time evolution equation is expressed as follows.

where, $\bar J = \int_D {\left( {j\left( x \right) + \lambda } \right)} {\chi _\phi }\left( \phi \right)d\Omega - \lambda {V_{max}}$. Equation (6) is also called the reaction-diffusion equation (RDE).

To solve the problem conveniently, on the basis of Eq (6), the normalization coefficient $C$ is introduced to normalize the reaction terms in RDE. Since the reaction term $\frac{{d\bar J}}{{d\phi }}$ is equivalent to the topological derivative ${d_t}\bar J$ numerically, the topological derivative ${d_t}\bar J$ is used for numerical calculation. The level set function is independent of the outside of the fixed design domain, so it is assumed that the boundary conditions of the level set function are Dirichlet boundary conditions on the non-design boundary and Neumann boundary conditions on other boundaries. Therefore, the RDE with boundary conditions can be adjusted to the formulation of Eq (7).

where, the value of ${K_0}$ in Eq (6) is set as 1, referring to the literature (Otomori et al. ^[28]). The sign of topological derivative ${d_t}\bar J$ is opposite to that of $\frac{{d\bar J}}{{d\phi }}$, because the sign definitions of topological derivative ${d_t}\bar J$ and level set function $\phi$ are different. The solving process of topological derivative will be detailed in Section 2.4. The formula for calculating the normalization coefficient $C$ is expressed as follows.

2.2.   Level set method based on RDE combined with CS-RBF

Through the above derivation, solving Eq (7) can update the level set function and realize the topology optimization of the structure. When the level set method based on RDE is used for structural optimization, the finite element shape function is used to interpolate the level set equation. In order to move the structure boundary to an ideal position, a large number of iterations are needed in the actual operation.

To obtain better interpolation effect, a continuous and differentiable CS-RBF model is introduced into the structural topology optimization. CS-RBF is used to interpolate the level set equation instead of the finite element shape function to improve the solving process of RDE, which can not only retain the advantages of the level set method based on RDE, but also improve the computational efficiency.

The second-order continuous CS-RBF proposed by Wendland ^[13] is selected as the interpolation function to parameterize the level set equation. Taking the 2D structure as an example, the mathematical expression of CS-RBF is given as follows.

The support radius $r$ is defined as

where, ${d_i}$ represents the distance from any point $\left( {x, y} \right)$ within the support radius $r$ to the fixed reference point $\left( {x_i^c, y_i^c} \right)$. The search radius of basis function ${d_m}$ represents the maximum control range of CS-RBF centered on reference point $\left( {x_i^c, y_i^c} \right)$. The selection of ${d_m}$ is very important. If the value of ${d_m}$ is too large, the number of nodes in the interpolation radius will increase and the amount of calculation in the interpolation process will increase. If the value of ${d_m}$ is too small, the interpolation matrix will be singular. In other words, the value of ${d_m}$ can lead to mesh-dependent optimal solutions. To ensure the non-singularity of the function interpolation and simultaneously guarantee the computational efficiency, an experiential criterion in selecting a radius of support is often needed (Luo et al. ^[57]). In this paper, the value of ${d_m}$ is selected as 2 times of the finite element mesh size to simplify the calculation.

By using CS-RBF to parameterize the level set equation, N fixed reference points are selected, and the level set equation is decoupled into a series of products of CS-RBF and related expansion coefficients as follows.

where, ${g_i}\left( x \right)$ represents the ith set of CS-RBF related to spatial variables, and ${\alpha _i}\left( t \right)$ represent the expansion coefficient corresponding to the ith CS-RBF, which is closely related to virtual time $t$ in the iterative process. Since the CS-RBF matrix is only related to the type of RBF and the position of the reference point, ${g_i}\left( x \right)$ only needs to be calculated once in the whole iterative process when the reference point is fixed, and the level set equation is only affected by the expansion coefficient ${\alpha _i}\left( t \right)$. Therefore, the expansion coefficient ${\alpha _i}\left( t \right)$ is selected as the design variable of PLSM.

By substituting the level set function interpolated by the CS-RBF into the RDE, the following equations can be obtained

By solving Eq (12), the following equations can be obtained

where, $\Delta t$ is the time step of the equation evolution. According to Eq (13), the updated form of the expansion coefficient ${\alpha _i}\left( t \right)$ is shown as follows.

In Eq (3), the upper and lower limit constraints are imposed on the level set function to smooth the points near the structure boundary. Therefore, in the actual calculation, the level set function update formula after CS-RBF parameterization is shown as Eq (15).

2.3.   Volume constraint processing via fuzzy PID control algorithm

Fuzzy PID controller (Zadeh ^[23]) is composed of fuzzy controller and PID controller. The block diagram of the fuzzy PID controller is shown in Figure 3.

In Figure 3, $r\left( t \right)$ is the target setting value of the control system, and $y\left( t \right)$ is the measured value. Firstly, the input variables of $ER\left( t \right)$ and $EC\left( t \right)$ are transferred to the fuzzy controller for processing, and the output variables of $\Delta kp, \Delta ki, \Delta kd$ are obtained. Then, the output variables are transferred to the PID controller to modify the parameters $kp0, ki0, kd0$, and the modified PID parameters are used to control the system.

In structural topology optimization, the augmented Lagrange method is usually used to deal with volume constraints. However, at the end of the iteration process, the value of Lagrange multiplier cannot continue to increase. If the volume constraint still does not meet the convergence conditions, the volume fraction cannot be effectively reduced, and the iteration process will oscillate or even fail to converge.

Therefore, in order to obtain better convergence effect of volume constraint and solve the problem that Lagrange multiplier cannot be effectively increased, fuzzy PID control algorithm is introduced into structural topology optimization. Fuzzy PID controller can calculate the value of Lagrange multiplier according to the difference between the current volume fraction and the target volume fraction, thereby obtaining good convergence effect of volume constraint.

Fuzzy PID controller uses Mamdani fuzzy controller, and the triangular membership function is used to fuzzify the input variables because of its uniform distribution and high sensitivity in the universe. The Max-Min rule is employed for the fuzzy control decision-making, and the gravity center method is employed for the de-blurring.

The input variables of the fuzzy controller are: volume deviation $ER$ between the set volume fraction value and the current volume fraction, and volume deviation change rate $EC$. And the calculation formulas of input variables $ER$ and $EC$ are expressed respectively as follows.

where, $E{R_{iT}}$ is the volume deviation value in the $iT{\text{th}}$ iteration, and $E{C_{iT}}$ is the volume deviation change rate in the $iT{\text{th}}$ iteration.

The vector, consisting of the output variables $\Delta kp, \Delta ki, \Delta kd$ of the fuzzy controller, is defined as the correction value of PID controller, and the calculating formula of the output variables $\Delta kp, \Delta ki, \Delta kd$ of the fuzzy controller is expressed as follows.

where, $\Delta pid{\text{ = }}{\left[ {\begin{array}{*{20}{c}} {\Delta kp}&{\Delta ki}&{\Delta kd} \end{array}} \right]^T}$, $\alpha$ is the fuzzy rules defined according to Table 1, $evalfis$ is the MATLAB function, which is used to calculate the output variables.

In Table 1, $E{R_{iT}}$ and $E{C_{iT}}$ are the input variables of the fuzzy controller. $\Delta kp$, $\Delta ki$, and $\Delta kd$ are the output variables of the fuzzy controller. NS, NM, NB, Z, PS, PM, and PB indicate the seven values of triangular membership functions used to fuzzify the input and output variables following the fuzzy rules, which are the abbreviated forms of negative small, negative middle, negative big, zero, positive small, positive middle, positive big, respectively.

Then the output variables are used to modify the PID initial parameters $kp0, ki0, kd0$. The calculation formula of the modified PID initial parameters is formulated as follows.

The corrected PID parameters can be used to calculate the Lagrange multiplier and deal with the volume constraint. The calculation formula of the improved Lagrange multiplier is formulated as follows.

where, ${\lambda _{iT}}$ is the Lagrange multiplier at the $iT{\text{th}}$ iteration. By substituting the Lagrange multiplier updating formula proposed in Eq (19) into Eq (5), the volume constraint can be treated.

2.4.   Calculation of topological derivatives

The calculation of topological derivative plays an important role in the process of topology optimization. To simplify the calculating process of topological derivatives and improve the calculating efficiency of topology optimization, the topological derivative is calculated by using the variation of the objective function. The finite difference has the advantages of relatively simple equations and convenient calculation procedures, by which the grid can be easily divided and time can be saved. Moreover, it can make use of the structural topological meshes to easily construct high-precision schemes. Therefore, the finite difference is adopted for sensitivity calculation.

When only the fixed external load is considered, the objective function, the stiffness matrix and the global displacement vector of structure will change when the internal material is added or deleted. Therefore, according to the variation of the objective function before and after the insertion of hole, the calculating formula of objective function change $\Delta J$ can be obtained as follows.

where, $\Delta \boldsymbol{J}$ is the variation of the objective function; $\boldsymbol{K}$ is the stiffness matrix of structure; $\boldsymbol{U}$ is the global displacement vector; $\Delta \boldsymbol{K}$ is the variation of the stiffness matrix; $\Delta \boldsymbol{U}$ is the variation of the global displacement vector; $\Delta {\boldsymbol{U}^T}\left( {\boldsymbol{K} + \Delta \boldsymbol{K}} \right)\Delta \boldsymbol{U}$ is the high-order infinitesimal term, which can be ignored. Therefore, the calculating formula of the variation of objective function expressed by Eq (20), namely the calculating formula of the topological derivative, can be expressed as follows.

3.   Formulation of topology optimization problem

3.1.   Topology optimization for minimizing the structural compliance

The mathematical expression of the minimum compliance problem is formulated as follows.

where, $J$ represents the overall compliance value of the structure; $u$ is the actual displacement; $U$ is the allowable displacement in the displacement field; $v$ is the virtual displacement; $a\left( {u, v} \right) = l\left( v \right)$ is the weak form of the elastic balance equation; $a\left( {u, v} \right)$ is the internal force virtual work of the elastic body, and $l\left( v \right)$ is the virtual work of external force, which are expressed respectively as follows.

where, $E$ is the elastic modulus of the material, $\varepsilon$ is the strain tensor.

In this work, we mainly study the optimization problem in the static load equilibrium state where the virtual displacement $v$ and the actual displacement $u$ of the elastic equilibrium equation are equal. The elastic equilibrium equation always satisfies the constraint. Therefore, the mathematical model of Eq (23) can be transformed into a minimum compliance problem with only one volume constraint as follows.

Here, the Lagrange multiplier method is employed to transform Eq (24) into an unconstrained structural topology optimization problem, as shown in Eq (25).

where, $\lambda$ denotes the Lagrange multiplier.

3.2.   Implementation of the numerical algorithm

Figure 4 shows the flow chart of the PLSM based on RDE and fuzzy PID control algorithm.

The implementation process of the optimization algorithm shown in Figure 4 can be also described as the following iterative steps.

Step 1: Discretize the structure design domain into finite element meshes, define the initial state of level set function $\phi$, calculate the CS-RBF matrix $\boldsymbol{g}\left( x \right)$ at the reference point, and calculate the initial values of design variables $\boldsymbol{\alpha} \left( t \right)$.

Step 2: Initialize fuzzy PID controller: define the input variables $ER$ and $EC$, the output variables $\Delta kp, \Delta ki, \Delta kd$ of the fuzzy adaptive PID controller according to Eqs (16) and (17). Formulate the fuzzy rules, and set the PID initial parameters $kp0, ki0, kd0$.

Step 3: Prescribe the boundary constraints and load conditions of the structure, and calculate the element stiffness matrix.

Step 4: Start iteration, and set the number of iterations as $iT = 1:200$.

Step 5: Calculate the structural displacement matrix $\boldsymbol{U}$, the structural stiffness matrix $\boldsymbol{K}$, the value of structural compliance $J$, and the values of topological derivative according to Eq (21).

Step 6: Update the value of Lagrange multiplier ${\lambda _{iT}}$ according to Eq (19), and deal with the volume constraint.

Step 7: Update the design variables $\boldsymbol{\alpha} \left( t \right)$ and the level set functions $\boldsymbol{\phi}$ according to Eq (15).

Step 8: Set the convergence value as $\theta = 1e - 3$, and judge whether the optimization results satisfy the convergence condition $\int_D {{\text{d}\Omega } \leqslant {V_{max}}} \;$ and ${J_{iT}} - {J_{iT - 1}} \leqslant \theta$. If the convergence condition is not satisfied, repeat steps (5) to (8); if the convergence condition is satisfied, stop the iteration and output the optimization results, the corresponding optimization curves of objective function and volume fraction.

3.3.   Case study on two methods for processing volume constraint

In this section, to verify that using fuzzy PID control algorithm to calculate Lagrange multiplier can get good effect on processing volume constraint, topology optimization of the Michell beam structure by the PLSM based on the RDE is used for case study. Where, two methods, namely the fuzzy PID control algorithm and the augmented Lagrange method, are used for processing volume constraint, respectively.

The design domain and boundary conditions of the Michell beam structure are shown in Figure 5. The aspect ratio of the design domain is 5:2, which is divided into 100 times 40 quadrilateral meshes. The lower left boundary is fixed, whereas the lower right boundary is simply supported. The concentrated load of $F = 100$ is applied at the midpoint of the lower side. Elastic modulus of material is set as ${E_0} = 1$. To avoid the singularity of structural stiffness during calculation, elastic modulus of the void is set as ${E_{min}} = 1e - 9$. Poisson's ratio is set as $\mu = 0.3$. The volume fraction constraint is set as ${V_{max}} = 0.5$. The regularization parameter is set as $\tau = 2e - 4$. The initial values of fuzzy PID controller are set as $kp0 = 10, ki0 = 10, kd0 = 5$, respectively. Figure 5 shows the history curves of Lagrange multiplier and volume fraction obtained by the two different methods during topology optimization of the Michell beam structure. Where, Figure 6(a) shows the history curves of Lagrange multiplier and volume fraction obtained by the fuzzy PID control algorithm; and Figure 6(b) shows the history curves of Lagrange multiplier and volume fraction obtained by the augmented Lagrange method.

From Figure 6(a), it can be seen that: at the beginning of the iteration process, the Lagrange multiplier obtained by fuzzy PID control algorithm increases greatly; at the end of the iteration process, the growth rate of Lagrange multiplier gradually decreases; when the value of Lagrange multiplier remains stable, the volume fraction is optimized to its target value and remains stable.

The iterative format of the augmented Lagrange multiplier with the curves shown in Figure 6(b) is expressed by Eq (26).

According to Eq (26), the value of $\lambda$ is not zero at the beginning of iteration. During the iterative optimization process, the topology optimization of the structure is realized by continuously deleting the elements with small topological derivative.

From Figure 6(b), it can be seen that: at the beginning of the iteration process, the value of Lagrange multiplier gradually increases; at the end of the iteration process, the value of Lagrange multiplier oscillates obviously, which leads to the corresponding fluctuation in volume fraction. The fluctuation in volume fraction is not conducive to the convergence of topology optimization. Moreover, the volume fraction is finally optimized to 0.46, which is lower than its target value. According to the curves of Lagrange multiplier and volume fraction in Figure 6, it can be concluded that: compared with the augmented Lagrange method, fuzzy PID control algorithm has obviously good numerical stability and high optimization convergence efficiency. In other words, using fuzzy PID control algorithm to deal with volume constraints can get relatively stable Lagrange multiplier optimization curve, and make the optimization process of volume constraints relatively stable.

The final topology optimization results of the Michell beam structure obtained by the two different methods for processing volume constraint are shown in Figure 7. Where, Figure 7 (a) shows the final optimization results obtained by the fuzzy PID control algorithm, and Figure 7(b) shows the final optimization results by the augmented Lagrange method. From Figure 7, it can be concluded that the correct and reasonable structure can be obtained by both of the two methods for processing volume constraint.

3.4.   Selection of initial PID control parameters

In fuzzy PID control algorithm, the suitable initial PID control parameters can greatly reduce the calculating burden, accelerate the convergence speed and improve the convergence stability. In this paper, the initial PID control parameters are selected by the trial method. In other words, the initial PID control parameters are determined according to the specific optimization requirements and effect. To illustrate the influence of different values of $kp0, ki0, kd0$ on the optimization results, topology optimization of the cantilever beam structure by the proposed method is used as a case in this section.

The design domain and boundary condition of the cantilever beam structure are shown in Figure 8. The aspect ratio of the design domain is 5:3, which is divided into 100 times 60 quadrilateral meshes. The left side of the structure is fully fixed, whereas the other three sides are free. The concentrated load of $F = 100$ is applied at the midpoint of the right side. Elastic modulus of material is set as ${E_0} = 1$. To avoid the singularity of structural stiffness during calculation, elastic modulus of the void is set as ${E_{min}} = 1e - 9$. Poisson's ratio is set as $\mu = 0.3$. The regularization parameter is set as $\tau = 2e - 4$. The volume fraction constraint is set as ${V_{max}} = 0.5$.

To verify the influence of initial PID control parameters on the optimization results, the initial PID control parameters are taken as ($kp0{\text{ = }}5, ki0{\text{ = }}5, kd0{\text{ = }}5$), ($kp0{\text{ = }}5, ki0{\text{ = }}5, kd0{\text{ = 1}}5$), ($kp0{\text{ = }}5, ki0{\text{ = }}8, kd0{\text{ = 1}}5$) and ($kp0{\text{ = }}15, ki0{\text{ = }}8, kd0{\text{ = 1}}5$), respectively. Table 2 lists the numerical results during topology optimization of the cantilever beam structure with the four groups of different initial PID control parameters for comparison.

Figures 9, 10 and 11 show the history curves of Lagrange multiplier, objective function and volume fraction during topology optimization of the cantilever beam structure with the four groups of different initial PID control parameters, respectively. Where, the red curves a correspond to the initial PID control parameters $kp0{\text{ = }}5, ki0{\text{ = }}5, kd0{\text{ = }}5$, the blue curves b correspond to the initial PID control parameters $kp0{\text{ = }}5, ki0{\text{ = }}5, kd0{\text{ = 1}}5$, the green curves c correspond to the initial PID control parameters $kp0{\text{ = }}5, ki0{\text{ = }}8, kd0{\text{ = 1}}5$ and the black curves d correspond to the initial PID control parameters $kp0{\text{ = }}15, ki0{\text{ = }}8, kd0{\text{ = 1}}5$.

By comparing curve a and curve b in Figure 9, it can be found that the optimization history curve of Lagrange multiplier does not change obviously when only increasing the value of initial parameter $kd0$. By comparing curve b and curve c in Figure 9, it can be found that the growth rate of Lagrange multiplier is obviously accelerated when only increasing the value of initial parameter $ki0$. By comparing curve c and curve d in Figure 9, it can be found that the growth rate of Lagrange multiplier will also increase by a relatively small margin when only increasing the value of initial parameter $kp0$. Accordingly, it can be concluded that the value of Lagrange multiplier can be adjusted and controlled by changing the values of initial PID control parameters. Specifically, increasing the value of initial parameters $kp0$ and $ki0$ can accelerate the growth rate of Lagrange multiplier, and the change of initial parameter $ki0$ can control the growth rate of Lagrange multiplier more obviously.

From Figures 10 and 11, it can be found that: when the initial PID parameters are taken as $kp0{\text{ = }}5, ki0{\text{ = }}5, kd0{\text{ = }}5$, about 150 iterations are needed to make the optimization process converge, the final objective function value is 3.96 times 10⁵, and the optimization history curves of objective function and volume fraction have small fluctuations, resulting in long adjustment time and poor stability of optimization process.

By comparing curve a and curve b in Figures 10 and 11, it can be found that: when keeping the values of initial parameters $kp0$ and $ki0$ unchanged and only increasing the value of initial parameter $kd0$, it also takes about 150 iterations to make the optimization process converge, and the final objective function value is 3.99 times 10⁵. It is certified that increasing the value of initial parameter $kd0$ can make the optimization history curves of objective function and volume fraction more smooth and improve the stability of optimization process.

By comparing curve b and curve c in Figures 10 and 11, it can be found that when keeping the values of initial parameters $kp0$ and $kd0$ unchanged and only increasing the value of initial parameter $ki0$, the increasing rate of objective function and the decreasing rate of volume fraction are both obviously accelerated, it takes about 80 iterations to make the optimization process converge, the final objective function value is 3.99 times 10⁵, and the optimization history curves of objective function and volume fraction are relatively smooth, without wild fluctuations and oscillations. It is certified that increasing the value of initial parameter $ki0$ can reduce the number of iterations and improve the optimization efficiency.

By comparing curve c and curve d in Figures 10 and 11, it can be found that when keeping the values of initial parameters $ki0$ and $kd0$ unchanged and only increasing the value of initial parameter $kp0$, the increasing rate of objective function and the decreasing rate of volume fraction are both further accelerated, it takes about 50 iterations to make the optimization process converge, the final objective function value is 4.00 times 10⁵, and the optimization history curves of objective function and volume fraction are smooth, implying that the optimization process is stable. It is proved that increasing the value of initial parameter $kp0$ can also improve the optimization efficiency and obtain satisfactory results. The maximum difference between the final objective function values obtained by the four groups of initial PID control parameters is less than 2%, which proves that changing the values of the initial PID control parameters will not affect the objective function significantly.

Therefore, according to the specific effects of the initial PID control parameters on the optimization results, the values of the initial PID control parameters can be selected appropriately, so as to obtain stable and fast topology optimization process and accurate topology optimization results.

4.   Numerical examples

In this section, several typical 2D and 3D structures are optimized as numerical examples to illustrate the effectiveness of the proposed method. Taking the minimum compliance as the objective function and the volume fraction as the constraint condition, the optimization results obtained by the proposed method are analyzed and compared with those obtained by other methods to demonstrate the stability and effectiveness of this method.

4.1.   The L-beam structure

Figure 12 shows the design domain and boundary conditions of the L-beam structure, which is divided into 3200 quadrilateral meshes. The upper boundary is fixed, whereas a vertical downward force $F = 100$ is applied at the midpoint of the right boundary. The design domain is divided into 3200 elements. Elastic modulus of material is set as ${E_0} = 1$. To avoid the singularity of structural stiffness during calculation, elastic modulus of the void is set as ${E_{min}} = 1e - 9$. Poisson's ratio is set as $\mu = 0.3$. The volume fraction constraint is set as ${V_{max}} = 0.5$. The regularization parameter is set as $\tau = 2e - 4$. The initial values of fuzzy PID controller are set as $kp0 = 15, ki0 = 12, kd0 = 15$, respectively.

Figure 13 shows the topology optimization result of the L-beam structure obtained by the proposed method. Figure 14 shows the history curves of objective function and volume fraction during topology optimization of the L-beam structure by the proposed method.

From Figure 13, the reasonable topology configuration of the L-beam structure is obtained by the proposed method, and the structure boundary is smooth and clear. From Figure 14, the L-beam structure from initial volume fraction of 0.75 to target volume fraction of 0.4, it takes about 50 iterations, and the final value of objective function is 6.40 times 10⁵. In the optimization process, the optimization curves of objective function and volume fraction are relatively smooth, almost no fluctuations, and the optimization process is stable.

4.2.   The multi-load bridge structure

Figure 15 shows the design domain and boundary conditions of the multi-load bridge structure. The aspect ratio of the design domain is 8:3, which is divided into 80 times 30 quadrilateral meshes. The lower left boundary is fixed, whereas the lower right boundary is simply supported. A concentrated load of $F = 60$ is applied at the midpoint of the bottom of bridge, and two concentrated loads of $F = 30$ are applied at the one-fourth point and the three-fourth point of the bottom of bridge, respectively. Elastic modulus of material is set as ${E_0} = 1$. To avoid the singularity of structural stiffness during calculation, elastic modulus of the void is set as ${E_{min}} = 1e - 9$. Poisson's ratio is set as $\mu = 0.3$. The volume fraction constraint is set as ${V_{max}} = 0.5$. The regularization parameter is set as $\tau = 2e - 4$. The initial values of fuzzy PID controller are set as $kp0 = 15, ki0 = 12, kd0 = 15$, respectively.

To further verify the effectiveness of the proposed method, the proposed method, level set method based on RDE (Otomori et al. ^[28]) and PLSM based on CS-RBF (Cui et al. ^[33]) are respectively utilized to address topology optimization of the multi-load bridge structure, and the optimization results are analyzed and compared in the meanwhile. Figure 16 shows the topology optimization results of the multi-load bridge structure obtained by the above-mentioned three methods. Figure 17 shows the history curves of objective function during topology optimization of the multi-load bridge structure by the above-mentioned three methods. Figure 18 shows the history curves of volume fraction during topology optimization of the multi-load bridge structure by the above-mentioned three methods. In Figures 17 and 18, the red curve b represents the curve obtained by the proposed method, the blue curve c represents the curve obtained by the level set method based on RDE, and the green curve d represents the curve obtained by the PLSM based on CS-RBF. Table 3 lists the numerical results in topology optimization of the multi-load bridge by the three methods for comparison.

From the image optimization effect in Figure 16, the topological structures obtained by the three optimization methods are basically the same, and the correct and reasonable optimal structures can be obtained. However, compared with the level set method based on RDE, the structure boundary obtained by the proposed method is smoother and clearer, and compared with the PLSM based on CS-RBF has simpler topological configuration. According to the optimization curves in Figure 17, Figure 18 and the optimization results in Table 3, the proposed method needs about 50 iterations when the volume constraint converges, and the objective function value is 1.81 times 10⁵. Compared with PLSM based on CS-RBF, the objective function value is reduced by 4.7%. Compared with the level set method based on RDE, the objective function value is reduced by 11.1%, and better structural performance is obtained. It can be found that the introduction of virtual interface energy fully relaxes the topology optimization problem, obtains the optimal structure with smooth boundary, and improves the numerical stability of the optimization problem. Therefore, by comparing the three optimization methods, it can be concluded that the proposed method can achieve satisfactory optimization results in terms of optimization effect and computational efficiency.

4.3.   The 3D cantilever beam structure

Figure 19 shows the design domain and boundary conditions of the 3D cantilever beam structure. The ratio of the three dimensions of the design domain is 15:10:1, which is divided into 60 times 40 times 4 elements. The left boundary is fixed, whereas the right midpoint is applied with vertical downward load $F = 100$. Elastic modulus of material is set as ${E_0} = 1$. To avoid the singularity of structural stiffness during calculation, elastic modulus of the void is set as ${E_{min}} = 1e - 9$. Poisson's ratio is set as $\mu = 0.3$. The target volume fraction is set as ${V_{max}} = 0.5$. The regularization parameter is set as $\tau = 2e - 4$. The initial PID parameter are set as $kp0 = 40, ki0 = 20, kd0 = 15$, respectively.

The topology optimization result of the 3D cantilever structure obtained by the proposed method is shown in Figure 20. The history curves of the objective function and volume fraction during topology optimization of the 3D cantilever structure by the proposed method are shown in Figure 21.

It can be seen from Figure 20, using the algorithm proposed in this paper to optimize the 3D cantilever structure, clear structure boundary can be obtained, and there is no bad numerical phenomenon such as checkerboard. From the optimization curves shown in Figure 21, the objective function and the volume fraction change smoothly in the whole optimization process, and there is almost no fluctuation and oscillation phenomenon. And the final objective function value is 9.88 times 10⁵. Therefore, it can be concluded that although the number of 3D structure meshing is several times larger than 2D structure, the proposed algorithm can improve the optimization efficiency, ensure the stability of the optimization process, and get a clear structure configuration.

4.4.   The 3D arch bridge structure

Figure 22 shows the design domain and boundary conditions of the 3D arch bridge structure. The ratio of the three dimensions of the design domain is 15:5:1, which is divided into 60 times 20 times 4 elements. The displacement of nodes within two-unit length at the left and the right ends of the lower boundary of the structure is fully constrained, whereas the plate area with one-unit thickness on the upper surface of the structure is non-design area. The vertical downward load $F = 5$ is applied to all nodes on the upper surface of the structure. The design domain Elastic modulus of material is set as ${E_0} = 1$. To avoid the singularity of structural stiffness during calculation, elastic modulus of the void is set as ${E_{min}} = 1e - 9$. Poisson's ratio is set as $\mu = 0.3$. The target volume fraction is set as ${V_{max}} = 0.4$. The regularization parameter is set as $\tau = 2e - 4$. The initial PID parameter are set as $kp0 = 80, ki0 = 45, kd0 = 40$, respectively.

Figure 23 shows the topology optimization results of the 3D arch bridge structure from different angles, in which (a) is the front view and (b) is the side view. Figure 24 shows the history curves of objective function and volume fraction during topology optimization of the 3D arch bridge structure by the proposed method.

From Figures 23 and 24, it can be seen that the final result obtained by using the algorithm proposed in this paper to optimize the arch bridge structure is similar to the arch bridge structure in practical application, which proves that the method proposed in this paper has certain reference value in practical engineering application. Due to the existence of non- design area in the optimization process, it takes about 70 iterations to reduce the volume fraction from the initial value of 0.95 to the target value of 0.4, and the final objective function value is 1.10 times 10⁶. The results show that using the proposed algorithm for topology optimization of 3D structure can get the correct structure configuration, and the optimization efficiency is high, the convergence speed is high, and the iterative process is stable, which has a certain practicality in engineering applications.

5.   Conclusions

In this paper, a PLSM based on the RDE and the fuzzy PID control algorithm is proposed to solve structural topology optimization problem of minimum compliance under volume constraints. Firstly, the RDE is used as the evolution equation of the level set function, meanwhile, the topological derivative of the material domain is used as the reaction term of the RDE to drive the update of level set function, which can generate holes in the material domain and ensure the stability of numerical calculation. Secondly, in the process of solving the RDE, the CS-RBF is used to interpolate the level set function to correct the diffusion term in RDE, which can improve the efficiency and accuracy of calculation because the CS-RBF owns the advantages of sparse coefficient matrix and local correlation of parameters. Thirdly, the fuzzy PID control algorithm is used to adjust the value of Lagrange multiplier so as to control the volume constraint of structural topology optimization. The comparative examples illustrate that a more stable optimization history curve of Lagrange multiplier is obtained by the fuzzy PID control algorithm, which improves the stability of the convergence process of volume constraint as well as topology optimization. Thus, in addition to inheriting the merits of the PLSM, the proposed method takes the advantages of the RDE, the CS-RBF and the fuzzy PID control algorithm. By this method, not only can clear and smooth boundary of structural topology optimization result be obtained, but also the stability and efficiency of structural topology optimization process can be improved. Furthermore, this method has been applied to 3D structural topology optimization, which broadens its application scope and confirms its effectiveness in the meanwhile.

Acknowledgments

The authors would like to express thanks to the reviewers for their valuable suggestions for improving the manuscript. This work was supported by Natural Science Foundation of Shaanxi Province (2020JM-200) and China Scholarship Council (201506965015).

Conflict of interest

The authors declare there is no conflict of interest.