## 2 The Topology Optimization Problem
### Linear Buckling Analysis  
The traditional compliance minimization problem considering volume constraints only requires the equilibrium equation to be solved: \(\textbf{K}\textbf{U}=\textbf{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:  
\[(\textbf{K}+\lambda_{i}\textbf{K}_{G})\boldsymbol{\varphi}_{i}=0 \tag{10}\]The eigenpairs \((\lambda_{i},\mathbf{\varphi})\) consist of the critical load factor \(\lambda_{i}\) and the corresponding buckling modes \(\varphi_{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:  
\[\mathbf{K}_{G}(\hat{x}_{e})=\int\limits_{\Omega}\mathbf{s}^{T}\mathbf{S}\mathbf{g}d\Omega \tag{11}\]  
where \(\mathbf{S}\) is the stress matrix, which can be defined as:  
\[\mathbf{S}=\mathbf{I_{2}}\otimes\mathbf{\sigma}_{e}=\begin{bmatrix}\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{bmatrix} \tag{12}\]  
where \(\mathbf{I_{2}}\) is a \(2\times 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:  
\[\mathbf{\sigma}_{e}=\begin{Bmatrix}\sigma_{x}\\ \sigma_{y}\\ \tau_{xy}\end{Bmatrix}=\mathbf{D}\mathbf{B}\mathbf{u} \tag{13}\]  
where \(\mathbf{D}\) represents the elastic matrix of the element (the two-dimensional plane stress elastic matrix), \(\mathbf{B}\) represents the strain displacement matrix of the element, and \(\mathbf{U}\) represents the displacement vector of the element node. A quadrilateral iso-parametric element is used in this problem. The matrix \(\mathbf{g}\) is defined as:  
\[\mathbf{g}=\begin{bmatrix}\frac{\partial N_{\text{i}}}{\partial\mathbf{x}}&0&\frac{ \partial N_{\text{i}}}{\partial\mathbf{x}}&0&\frac{\partial N_{\text{i}}}{ \partial\mathbf{x}}&0&\frac{\partial N_{\text{i}}}{\partial\mathbf{x}}&0\\ \frac{\partial N_{\text{i}}}{\partial\mathbf{y}}&0&\frac{\partial N_{\text{i}}}{ \partial\mathbf{y}}&0&\frac{\partial N_{\text{i}}}{\partial\mathbf{y}}&0&\frac{ \partial N_{\text{i}}}{\partial\mathbf{y}}&0\\ 0&\frac{\partial N_{\text{i}}}{\partial\mathbf{x}}&0&\frac{\partial N_{\text{i}}}{ \partial\mathbf{x}}&0&\frac{\partial N_{\text{i}}}{\partial\mathbf{x}}&0&\frac{ \partial N_{\text{i}}}{\partial\mathbf{x}}\\ 0&\frac{\partial N_{\text{i}}}{\partial\mathbf{y}}&0&\frac{\partial N_{\text{i}}}{ \partial\mathbf{y}}&0&\frac{\partial N_{\text{i}}}{\partial\mathbf{y}}&0&\frac{ \partial N_{\text{i}}}{\partial\mathbf{y}}\end{bmatrix} \tag{14}\]  
where \(\mathbf{N_{\text{i}}}\) is the shape function of the element. With matrix \(\mathbf{S}\) and matrix \(\mathbf{g}\), the geometric stiffness matrix of the element can be obtained with Gaussian quadrature:  
\[\mathbf{K}_{G}(\hat{x}_{e})=\sum_{i=1}^{4}\mathbf{g}^{T}\mathbf{S}\mathbf{g} \tag{15}\]