Topology Optimization for Minimum Compliance with Material Volume and Buckling Constraints under Design-Dependent Loads

Abstract

Stability is a critical factor in structural design. Although buckling-constrained topology optimization has been investigated in previous work, the problem has not been considered under design-dependent loads. In this study, a model of buckling constraints in topology optimization problems under design-dependent loads was proposed to solve the above problem. First, the Kreisselmeier–Steinhauser aggregation function was employed to reduce multiple constraints to a single constraint. Then, the problem was sequentially approximated using the optimality criteria method tailored to update the variables. After that, a gradient-based optimization algorithm was established based on finite element and sensitivity analyses for the topology optimization problem with design-dependent loads. Finally, four numerical examples with design-dependent loads were comparatively analyzed, with and without bucking-constrained solutions. The calculation results proved the effectiveness and reliability of the optimization algorithm. Therefore, in this study, it was suggested that the developed optimization algorithm gained improved applicability.

1. Introduction

Topology optimization is a method for solving design problems subject to multiple constraints (such as stiffness, mass, stress, stability, and weight) and providing optimal structural forms in a given design domain. In the past 30 years, topology optimization has been increasingly used as a structural design tool in engineering design problems, and increasing amounts of new methods are being developed [1,2,3]. However, previous topology optimization methods have focused on compliance minimization subjected to certain material volume constraints. A more realistic objective for a structure under pressure loads would be to minimize compliance while satisfying buckling constraints.

Incorporating buckling constraints into continuum topology optimization presents some challenges. Topology optimization problems typically have a large number of elements. After adding buckling constraints, large-scale eigenvalue equations need to be solved, which leads to a sharp increase in the amount of calculation. Furthermore, in eigenvalue-related topology optimization problems, there are usually numerical problems, such as artificial modes and local minima, that can seriously affect the convergence of the optimization problem.

There have been several studies on stability in the topology optimization design of continuum structures. Neves et al. [4] used the solid isotropic material with penalization (SIMP) interpolation model to study the problem of maximizing a buckling factor subjected to material volume constraints for the first time, and they pointed out that, due to the introduction of low-density elements to the SIMP model, there was an artificial mode problem when solving the eigenvalue problem. They proposed that the geometric stiffness matrix of the low-stress element should be set to zero. Neves [5] believed that this method would cause the algorithm to oscillate, and they proposed an improved SIMP interpolation model that interpolated the element stiffness matrix and the geometric stiffness matrix separately. Zhou [6] believed that the method required a reasonable selection of parameters to avoid artificial modes, so the same interpolation method was adopted for the element stiffness matrix and the geometric stiffness matrix. Sun et al. [7] used the ICM method to study the topology optimization problem with consideration of stability constraints. Based on the eigenvalue translation technique, Gao and Ma [8,9] proposed a new method for dealing with artificial modes. In combination with a two-stage algorithm, Gao and Ma studied buckling and material volume constrained compliance minimization problems. The coating topology optimization method was used for additively manufactured infill components by Clausen et al. [10]. The buckling performance of the optimized structure was significantly improved compared to that of the standard minimum compliance optimized structure, and it was pointed out that the improvement in the impact resistance of the porous structure was due to the increased number of structural components. In terms of stiffness and improving the specific bending stiffness, Ferrari and Sigmund [11] studied the problem of compliance minimization subjected to material volume constraints and buckling constraints and discussed the effects of different finite elements, sensitivities, and aggregation functions instead of a single buckling constraint on the optimization results. After that, Ferrari and Sigmund [12] proposed a method to approximate the calculation eigenvalue problem. Their method resulted in a remarkable reduction in the computational cost of eigenvalue calculation, and they applied the method applied to topology optimization problems subjected to buckling constraints of cantilever beams in 2D and 3D. Recently, Ferrari et al. [13] published an educational article considering buckling constraints, in which the MATLAB code was open source. Zhang et al. [14] adopted the feature-driven optimization method and the finite cell method to solve buckling-constrained topology optimization problems. Mendes et al. [15] adopted the topology optimization of binary structures method to solve structural optimization problems under buckling constraints and design-dependent loads.

This paper describes an investigation of the effects of buckling constraints on problems with design-dependent loads. For these problems, the structural loading changed in terms of magnitude, direction, and location as the design changed. In pressure loading, the most crucial task is to identify the material boundary location at which the pressure acts.

At present, research on the loading surface search algorithm is generally divided into two categories. The first category involves achieving boundary identification based on the element or node density. Hammer and Olhoff [16] proposed a boundary search method for pressure loading for the first time, and they obtained loading surfaces by identifying equal-density nodes. Du and Olhoff [17] improved their method by representing loading surfaces with straight-line segments instead of a B-spline, which made the loading surface identification process more robust. Zhang et al. [18] proposed a direct element loading method that uses a density threshold to classify boundary elements. In their method, the load was directly applied to the finite element boundary. The advantage of this method was that it was easy to implement, and load sensitivity analysis could be easily avoided. In addition, when using the above method for boundary search, it is necessary to determine the endpoint of the loading surface before optimization, which does not apply to the design of asymmetric structures under internal pressure. To solve this issue, Lee and Martins [19] introduced a void region to apply pressure, avoiding the pre-definition of endpoints. In addition, Lee’s work analytically deduced the sensitivity of the load. Recently, many researchers have solved mechanical problems through image processing [20,21]. In 2015, Zhao and Wang [22] adopted an image segmentation technique based on Distance Regularized Level Set Evolution [23] (DRLSE: Distance Regularized Level Set Evolution) for boundary recognition, and they successfully applied this method to the design of an underwater pressure hull. In 2018, Li et al. [24] proposed a simple algorithm based on digital image processing and area contour tracking technology that could generate suitable loading surfaces in the process of topology evolution. Since the endpoints of the loading boundary were determined semi-automatically, the scheme could handle both the loading from outside the domain and the pressure completely contained within the domain while avoiding the calculation of the load sensitivity.

The second category involves tracking changes in pressure loads by constructing multi-physics models. Chen and Kikuchi [25], as well as Bourdin and Chanbolle [26], applied a fluid-solid-void model to identify a loading surface. Sigmund and Clausen [27] suggested that a hybrid finite element could be used to simulate incompressible fluid structural regions. Some scholars have used the evolutionary structure method to study topology optimization problems with design-dependent loads. In 2014, Picelli et al. [28,29] extended the Bidirectional Evolutionary Structural Optimization (BESO) method using binary {0,1} design variables to solve design-dependent loads, acoustics, and fluid-structure coupling problems. In 2018, Sivapuram and Picelli [30] invented a binary structure topology optimization method using design variables 0 and 1 and integer mathematical programming, and they used this method to solve the fluid pressure load problem. The binary structure topology optimization method does not have the problem of blurred boundaries, which is its advantage over the SIMP method. Some scholars have also studied the loading surface search algorithm using the level set method. Zhao and Guo [31] and Jiang and Zhao [32] used level set evolution technology to propose a search method for designing the boundary of a load. When appropriating mathematical transformation, it is convenient to deal with the design-dependent loads imposed on the structure, avoiding tedious boundary extraction work in the density method. In 2015, Xia et al. [33] used two-level set functions to represent the free boundary and pressure. In 2018, Emmendoerfer et al. [34] adopted the level set method based on the reaction-diffusion equation to mark the pressure load boundary by tracking the level set function of the pressure load boundary. In 2019, based on the level-set method, Picelli et al. [35] adopted the Laplace equation to solve the pressure field of a fluid. Recently, Kumar [36] used Darcy’s law in rock mechanics to deal with design-dependent loads. Osezua et al. [37] proposed a new simple but practical boundary identification-load evolution (BILE) model to deal with structural optimization problems under pressure loads. Using this method, the design of an underwater pressure hull was explored. The considerable advantage of the topology optimization method involving pressure loads based on the BILE model is that it can complete the iteration within 80 to 100 iteration steps. The computational efficiency is greatly improved compared with the previous method, and the optimal configuration of the result is similar to that of the previous method. The density gradient method proposed by Wang and Qian [38] to track a loading surface was similar to a multi-physics model. Design-dependent loads were applied through the area integration of the density gradient without the need for a loading surface search.

The SIMP method is an optimization solver that is more mature than other methods, and it can solve a number of topology optimization problems. Although buckling-constrained topology optimization has been studied in previous work, the problem has not been considered under design-dependent loads using the SIMP method. Likewise, research has been conducted on topology optimization under design-dependent loads. Based on the above discussion, the boundary searching method plays a significant role in solving this problem, and the DRLSE loading surface search algorithm based on image segmentation, which is easier to apply, is used in this study. However, the objective has always been to obtain the maximum stiffness structure subjected to volume constraints. The goal of this research was to solve structural topology optimization problems considering buckling and volume fraction constraints under design-dependent loads. A comparison with the compliance minimization results is also provided in this paper. In addition, the method used in this work is preliminarily applied to designing underwater pressure structures.

This paper is organized as follows. Section 2 gives the problem statements for buckling-constrained compliance minimization. Section 3 describes the derivations of the expressions for the sensitivity of the constraint functions and objective functions. The methods for pressure loads are described in Section 4. The optimization algorithm is given in Section 5. In Section 6, numerical examples are presented, as well as an analysis of the differences between with and without bucking-constrained solutions. Section 7 provides concluding remarks.

2. The Topology Optimization Problem

2.1. Problem Formulation

The discrete mathematical model of the topology optimization formulation of a compliance minimization problem with consideration of the material volume and buckling constraints can be expressed as

$$\begin{array}{l}Find:\mathit{x}=\{x_1,x_2,\cdots ,x_n\}\\ min:C={\mathit{F}}^T\mathit{U}\\ s.t.:\{\begin{array}{c}\begin{array}{c}\mathit{K}\mathit{U}=\mathit{F}\\ \min \left|{\lambda }_j\right|\ge \underset{\_}{\lambda }>0\end{array}\\ \begin{array}{c}{\displaystyle {\sum }_{e=1}^{n}V(x_e)\le f{V}_0}\\ 0<x_{min}\le x_e\le 1\end{array}\end{array}\end{array}$$

(1)

where x_i (i = 1, 2, 3, …, n) are the design variables of the relative material density, n is the number of design variables, C is the compliance, F represents the nodal force vectors obtained after the boundary search, U is the corresponding global displacement vector, λ_j is the j-th buckling factor of the structure in the optimization process, λ is the lower bound of the first-order buckling factor of the structure, V is the total volume of the structure in the optimization process, V₀ is the volume of the initial design domain, f is the upper bound of the given volume fraction, and x_min is a non-zero lower bound introduced to ensure the positive definiteness of the total stiffness matrix of the structure. The model is the same as the models used in the literature [14,15].

For convenience of calculation, the Kresselmeier–Steinhauser (K-S) aggregation function is used to convert the multi-constraint problem of Equation (2) into a single-constraint problem13. The new topology optimization formulation can be expressed as:

$$\begin{array}{l}Find:\mathit{x}=\{x_1,x_2,\cdots ,x_n\}\\ min:C={\mathit{F}}^T\mathit{U}\\ s.t.:.\{\begin{array}{c}\begin{array}{c}\mathit{K}\mathit{U}=\mathit{F}\\ J_1^{KS}[g_V,g_{\lambda }](\mathit{x})\le 0\end{array}\\ 0<x_{min}\le x_e\le 1\end{array}\end{array}$$

(2)

where we used the K-S function to aggregate multiple constraints [13]. Explicitly, the constraint functions that are aggregated in Equations (1) and (2) can be expressed as:

$$\begin{array}{l}{g}_V(\mathit{x})=f(\mathit{x})/\overline{f}-1\\ g_{\lambda }(\mathit{x})=1-\underset{\_}{\lambda }{J}^{KS}[{\mu }_i](\mathit{x})\end{array}$$

(3)

where f(x) is the volume fraction of the optimized structure, f is the upper bound of the given volume fraction, and $J^{KS}\left[{\mu }_i\right]\left(x\right)$ is a K-S aggregation function defined below:

$$J^{KS}[{\mu }_i](\mathit{x})={\mu }_1(\mathit{x})+\frac{1}{\rho }\ln ({\displaystyle \sum _{i=2}^q{e}^{\rho ({\mu }_i(x)-{\mu }_1(x))}})$$

(4)

where $\rho$ is the aggregation parameter, $\rho \in \left[1,\infty \right]$, and ${\mu }_i$ is expressed as the $i$-$\mathrm{th}$ order eigenvalue. Equation (4) produces an upper bound and a smooth approximation for ${\mu }_1=\underset{i\in [1,12]}{\max }\left|{\mu }_i\right|$. Therefore, a smooth lower bound for λ is indirectly obtained.

2.2. Material Model

When using the SIMP model for topology optimization problems involving eigenvalue optimization, such as dynamic topology optimization problems and topology optimization problems with consideration of the stability constraints, artificial modes (also known as local modes) may appear in low-density areas, and these artificial modes will have an adverse impact on the optimization process. Many scholars have studied the problem of false modes and put forward corresponding solutions. To avoid artificial modes, in this research, the method proposed by Bendsøe and Sigmund 5 was adopted to calculate the element stiffness matrix and geometric stiffness matrix using different material interpolation models. When calculating the element stiffness matrix, the expression for the material interpolation model used is as follows

$$E(x)=E_{\min }+(E_0+E_{\min })x_{{}_e}^P$$

(5)

where E_min is a non-zero lower bound introduced to ensure the positive definition of the total stiffness matrix of the structure. In the optimization problem considering stability constraints, Bendsøe and Sigmund 5 suggested that this should be ${10}^{-6}$. E₀ represents the elastic modulus of solid material, x_e represents the density of the element, and P is the penalty factor, which is taken as three when interpolating the element stiffness.

When calculating the geometric stiffness matrix, the expression of the material interpolation model is as follows:

$$E_G=E_0{x}_{{}_e}^P$$

(6)

where P is the penalty factor, which is taken as three when interpolating the element stiffness.

In this research, three density fields were used to achieve the final optimal design. The density filtering method proposed by Bourdin [39] was used to avoid numerical instability, and the projection filter proposed by Wang et al. [40] was used to obtain a black and white design. In the expression below, x_e is the initial density field, ${\tilde{x}}_e$ is the density field filtered by the density filter, and ${\widehat{x}}_e$ is the density field filtered by the projection. The relationship between x_e and${\tilde{x}}_e$ can be expressed as follows:

$${\tilde{x}}_e=\frac{{\displaystyle {\sum }_{i=1}^m{h}_{e,i}{x}_i}}{{\displaystyle {\sum }_{i=1}^m{h}_{e,i}}}$$

(7)

where $h_{e,i}=\max \left(0,r_{min}-dist\left({\Omega }_i, {\Omega }_e\right)\right)$, $r_{min}$ is the filter radius, and $dist\left({\Omega }_i, {\Omega }_e\right)$ is expressed as the distance between the center of the element ${\Omega }_i$ and the element ${\Omega }_e$, and ${\widehat{x}}_e$ is obtained with a projection filter of ${\tilde{x}}_e$:

$${\widehat{x}}_e=\frac{\tanh (\beta \eta )+\tanh (\beta ({\tilde{x}}_e-\eta ))}{\tanh (\beta \eta )+\tanh (\beta (1-\eta ))}$$

(8)

where η and β are projection parameters, and $\eta \in \left[0,1\right]$, $\eta \in \left[1,\infty \right]$. Therefore, when calculating the objective function and the sensitivity of each constraint to the design variables, the chain rule needs to be adopted in the following form:

$$\frac{\partial C}{\partial x_e}=\frac{\partial C}{\partial {\widehat{x}}_e}\frac{\partial {\widehat{x}}_e}{\partial {\tilde{x}}_e}\frac{\partial {\tilde{x}}_e}{\partial x_e}$$

(9)

2.3. Linear Buckling Analysis

The traditional compliance minimization problem considering volume constraints only requires the equilibrium equation to be solved: $\mathit{K}\mathit{U}=\mathit{F}$. In addition to solving the equilibrium equation, the topology optimization problem, with consideration of the stability constraints, also requires the calculation of the buckling factor according to the eigenvalue equation. For linear buckling analysis, the eigenvalue equation of the structure can be expressed as:

$$(\mathit{K}+{\lambda }_i{\mathit{K}}_G){\mathsf{\phi }}_i=0$$

(10)

The eigenpairs $\left({\lambda }_i,\mathsf{\phi }\right)$ consist of the critical load factor ${\lambda }_i$ and the corresponding buckling modes ${\phi }_i$. To calculate the eigenvalue problem shown in Equation (10), the geometric stiffness matrix of the element needs to be solved, and its calculation formula is:

$${\mathit{K}}_G({\widehat{x}}_e)={\displaystyle \underset{\Omega }{\int }{\mathit{g}}^T\mathit{S}\mathit{g}}d\Omega$$

(11)

where S is the stress matrix, which can be defined as:

$$\mathit{S}={\mathit{I}}_{\mathbf{2}}\otimes {\mathsf{\sigma }}_e=\left[\begin{array}{cccc}{\sigma }_x& {\tau }_{xy}& 0& 0\\ {\tau }_{xy}& {\sigma }_y& 0& 0\\ 0& 0& {\sigma }_x& {\tau }_{xy}\\ 0& 0& {\tau }_{xy}& {\sigma }_y\end{array}\right]$$

(12)

where I₂ is a 2 × 2 unit matrix, and the calculation of the stress matrix depends on the solution of the element stress. The stress of the element can be expressed as:

$${\mathsf{\sigma }}_e=\left\{\begin{array}{c}\begin{array}{c}{\sigma }_x\\ {\sigma }_y\end{array}\\ {\tau }_{xy}\end{array}\right\}=\mathit{D}\mathit{B}\mathit{u}$$

(13)

where D represents the elastic matrix of the element (the two-dimensional plane stress elastic matrix), B represents the strain displacement matrix of the element, and $\mathit{U}$ represents the displacement vector of the element node. A quadrilateral iso-parametric element is used in this problem. The matrix g is defined as:

$$\mathit{g}=\left[\begin{array}{cc}\begin{array}{cccc}\frac{\partial N_1}{\partial x}& 0& \frac{\partial N_2}{\partial x}& 0\\ \frac{\partial N_1}{\partial y}& 0& \frac{\partial N_2}{\partial y}& 0\\ 0& \frac{\partial N_1}{\partial x}& 0& \frac{\partial N_2}{\partial x}\\ 0& \frac{\partial N_1}{\partial y}& 0& \frac{\partial N_2}{\partial y}\end{array}& \begin{array}{cccc}\frac{\partial N_3}{\partial x}& 0& \frac{\partial N_4}{\partial x}& 0\\ \frac{\partial N_3}{\partial y}& 0& \frac{\partial N_4}{\partial y}& 0\\ 0& \frac{\partial N_3}{\partial x}& 0& \frac{\partial N_4}{\partial x}\\ 0& \frac{\partial N_3}{\partial y}& 0& \frac{\partial N_4}{\partial y}\end{array}\end{array}\right]$$

(14)

where N_i is the shape function of the element. With matrix S and matrix g, the geometric stiffness matrix of the element can be obtained with Gaussian quadrature:

$${\mathit{K}}_G({\widehat{x}}_e)={\displaystyle \sum _{i=1}^4{\mathit{g}}^T\mathit{S}\mathit{g}}$$

(15)

3. Sensitivity Analysis

This section discusses the derivation of the sensitivities of the objective function, the K-S functions for buckling constraints, and the material volume constraint.

3.1. Sensitivity of Compliance

The sensitivity of the structural compliance to the design variables can be expressed as:

$$\frac{\partial C}{\partial {\widehat{x}}_e}=\frac{\partial {\mathit{F}}^T\mathit{U}}{\partial {\widehat{x}}_e}=\frac{\partial {\mathit{U}}^T\mathit{K}\mathit{U}}{\partial {\widehat{x}}_e}=\frac{\partial {\mathit{U}}^T}{\partial {\widehat{x}}_e}\mathit{K}\mathit{U}+{\mathit{U}}^T\frac{\partial \mathit{K}}{\partial {\widehat{x}}_e}\mathit{U}+{\mathit{U}}^T\mathit{K}\frac{\partial \mathit{U}}{\partial {\widehat{x}}_e}$$

(16)

According to KU = F, the following expression can be obtained:

$$\frac{\partial {\mathit{U}}^T}{\partial {\widehat{x}}_e}=(\frac{\partial {\mathit{F}}^T}{\partial {\widehat{x}}_e}-{\mathit{U}}^T\frac{\partial \mathit{K}}{\partial {\widehat{x}}_e}){\mathit{K}}^{-1}$$

(17)

$$\frac{\partial \mathit{U}}{\partial {\widehat{x}}_e}={\mathit{K}}^{-1}(\frac{\partial \mathit{F}}{\partial {\widehat{x}}_e}-{\mathit{U}}^T\frac{\partial \mathit{K}}{\partial {\widehat{x}}_e})$$

(18)

Substituting Equations (16)–(18), the following expression can be obtained:

$$\frac{\partial C}{\partial {\widehat{x}}_e}=2{\mathit{U}}^T\frac{\partial \mathit{F}}{\partial {\widehat{x}}_e}-{\mathit{U}}^T\frac{\partial \mathit{K}}{\partial {\widehat{x}}_e}\mathit{U}$$

(19)

For the DRLSE loading surface search algorithm based on image segmentation, it has been proven in the research of Wang et al. [40] that the sensitivity of the load term can be set to zero:

$$\frac{\partial \mathit{F}}{\partial {\widehat{x}}_e}=0$$

(20)

Therefore, the sensitivity of structural compliance to the design variables can be expressed as:

$$\frac{\partial C}{\partial {\widehat{x}}_e}=-{\mathit{U}}^T\frac{\partial \mathit{K}}{\partial {\widehat{x}}_e}\mathit{U}$$

(21)

3.2. Sensitivity of Buckling Constraint Functions

For convenience of calculation, this paper introduces the auxiliary variable μ = 1/λ. Then, the eigenvalue equation of the linear buckling analysis is:

$$({\mathit{K}}_G+{\mu }_i\mathit{K}){\mathsf{\phi }}_i=0$$

(22)

Referring to the paper by Rodrigues et al. [41], the sensitivity of a single critical buckling factor to design variables can be derived from the eigenvalue equation:

$$\frac{\partial {\mu }_i}{\partial {\widehat{x}}_e}=-{\mathsf{\phi }}_i^T(\frac{\partial {\mathit{K}}_G}{\partial {\widehat{x}}_e}-{\mu }_i\frac{\partial \mathit{K}}{\partial {\widehat{x}}_e}){\mathsf{\phi }}_i+{\mathit{v}}^T\frac{\partial \mathit{K}}{\partial {\widehat{x}}_e}{\mathit{u}}_0$$

(23)

where v^T is the adjoint displacement vector, which can be solved with the following adjoint equation:

$$\mathit{K}\mathit{v}={\mathsf{\phi }}_i^T\frac{\partial {\mathit{K}}_G}{\partial U}{\mathsf{\phi }}_i=\left\{\begin{array}{c}\begin{array}{c}{\mathsf{\phi }}_i^T\frac{\partial {\mathit{K}}_G}{\partial u_1}{\mathsf{\phi }}_i\\ ⋮\end{array}\\ {\mathsf{\phi }}_i^T\frac{\partial {\mathit{K}}_G}{\partial u_n}{\mathsf{\phi }}_i\end{array}\right\}$$

(24)

where n is the total number of degrees of freedom in the structural system, and the sensitivity of the geometric stiffness matrix to the design variables is:

$$\frac{\partial {\mathit{K}}_G}{\partial {\widehat{x}}_e}=\frac{\partial E_G}{\partial {\widehat{x}}_e}{\mathit{K}}_{\mathit{G}\mathbf{0}}$$

(25)

where K_G0 is the geometric stiffness matrix and ${\mathit{K}}_{\mathit{G}}$ is defined as $K_G=E_G\left(x\right)K_{G0}$, and the sensitivity of the geometric stiffness matrix to the element node displacement is:

$$\begin{array}{l}\frac{\partial {\mathit{K}}_G}{\partial u_i}(\mathit{u})=\frac{\partial }{\partial u_i}{\displaystyle \underset{{\Omega }_e}{\int }{\mathit{g}}^T\mathit{S}\mathit{g}}d{\Omega }_e={\displaystyle \underset{{\Omega }_e}{\int }{\mathit{g}}^T({\mathit{I}}_2\otimes \frac{\partial {\mathsf{\sigma }}_e}{\partial u_i})\mathit{g}}d{\Omega }_e\\ ={\displaystyle \underset{{\Omega }_e}{\int }{\mathit{g}}^T({\mathit{I}}_2\otimes \mathit{D}\mathit{B}\frac{\partial {\mathit{u}}_e}{\partial u_i})\mathit{g}}d{\Omega }_e={\displaystyle \underset{{\Omega }_e}{\int }{\mathit{g}}^T({\mathit{I}}_2\otimes \mathit{D}\mathit{B}{\delta }_{e,i})\mathit{g}}d{\Omega }_e\end{array}$$

(26)

If ${\mathit{u}}_{\mathit{i}}$ is the nodal displacement component of the i-th element, then it takes a value of 1; otherwise, it takes a value of 0. The rest of the vectors or matrices have already been described above and will not be repeated here.

Referring to the paper of Ferrari et al. [13], the first 12-order buckling factors of the structure are aggregated by the K-S aggregation function. The sensitivity of the buckling factor to the design variables after aggregation is:

$$\frac{\partial J^{KS}[{\mu }_i]}{\partial {\widehat{x}}_e}(\widehat{\mathit{x}})=\frac{{\displaystyle \sum _{i=1}^q{e}^{\rho ({\mu }_i(\widehat{x})-{\mu }_1(\widehat{x}))}}\frac{\partial {\mu }_i}{\partial {\widehat{x}}_e}(\widehat{\mathit{x}})}{{\displaystyle \sum _{i=1}^q{e}^{\rho ({\mu }_i(\widehat{x})-{\mu }_1(\widehat{x}))}}}$$

(27)

3.3. Sensitivity of Total Material Volume

The volume of the structure can be expressed as:

$$V={\displaystyle \sum _{e=1}^n{\widehat{x}}_e{V}_e}$$

(28)

where ${\widehat{x}}_e$ is the projection density field, V_e is the volume of the $e$-$\mathrm{th}$ element, and n is the total number of elements. The sensitivity of the structure volume to the design variables is:

$$\frac{\partial V}{\partial {\widehat{x}}_e}=V_e$$

(29)

4. Design-Dependent Loads

For pressure to act on a structure, a method of determining the boundary of the material surface is required, and the surface must be tracked as the design changes. In this section, we introduce the loading surface search algorithm based on DRLSE image segmentation technology, which was proposed by Li et al. [23].

The DRLSE method replaced the traditional level set formulation (LSF) with an essential distance regularization term to eliminate the irregularity of the level set function in the evolution process. Because the DRLSE method is very stable, it has been used in medical image processing, boundary recognition, and other fields. Zhao and Wang [22] introduced the DRLSE method for the topology optimization problem for design-dependent loads to obtain the boundary information of the structure during the optimization process. The steps of the algorithm are as follows:

(1)

Transfer the density matrix x to an 8-bit grayscale image.

(2)

Identify the boundary using the DRLSE method.

(3)

Extract boundary information and compute F.

4.1. Identify the Boundary Using the DRLSE Method

In the DRLSE model, the boundary of the gray image is represented as the zero-level contour of the final level set function. The grayscale image (0–1) obtained in topology optimization is actually a 1-bit grayscale image. To use the DRLSE algorithm, the image must first be turned into an 8-bit grayscale image:

I = 255 − 255*x

(30)

where $I$ is the 8-bit grayscale image data and x is the element node density matrix.

Since the boundary of the grayscale image is not obvious, the DRLSE algorithm does not necessarily find a suitable loading surface, which also makes it difficult to control the expansion of the loading surface during the optimization process. To overcome these types of problems, a density threshold is introduced, and the density of elements that is less than the density threshold is reset to zero during the image search. In this way, the obtained grayscale image has a clear boundary, and the expansion of the loading surface can also be controlled by adjusting the density threshold. The setting of the density threshold is similar to that in reference [19].

4.2. Design-Dependent Load Computing

In the process of finite element calculation, the load acting on the pressure-loading surface in the element must be transformed into the corresponding node. The pressure loading surface obtained with DRLSE is shown in Figure 1a. To simplify the calculation, in this research, the method recommended by Du and Olhoff [17] was adopted, and the line segment between the two intersections of the loaded surface and the element boundary was used as the new loading surface (Figure 1b). The process of converting the ballast into the nodal force in the pressure action element is shown in Figure 2, Figure 3 and Figure 4.

Assuming that the coordinates of the two intersection points of the pressure loading surface and the element boundary are $\left(x_1,y_1\right)$ and $\left(x_2,y_2\right)$, the equivalent nodal load matrix $f_e$ acting on the element can be determined with the element shape function N and the load matrix $ℱ$. The integral over the action boundary is obtained [19]. Then, the following expression is obtained:

$$f_e=\underset{L}{{\displaystyle \int }}{N}^Tℱ\mathrm{d}l=P{{\displaystyle \int }}_{x_1}^{x_2}{N}^T\left\{\begin{array}{c}-\frac{\mathsf{\Delta }y}{\mathsf{\Delta }x}\\ 1\end{array}\right\}dx$$

(31)

Converting the integral interval from $\left[x_1,x_2\right]$ to [−1, 1], the above formula can be calculated using the Gaussian integral formula, and the following expression is obtained:

$${\mathit{f}}_e=\frac{L}{2}{\displaystyle {\int }_{-1}^1{\mathit{N}}^{T}d\xi \left\{\begin{array}{c}-\Delta y/L\\ \Delta x/L\end{array}\right\}}=\frac{1}{2}{{\displaystyle \sum _i{\omega }_i\mathit{N}({\mathit{g}}_i^{’})}}^T\left\{\begin{array}{c}-\Delta y/L\\ \Delta x/L\end{array}\right\}$$

(32)

where P represents the magnitude of the pressure load, and ω_i and g_i are the i-th weight function and the corrected Gaussian quadrature point, respectively. Letting $\Delta x=x_1-x_2,\overline{x}=\frac{x_2-x_1}{2}$, and similarly for $y$, then:

$${\mathit{g}}_i^{’}=\infty \left\{\begin{array}{c}\frac{1}{2}{g}_i\Delta x+\overline{x}\\ \frac{1}{2}{g}_i\Delta y+\overline{y}\end{array}\right\}$$

(33)

After the above calculation, the pressure on the loading surface in the element is equivalent to the nodal load, and finite element analysis can then be carried out.

5. Optimization Algorithm

Based on the above finite element and sensitivity analyses, the topology optimization problem in Equation (1) can now be solved using a gradient-based optimization algorithm. In this research, the optimality criteria update optimization solver proposed by Ferrari et al. [13] and the parameters recommended by Ferrari et al. [13] were adopted. The flowchart of the proposed iterative solution algorithm is shown in Figure 5.

Step 1 Parameter setting: Define the number of finite elements in the design domain and the material properties, such as the Poisson’s ratio, elastic modulus, β, η, filter radius, volume constraints, and buckling constraints.

Step 2 Finite element preparation: Prepare the element stiffness matrix, node number of the geometric stiffness matrix, preparation matrix for density filtering, etc.

Step 3 Boundary condition setting: Set the boundary, load, and other conditions.

Step 4 Initialize the density field: Calculate the physical density by using linear density filters and projection filters.

Step 5 Boundary search: Based on the DRLSE loading surface search algorithm, search for the position where the force is applied and obtain the force vector.

Step 6 Finite element analysis: Perform finite element analysis.

Step 7 Sensitivity analysis: Calculate the sensitivities of the objective function, the K-S functions for the buckling constraints, and the material volume constraint.

Step 8 Update of the density field: Update the density field.

Step 9 Judgment of convergence condition: After judging whether the convergence condition is reached, the convergence automatically changes the buckling constraints to continue the optimization.

Step 10 Repeat or End: Repeat steps 4–10 and continue to iterate if not converged. Otherwise, end the loop and post-process the optimization results.

6. Numerical Examples

In this section, the effectiveness of the topology optimization model proposed in this paper is illustrated with numerical examples, such as arch structures subjected to external pressure. First, the initial solution was obtained with the compliance minimization model, which only considered volume constraints. Then, the buckling factor of the optimization result was used as a reference to determine the initial value of the buckling constraints when the topology optimization problem, with consideration of the buckling constraints, was performed. In the following numerical examples, Young’s modulus was set to 1, and Poisson’s ratio was set to 0.3. The optimization model was realized with MATLAB programming. The MATLAB version was R2020b. The eigs function in MATLAB was used to solve the structural eigenvalues and eigenvectors. The Lanczos algorithm was used to solve the symmetric matrix. The computer CPU used was an Intel(R) Xeon(R) E5-2678, the main frequency was 2.6 GHz, and the computer memory was 64 GB.

6.1. Arch Structure Subjected to External Pressure

First, a topology optimization design for the compliance minimization problem subjected to volume constraints was carried out. The final optimization result is shown in Figure 4. The area was divided into a 100 × 50 grid. The objective function value of the optimized structure was 2.44 × 10⁻², and the first-order buckling factor was 0.74. The corresponding convergence curve is shown in Figure 5.

Then, the topology optimization design subjected to buckling constraints was carried out, and the buckling constraints were taken as 0.9, 1.0, 1.2, and 1.3. The results are shown in Figure 6. The corresponding first-order buckling factors and compliance values were 0.91 and 2.30 × 10⁻², 1.01 and 2.35 × 10⁻², 1.21 and 2.58 × 10⁻², and 1.30 and 2.83 × 10⁻². In this numerical example, each iteration step lasted 0.8 s. To obtain the relationship between the first-order buckling factor and compliance,

Table 1. Table for buckling constraints and structural compliance.

The Value of the Buckling Restraint	The First Order Buckling Factor of the Optimization Result	The Compliance of the Optimization Result
0.90	0.91	2.30 × 10−2
1.00	1.01	2.35 × 10−2
1.20	1.21	2.58 × 10−2
1.30	1.30	2.83 × 10−2

was established.

Table 1 shows that when the volume fraction of the structural was constant, given the different buckling restraint values, the compliance of the structure increased with the increase of the buckling restraint value, which indicated that the improvement of the structural stability was at the expense of a small amount of structural stiffness. The above conclusion was the same as that in the literature 8. Figure 6 presents the intermediate topologies when the buckling constraint was 1.0. Figure 7 shows that when the buckling constraint was considered, the structure no longer maintained a complete arch structure, and the left and right supports expanded outward. When the buckling constraint was small, holes were generated inside. As the buckling constraint increased, the left and right supports on both sides changed into a porous structure, and there were more support rods to improve the stability of the structure, which was similar to Clausen et al. [10] using a porous filling structure to improve stability.

Figure 8 shows the optimization results for the compliance minimization subjected to the volume constraint and the first mode of the optimization results when the buckling constraint was 1.0. Since the low-order mode accounted for a large proportion of the response, and the low-order mode would be of more concern in engineering practice, only the first-order mode is given in this paper. Comparing the first-order modes of the two optimization results, it could be found that the stability was improved by expanding the left and right supports outward and adding holes when considering buckling constraints. Although the deformation region of the results considering buckling constraints was similar to the compliance minimization result, the buckling factor values were greatly increased, and the structure was more stable.

The final optimization results in this paper are consistent with the results of Mendes et al. [15], which are compared in Figure 9. Figure 10 and Figure 11 show the convergence history of the objective function, the structural volume fraction, and the first-order buckling factor when the buckling constraint was taken as 1.0. It can be seen from the convergence history that the proposed topology optimization model and the adopted optimization algorithm could stably converge when the topology optimization problem was solved with consideration of the stability constraints, which verified the feasibility of the model and algorithm.

6.2. Piston Design Subjected to Hydrostatic Pressure

A dam, also known as a piston structure, is a common structure under design-dependent loads. The initial design domain, initial pressure, and boundary conditions in this research are shown in Figure 12. The design area was a rectangular area of 6 m × 2 m. The structure was fixed in place at the midpoint of the bottom and horizontally fixed in place on the left and right sides. The value of the hydrostatic pressure was 6 × 10⁻³. The area was divided into a 120 × 40 grid. The design material volume fraction was 0.3 of the design area.

First, the topology optimization design of the compliance minimization problem subjected to volume constraints was examined. The final optimization result is shown in Figure 13. The objective function value of the optimized structure was 2.36 × 10⁻², and the first-order buckling factor was 0.84. The corresponding convergence curve is shown in Figure 14. In this numerical example, each iteration step lasted 1.8 s.

Then, the topology optimization design subjected to buckling constraints was examined. The results are shown in Figure 15. The corresponding first-order buckling factors and compliance values were 1.23 and 2.41 × 10⁻², 1.42 and 2.46 × 10⁻², 1.62 and 2.47 × 10⁻², and 1.82 and 2.53 × 10⁻², as shown in

Table 2. Table for the buckling constraint and structural compliance.

The Value of Buckling Restraint	The First Order Buckling Factor of Optimization Result	The Compliance of Optimization Result
1.20	1.23	2.41 × 10−2
1.40	1.42	2.46 × 10−2
1.60	1.62	2.47 × 10−2
1.80	1.82	2.53 × 10−2

. From Table 2, one can obtain the same conclusion as that discussed in Section 6.1.

Figure 15 shows that the optimization results with consideration of the buckling constraint had a significant difference from the compliance minimization optimization results with consideration of the volume constraints. When the buckling constraint was 1.2, holes were generated at the left and right ends of the bottom, and the number of support rods increased. With the increase of the buckling constraint, the material at the rod connecting the bottom and the left and right ends expanded outward, the amount of base material decreased, the middle holes became wider, and the rods were slender. When the constraint value was increased to 1.6, new holes and small rod supports were generated. If it was further increased, more holes and small rods were generated on the left and right tie rods. This was the same as the change rule of the arched structure subjected to external pressure.

The final optimization results in this paper are consistent with the results of Mendes et al. [15], which are compared in Figure 16. Figure 17 shows the optimization results for compliance minimization subjected to the volume constraint and the first mode of optimization results when the buckling constraint was 1.6. Comparing the first-order modes of the two optimization results, we obtained the same conclusion as that discussed in Section 5.

6.3. Structure Design Subjected to Hydrostatic Pressure

A ring-stiffened cylindrical hull is a commonly used pressure hull. Since the shape along the length of the cylinder is basically unchanged, in this research, a simplified plane strain model was adopted to study the optimal design of an underwater pressure hull. The initial design domain, initial pressure constraints, and boundary conditions are shown in Figure 18. To improve calculation efficiency, a 1/4 model was used for the calculation. The design area was a square area with dimensions of 1 m × 1 m. The value of the hydrostatic pressure was 6.5 × 10⁻³. The area was first divided into an 80 × 80 grid. The design material volume fraction was 0.1 of the design area.

First, a topology optimization design for the compliance minimization problem subjected to volume constraints was carried out. The final optimization result is shown in Figure 19. The first-order buckling factor of the optimized structure was 0.70. The optimization results subjected to different buckling constraints are shown in Figure 18.

Figure 20 shows that when the buckling constraints were considered, the result was significantly different from the structure for which compliance was minimized with only the volume constraints considered. When the buckling constraint was 0.80, the material moved toward the left and bottom constraints and toward the external expansion. When the buckling constraint increased to 1.7, a porous structure was formed on the left and right supports, and support rods were created to improve the stability of the structure. There were similarities in the arched structure subjected to external pressure.

The above final optimization structures were solved with 80 × 80 mesh grids, and the optimization results with large-scale mesh grids (300 × 300) subjected to different buckling constraints are shown in Figure 21. The results with the same buckling factor limit λ = 1.5 under different mesh resolutions were compared, and the four mesh grid divisions are 80 × 80, 200 × 200, 300 × 300, and 400 × 400, respectively. The optimization results are shown in Figure 22 and could prove the validity of the method under large-scale mesh grids.

6.4. Multiple Intersecting Spherical Pressure Hulls

In practical applications, a pressure hull requires a certain intermediate blank area to be reserved for storing equipment and liquids [42]. Therefore, in this example, the optimal design of the axisymmetric model of the pressure hull was carried out. The initial design domain, initial pressure constraints, and boundary conditions are shown in Figure 23. To improve the calculation efficiency, the 1/4 model was used for the calculation. The design area was a rectangular area with dimensions of 2.5 m × 1 m. The value of the hydrostatic pressure was 2 × 10⁻². The area was divided into a 90 × 36 grid. The design material volume fraction was 0.3 of the design area.

First, the topology optimization design for the compliance minimization problem subjected to volume constraints was carried out. The final optimization result is shown in Figure 24. The first-order buckling factor of the optimized structure was 0.80. The optimization results subjected to different buckling constraints are shown in Figure 21.

The optimization results with different constraints are shown in Figure 25. When the buckling constraint was 0.85, the material moved to the constraint at the bottom. As the buckling constraint increased, the material gathered at the left constraint and generated holes or support rods. When the buckling constraint was increased to 1.00, two support rods were created, and a hole and a support rod were created at the bottom-right constraint, which was similar to the three previous examples.

The above final optimization structures were solved with 90 × 36 mesh grids, and the optimization results with large-scale mesh grids (270 × 108) subjected to different buckling constraints are shown in Figure 26. Because of the high concentrations of stress in the corners of the structures, with the increase in the scale of the mesh grid division, the stress in the corners increases. This results in an increase in the value of the determinant of the geometric stiffness matrix, so the values of the buckling factor become smaller. For a certain pressure structure, the critical buckling loads were constant and equal to the buckling factor, timing the hydrostatic pressure. Therefore, to highlight the clear change in the structures with increasing buckling factor values, hydrostatic pressure was set to 2 × 10⁻⁴ in this numerical example, and the buckling factor was set to 25–29.5.

7. Conclusions

In this paper, an algorithm for topology optimization for minimum compliance with material volume and buckling constraints under design-dependent loads was developed. The structural topology optimization problems were solved using the SIMP method along with a pressure boundary searching scheme based on the DRLSE model to deal with design-dependent loads. The K-S aggregation function was used to reduce multiple constraints to a single constraint and replace the original buckling constraints, which not only improved the smoothness of the original buckling constraints and caused the solution to converge to a certain extent but also reduced the number of constraints.

The numerical examples demonstrated the effectiveness and reliability of the algorithm. The results of the arch structure and piston examples showed that the compliance minimization problems subject to buckling and volume constraints and the compliance minimization problems subject to volume constraints could result in considerably different topologies. The compliance minimization subject to buckling and volume constraints generally resulted in a more fully stable design when compared to the compliance minimization with the same volume constraints for which buckling constraints were not considered. In engineering structure design, stability is a critical factor that cannot be ignored. Hence, the buckling constrained approach is of importance in topology optimization.

The application of the algorithm in practical engineering design is the design of underwater structures, and the optimal results with different buckling factor limits are shown in this work. The results showed that more sandwich support structures were produced with increasing buckling factor limits. In addition, the method could be used for designing multiple intersecting spherical pressure hulls that have complex boundary conditions. Therefore, the developed algorithm has better applicability.

The future direction of this research is to simultaneously consider more design constraints, such as stress constraints and buoyancy constraints. In addition, the topology optimization problem will be extended into three dimensions in future work.