Topology Optimization of Geometrically Nonlinear Structures Based on a Self-Adaptive Material Interpolation Scheme

Abstract

In this paper, a simple and effective self-adaptive material interpolation scheme is proposed to solve the numerical instability problem, which may occur in topology optimization considering geometrical nonlinearity when using density-based method. The primary concept of the proposed method revolves around enhancing the deformation resistance of minimum-density or intermediatedensity elements, thus avoiding numerical instability due to excessive distortion of these elements. The proposed self-adaptive material interpolation scheme is based on the power law method, and the stiffness of minimum-density or intermediate-density elements can be adjusted by a single parameter, α. During the optimization process, the parameter α will be changed according to an adaptive adjustment strategy to ensure that elements within the design domain are not excessively distorted, while the mechanical behavior of the structure can be approximated with acceptable accuracy. Numerical examples of minimizing compliance and maximizing displacement of structure are given to prove the validity of the proposed self-adaptive material interpolation scheme.

1. Introduction

Topology optimization is a well-known structural design method that aims to seek an optimized material layout within a certain design domain under given constraints and boundary conditions. Topology optimization is based on the finite element method and numerical optimization method, which enable it to provide more novel and efficient designs without relying on the designer’s experience. In 1988, Bendsøe and Kikuchi [1] introduced the homogenization method and used it to realize the structural topology optimization. Since then, the topology optimization method has attracted a lot of attention from researchers. After years of development, researchers have proposed various kinds of topology optimization approaches, such as solid isotropic material with penalization (SIMP) approach [2], the level-set method (LSM) [3], moving morphable components/void (MMC/V) method [4,5], evolutionary structural optimization (ESO) and bi-directional evolutionary structural optimization (BESO) method [6,7]. Sigmund and Maute [8], Xia et al. [9] and Zhu et al. [10] have provided detailed reviews of existing topology optimization methods.

Many research works on topology optimization do not take nonlinearity into consideration which, however, are only applicable to small deformation cases. In most practical engineering cases, such as compliant mechanism design, the structure tends to have large deformations. Under these circumstances, the finite element analysis with linear assumption can not provide an accurate prediction of the structural response, which would provide incorrect sensitivity information and thus lead to unreasonable optimization results. Therefore, considering the nonlinear effects during the topology optimization process is necessary.

Tracing back to the last decade of 20th century, Bendsoe [11], Jog [12], and Bruns and Tortorelli [13] have introduced the nonlinear assumption into the topology optimization. Buhl et al. [14] formulated topology optimization model based on the SIMP approach while considering geometrically nonlinearity. Their numerical experiment result revealed that in the cases of small displacements, the differences in optimization results due to nonlinearity are small. However, in some situations, such as when buckling effects occur, the effect of nonlinearity may be significant.

In recent years, researchers have devoted more and more efforts to studying nonlinear topology optimization problems. Bruns et al. [15] and Gea et al. [16] proposed topology optimization frame works of elastic structures and compliant mechanisms design with geometrically nonlinearities based on the SIMP approach. Based on the level set method, Kwak and Cho [17] introduced a topology optimization method for geometrically nonlinear structures. Huang and Xie [18] applied the BESO method to structure stiffness optimization problem with a combination of geometrically and materially nonlinear modeling. Xue et al. [19] considered nonlinear structural response analysis and developed an MMV-based method with an explicit B-spline description.

At present, the power-law penalization in SIMP is a material interpolation scheme that is widely used in addressing the topology optimization problem [5,20]. With this method, the material properties are formulated as a function of relative density, which is also known as the design variable. In general, the relative material density take values ranging from 0 to 1. Elements with a material density of 1 represent solid elements and elements with a material density of 0 represent void elements. The elements with intermediate-density materials would be penalized so that the most elements in the design domain could evolve to solid or void elements [2]. But, during or after the topology optimization procedure, it is found that elements with intermediate-density cannot be completely removed from the design domain. When performing nonlinear topology optimization, the “numerical instability” problem tends to occur in the region where the elements have low density, especially in cases of large deformation. This phenomenon results from excessive element distortion. The excessive distortion would make the elemental tangent stiffness matrix turn to indefinite or even negative definite, which result in failure of finite element analysis [21,22,23].

Researchers have put a lot of effort into solving the numerical instability problem. As the low density elements are easy to distort excessively, Bruns et al. [24] proposed method to remove them from the mesh that the finite element analysis is applied. Meanwhile, the low density elements would be added back according to the proposed strategy. Yoon et al. [25] presented a element connectivity parameterization method, in which the connectivity between elements is regarded as design variables instead of element density. And all the elements within design domain are defined as solid elements. Wang et al. [23] tried to address the problem from the aspect of energy and proposed a new energy interpolation scheme. In this method, the numerical instability is alleviated by modeling the low density elements with linear theory. Luo et al. [22] proposed a additive hyperelasticity technique. Its main idea is to add an additional layer of hyperelastic material on the top of original design domain to increase its stiffness, thus preventing low-density elements from distorting excessively. Chen et al. [21] present an educational nonlinear topology optimization Matlab codes that combines additive hyperelasticity techniques with the well-known finite element analysis software, ANSYS 17.1, to make it more suitable for beginners.

In this research, the numerical instability problem in the topology optimization considering geometrical nonlinearity is studied. The St. Venant Kirchhoff model can be used to model nonlinear materials, and has been widely applied in the nonlinear topology optimization. However, the St. Venant Kirchhoff model is based on a linear material model and is not well suited for situations where large deformations of the material occur. Therefore, the hyperelastic material model is applied in the following research instead of St. Venant Kirchhoff model [21]. Then, a simple but effective self-adaptive material interpolation scheme is constructed to address the problem of numerical instability. The main idea is to dynamically adjust the properties of minimum-density or intermediate-density elements within the design domain to prevent excessive distortion. Based on the proposed scheme, the adjoint method is employed to calculate the sensitivity numbers. Finally, the Method of Moving Asymptotes (MMA) is used to find the solution.

The organization of the following content is shown as follows: In Section 2, the general topology optimization model is illustrated. In Section 3, an effective material interpolation scheme for nonlinear topology optimization is proposed. In Section 4 and Section 5, the formulations of the finite element analysis and sensitivity analysis are constructed. Numerical examples concerning structural stiffness design and compliant mechanism design are given for illustrating the effectiveness of the material interpolation scheme in Section 6. Finally, the conclusions are drawn in Section 7.

2. Topology Optimization Formulation

In the research of topology optimization, the displacement maximization problem and compliance minimization problem are two typical topology optimization problems. Therefore, these two problems are considered in this research. One of the important parts of topology optimization is the finite element theory, which provides the basis for optimization. Before the optimization, the design domain is meshed into elements. The density of elements is regarded as variable and would change during the optimization. The compliance minimization problem and the displacement maximization problem can be formulated as the same form, which is presented as follows:

$$\begin{array}{cc}\hfill min& :\mathit{C}={\mathit{l}}^T\mathit{U}\hfill \\ \hfill & \hfill \\ \hfill \mathrm{s}.\mathrm{t}.& :\mathit{R}(\mathit{x},\mathit{U})=0\hfill \\ \hfill & :V\le V^{*}\hfill \\ \hfill & :0<x_{min}\le x_i\le 1\hfill \end{array}$$

(1)

where $\mathit{C}$ denotes the objective value. $\mathit{x}$ denotes the design variable vector, which represents the density of elements. $\mathit{R}(\mathit{x},\mathit{U})$ represents the residual force vector of the design domain derived by nonlinear finite element analysis, and its value being equal to 0 means the balance of force. $\mathit{U}$ is a vector which contains the nodal displacement. V represents the volume fraction and $V^{*}$ represents the target volume fraction. $x_i$ is the ith element in the vector $\mathit{x}$, where $x_i=1$ stands for a solid element and $x_i=x_{min}$ stands for a void element. The minimum value of the design variable, $x_{min}$, is defined as a small positive value (such as 0.001). For compliance minimization problem, $\mathit{l}$ equals to the external force vector, ${\mathit{F}}_{ext}$. As for the displacement maximization problem, $\mathit{l}$ represents a vector with its element corresponding to the direction of output displacement equals to −1.

3. Self-Adaptive Material Interpolation Scheme

Some compliant mechanisms, such as bi-stable mechanism and constant force mechanism, need to experience large deformations during operation. Thus, geometrical nonlinearity need to be considered in finite element analysis during topology optimization. However, topology optimization based on density interpolation often fails due to excessive distortion of some elements.

The main reason is that low-density elements are relatively softer and easier to deform. Therefore, even if the external force of the design domain is small, the low-density elements has a tendency of distorting excessively. Then, their tangent stiffness matrices may turn to indefinite or even negative definite. As a result, the displacement vector can not be calculated and the finite element analysis are force to stop. Many researchers [25,26,27] have investigated the above numerical instability phenomenon and found that it occurs in the region of low-density elements (e.g., those with relative densities between 0.1 and 0.3).

Here, we use the C-shape cantilever beam to demonstrate the numerical instability problem, while the schema of C-shape cantilever beam is presented in Figure 1a. The C-shape beam is fixed by its left side and two concentrate forces of different direction is applied on two nodes. The values of forces are set as $f_1=0.2$ N and $f_2=0.3$ N. The overall structure is meshed by square elements, and the side length of elements is 1mm. The Young’s modulus of elements is given as $E_1={10}^8$ Pa and the Possion ratio is 0.3. The gray area, ${\Omega }_{void}$, represents the void element region. The void region is empty and has no material, but the void region is filled with void elements in the density-based approach. Without the influence of void elements that have minimum density, the deformation of the C-shape cantilever beam is regarded as the true deformation, which is illustrated in Figure 1b.

In the topology optimization using the density-based method, the Solid Isotropic Microstructure with Penalty (SIMP) approach is widely used, which provides an interpolation method for defining the physical properties of the intermediate density element. In this way, the elastic modulus of the element can be formulated as a function of elemental relative densities, which is shown below:

$$E_e=x_e^p{E}_0$$

(2)

where $E_0$ represents the elastic modulus of the solid element. p is used to penalized the element density value, and its value is set to be 3 by referring to previous studies [2,21]. The minimum density value of element, $x_{min}$, is 0.001. Therefore, by using the interpolation scheme in Equation (2), the void element has an elastic modulus of ${10}^{-9}{E}_0$. A small elastic modulus can alleviate the influence of void element on predicting the behavior of the overall structure. However, if such elements appear in the early stage of optimization, it may lead to high distortion that would interrupt the finite element analysis. Therefore, void elements with a larger modulus of elasticity values are preferred in some cases.

Under this circumstance, a self-adaptive material interpolation scheme is proposed, aiming to improve the optimization stability while ensuring the accuracy. In the proposed scheme, the expression of the elastic modulus is written as:

$$E_e=E_0\left[x_e^p\left(1-\alpha \right)+\alpha \right]$$

(3)

where $\alpha$ is a control parameter with its value raging between 0 and 1. According to Equation (3), changing the value of $\alpha$ would not affect the elastic modulus of solid elements. Both the minimum-density elements and intermediate-density elements would be affected by the value of $\alpha$. By increasing (or decreasing) the value of $\alpha$, the elastic modulus of the minimum-density and intermediate-density elements will increase (or decrease) simultaneously. In order to derive a solution with fewer intermediate density elements, the elements’ density would be penalized by a power law function at the same time. The comparison between the interpolation schemes of Equations (2) and (3) is shown in Figure 2.

As shown in Figure 2, the value of $\alpha$ has a influence on the elastic modulus of the minimum-density and intermediate-density elements at the same time. The elastic modulus of the minimum-density elements is proportional to the value of $\alpha$. The proposed interpolation scheme in Equation (3) degenerates into the classical interpolation scheme shown in Equation (2), when the parameter $\alpha$ is set to 0. When the value of parameter $\alpha$ is increased, the element would become stiffer, and the deformation of the elements can be alleviated, which promotes the smooth convergence of the finite element analysis; however, it will inevitably leave an impact on the precision of the finite element analysis results. Decreasing the value of parameter $\alpha$ would increase the precision of the finite element analysis, but reduce the stability of the optimization process.

To demonstrate the effect of $\alpha$, a case analysis for the C-shape beam is given. By applying the proposed interpolation scheme and setting different $\alpha$, the finite element analysis results are compared. Firstly, elements with density ${\rho }_e=0.001$ are uniformly added to the void region, ${\Omega }_{void}$, of the C-shape cantilever beam. It should be noticed that, using the material interpolation scheme in Equation (2), the nonlinear finite element analysis would fail. Applying the proposed material interpolation scheme, the values of $\alpha$ are set to $1\times {10}^{-2}$, $1\times {10}^{-4}$ and $1\times {10}^{-7}$ in different cases, respectively. The effect of $\alpha$ on the finite element analysis can be visualized through the nonlinear finite element analysis of the above cases. Figure 3a–c show the deformation of the C-shape cantilever beam and deformed meshes with $\alpha$ values of $1\times {10}^{-2}$, $1\times {10}^{-4}$ and $1\times {10}^{-7}$, respectively.

Figure 3a shows the deformation of the C-shape cantilever beam when the parameter $\alpha$ equals $1\times {10}^{-2}$. The deformation presented in Figure 3a has a significant difference with respect to the true deformation shown in Figure 1b. The displacement direction of the beam at the bottom is downward rather than upward, as shown in Figure 1b. Therefore, an improper value of parameter $\alpha$ can significantly affect the finite element analysis result.

As presented in Figure 3a,b, the deformation of the C-shape cantilever beam becomes more obvious and also closer to the true deformation when the value of $\alpha$ decreases to $1\times {10}^{-4}$ and $1\times {10}^{-7}$. However, when the value of parameter $\alpha$ continues to decrease to $1\times {10}^{-8}$ or below, the finite element analysis, however, cannot converge and obtain an effective result because of the excessive deformation. In conclusion, decreasing the value of parameter $\alpha$ can improve the finite element analysis accuracy, but also increase the risk of failing of its failure.

To quantitatively compare the influence of various $\alpha$ values on accuracy, the output displacement is determined by measuring the displacement of the upper right corner of the C-shape beam. By setting the parameter $\alpha$ from $1\times {10}^{-7}$ to $1\times {10}^{-1}$, the corresponding estimated values of output displacement are summarized in Figure 4. The true value of output displacement, $U_{out}$, is calculated after removing the void elements totally from the void region, which is also known as element removal method. The comparison between estimated value and true value of output displacement is shown in Figure 4.

As shown in Figure 4, when $\alpha$ is in a certain range ($1\times {10}^{-7}\sim 1\times {10}^{-5}$ ), the estimated results of $U_{out}$ are very close to its true value. Therefore, choosing a suitable value of $\alpha$ can not only benefit the nonlinear finite element analysis process, but also ensure enough accuracy. In order to achieve a harmonious equilibrium between the precision and stability of the finite element analysis in the realm of topology optimization, a strategy of adaptive adjustment is put forth: an adaptive adjustment strategy based on maximum strain, which is proposed and shown below:

$${\alpha }^{\left(k\right)}=\left\{\begin{array}{cc}max\left({\epsilon }_{\mathrm{von}}^{(k-1)}/{\epsilon }^{*},0.8\right)·{\alpha }^{(k-1)},\hfill & \mathrm{if}{\epsilon }_{\mathrm{von}}^{(k-1)}\le {\epsilon }^{*}\hfill \\ \left({\epsilon }_{\mathrm{von}}^{(k-1)}/{\epsilon }^{*}\right)·{\alpha }^{(k-1)},\hfill & \mathrm{if}{\epsilon }_{\mathrm{von}}^{(k-1)}>{\epsilon }^{*}\hfill \end{array}\right.$$

(4)

where ${\alpha }^{\left(k\right)}$ represents the value of $\alpha$ at the kth iteration. ${\epsilon }_{\mathrm{von}}^{\left(k\right)}$ denotes the maximum von’s strain in all elements within the design domain at k iterations. $\alpha$ is conservatively defined as a relatively large value at the beginning of the topology optimization, thus ensuring that the finite element analysis converges properly. If the strain of element in previous iteration is smaller than a given value, ${\epsilon }^{*}$, then it means that the deformation is relatively slight at this time, so the value of $\alpha$ will be decreased to improve the accuracy. If the strain of element in previous iteration is greater than ${\epsilon }^{*}$, then it means that the deformation is relatively large, so the value of $\alpha$ will be increase to avoid the risk of excessive distortion. It can be found that when ${\epsilon }^{*}$ is set to 0.6, satisfactory topology optimization results can be derived in most cases. A detailed discussion can be found in Appendix A.

4. Finite Element Analysis

4.1. Hyperelastic Material Model

In this paper, large deformations of the structure are considered, so the linear assumptions used in the small deformation problem is not able to predict the true behavior of structure. By using the total Lagrangian finite element formulation, the Lagrange strain tensor of the structure is written as:

$$\mathit{E}=\frac{1}{2}\left({\mathit{F}}^T·\mathit{F}-\mathit{I}\right)$$

(5)

where $\mathit{F}=\mathit{I}+\mathit{u}\otimes \nabla$ represents the deformation gradient tensor. $\mathit{I}$ represents the identity tensor, with components described by the Kronekor delta symbol, ${\delta }_{ik}$.

The hyperelastic material model is suitable for simulating materials such as rubber and polymers that exhibit an elastic response when subjected to large deformations. The hyperelastic material model is defined by the strain energy potential equation, which is

$$W\equiv \Phi \left(\mathit{C}\right)+\varphi \left(J\right)$$

(6)

where $\mathit{C}={\mathit{F}}^T·\mathit{F}$ represents the right Cauchy–Green deformation tensor and J denotes the determinant of deformation tensor, i.e., $J=det\left(F\right)$.

Several hyperelastic models, such as Neo-Hookean model, Mooney–Rivlin model, Ogden model and Yeoh model, have been proposed and used in different fields of engineering. In this paper, the Neo-Hookean model is applied, since it has a simple form [28]. For Neo-Hookean material [29], the components of Equation (6) are shown as below:

$$\begin{array}{c}\hfill \Phi \left(\mathit{C}\right)=\frac{\mu }{2}({\overline{I}}_1-3)\\ \hfill \varphi \left(J\right)=\frac{K}{2}{(J-1)}^2\end{array}$$

(7)

where ${\overline{I}}_1$ is the first deviatoric strain invariant defined as ${\overline{I}}_1=J^{\frac{-2}{3}}{I}_1$ and $I_1=trace\left(\mathit{C}\right)$. $\mu$ is the initial shear moduli of the Neo-Hookean material, while the K represents the initial bulk moduli. The definitions of them are shown below:

$$\begin{array}{c}\hfill \mu =\frac{E}{2(1+\nu )}\\ \hfill K=\frac{E}{3(1-2\nu )}\end{array}$$

(8)

where E and $\nu$ represents the elastic modulus and Poisson’s ratio, respectively.

Taking the derivative of Equation (6) with respect to $\mathit{C}$, the 2nd Piola Kirchhoff stress can be derived as:

$$\begin{array}{cc}\hfill \mathit{S}& =2\frac{\partial W}{\partial \mathit{C}}=2\frac{\partial \Phi }{\partial \mathit{C}}+\frac{\partial \varphi }{\partial J}J{\mathit{C}}^{-1}\hfill \\ \hfill & =\mu J^{-2/3}\left(I-\frac{I_C}{3}{C}^{-1}\right)+KJ(J-1)C^{-1}\hfill \end{array}$$

(9)

According to Holzapfel [30], the Cauchy stress is written as:

$$\begin{array}{cc}\hfill \mathbf{\sigma }& =\frac{1}{J}\tau =\frac{1}{J}\mathit{F}\mathit{S}{\mathit{F}}^T\hfill \\ \hfill & =\frac{1}{J}\left(\mu J^{-2/3}\left(\mathit{b}-\frac{I_1}{3}\mathit{I}\right)+KJ(J-1)\mathit{I}\right)\hfill \\ \hfill & =\frac{1}{J}\mu \left(\overline{\mathit{b}}-\frac{{\overline{I}}_1}{3}I\right)+K(J-1)\mathit{I}\hfill \end{array}$$

(10)

where $\tau$ represents the Kirchhoff stress and $\overline{\mathit{b}}=J^{\frac{-2}{3}}\mathit{b}$.

The left Cauchy Green deformation tensor $\mathit{b}$ is written as:

$$\mathit{b}=\mathit{F}{\mathit{F}}^T$$

(11)

The material stiffness tensor is

$$\mathcal{C}=2\frac{\partial \mathit{S}}{\partial \mathit{C}}=4\frac{{\partial }^2\Phi }{\partial C\partial C}+2\frac{\partial }{\partial C}\left(\frac{\partial \varphi }{\partial J}\right)\otimes J{\mathit{C}}^{-1}+2\frac{\partial \varphi }{\partial J}\frac{\partial }{\partial \mathit{C}}\left(J{\mathit{C}}^{-1}\right)$$

(12)

4.2. Finite Element Formulation

Since the design domain is meshed in a series of elements, the displacement within the ith element is specified by interpolating between n nodal values:

$${\mathit{u}}_i=\sum _{a=1}^n{N}^a\left(\mathbf{\xi }\right)u_i^a$$

(13)

where $\mathbf{\xi }$ denotes the local element coordinates of an arbitrary point. $u_i^a$ is the displacement vector at node a. $N^a\left(\mathbf{\xi }\right)$ denotes the shape functions described in local element coordinates. n is the number of nodes of element.

Then, deformation corresponding to a given displacement field can be given:

$$F_{ij}={\delta }_{ij}+\frac{\partial u_i}{\partial {\xi }_j}={\delta }_{ij}+\sum _{a=1}^n\frac{\partial N^a}{\partial {\xi }_j}{u}_i^a$$

(14)

The relation between strain and displacement is written as

$$\mathit{\epsilon }=\mathit{B}{\left(\mathit{u}\right)}^T\mathit{u}$$

(15)

where $\mathit{B}$ represents the strain–displacement matrix.

By using the finite element method, the design domain is meshed. The Newton–Raphson method is employed to determine a displacement vector, denoted as $\mathit{U}$, that satisfies the structure equilibrium, which is

$$\mathit{R}(\mathit{x},\mathit{U})={\mathit{P}}_{ext}-{\mathit{P}}_{int}=0$$

(16)

where $\mathit{R}(\mathit{x},\mathit{U})$ represents the residual force vector. The internal nodal force vector ${\mathit{P}}_{int}$ is defined as ${\mathit{P}}_{int}={\sum }_e{\mathit{P}}_{int}^e$. The definition of elemental internal force vector ${\mathit{P}}_{int}^e$ is shown below:

$${\mathit{F}}_{int}^e=\mathit{B}{\left({\mathit{u}}_e\right)}^T{\mathit{S}}_e\mathrm{d}V$$

(17)

where ${\mathit{S}}_e$ is the 2nd Piola Kirchhoff stress of eth element and $\mathrm{d}V$ is the volume of element.

A state of static equilibrium is deemed to have been achieved when the norm of the residual vector $\mathit{R}(\mathit{x},\mathit{U})$ is below the threshold of an small positive quantity. The Newton–Raphson procedure is employed to find a solution of this nonlinear problem based on the following iterative calculations for linear systems:

$$\begin{array}{c}\hfill {\mathit{K}}_T^{\left(k\right)}\Delta {\mathit{u}}^{\left(k\right)}={\mathit{R}}^{\left(k\right)}\\ \hfill \Delta {\mathit{u}}^{\left(k\right)}={\mathit{u}}^{\left(k+1\right)}-{\mathit{u}}^{\left(k\right)}\end{array}$$

(18)

where $\left(k\right)$ denotes the number of iteration steps.

5. Sensitivity Analysis

Here, the adjoint method is employed for calculation of sensitivity analysis. An augmented objective function is constructed by adding a zero term to the objective function in Equation (1), which is written as:

$$\mathit{C}={\mathit{l}}^T\mathit{U}+{\mathit{\lambda }}^T\mathit{R}(\mathit{x},\mathit{U})$$

(19)

where $\mathit{\lambda }$ is the adjoint vector. Since the term $\mathit{R}(\mathit{x},\mathit{U})$ is a $\mathbf{0}$ vector, the adjoint vector $\mathit{\lambda }$ can be any random vector.

Differentiating both sides of Equation (19) yields

$$\begin{array}{cc}\hfill \frac{\partial \mathit{C}}{\partial x_e}& ={\mathit{l}}^T\frac{\partial \mathit{U}}{\partial x_e}-{\lambda }^T\left(\frac{\partial \mathit{R}}{\partial \mathit{U}}\frac{\partial \mathit{U}}{\partial x_e}+\frac{\partial \mathit{R}}{\partial x_e}\right)\hfill \\ \hfill & =\left({\mathit{l}}^T-{\lambda }^T\frac{\partial \mathit{R}}{\partial \mathit{U}}\right)\frac{\partial \mathit{U}}{x_e}+{\lambda }^T\frac{\partial \mathit{R}}{\partial x_e}\hfill \end{array}$$

(20)

Based on the premise that $\mathit{\lambda }$ is an arbitrary vector, we assume that ${\mathit{l}}^T-{\mathit{\lambda }}^T\frac{\partial \mathit{R}}{\partial \mathit{U}}=\mathbf{0}$. Then, the term involving with $\frac{\partial \mathit{U}}{x_e}$ can be eliminated. Therefore, the value of $\mathit{\lambda }$ can be given as follows [21]:

$${\mathit{\lambda }}^T=-{\mathit{l}}^T{\left(\frac{\partial \mathit{R}}{\partial \mathit{U}}\right)}^{-1}={\mathit{l}}^T{\mathit{K}}_T^{-1}$$

(21)

where ${\mathit{K}}_T^{-1}$ is the inverse of global tangent stiffness matrix.

Then, the formulation of sensitivity number is written as:

$$\begin{array}{cc}\hfill \frac{\partial \mathit{C}}{\partial x_e}& =-{\mathit{\lambda }}^T\frac{\partial {\mathit{F}}_{int}}{\partial x_e}\hfill \end{array}$$

(22)

According to Equation (3), the material stiffness tensor of element is related to $x_e$ by

$${\mathcal{C}}_e={\mathcal{C}}_0[x_e^p(1-\alpha )+\alpha ]$$

(23)

where p is defined as 3. And the parameter $\alpha$ determine the stiffness of minimum-density and intermediate-density elements. Thereby, the sensitivity number of the eth element is written as:

$$\frac{\partial \mathit{C}}{\partial x_e}=-{\mathit{\lambda }}^T\frac{\partial {\mathit{F}}_{int}}{\partial x_e}=-\frac{p{x}_e^{p-1}(1-\alpha )}{x_e^p(1-\alpha )+\alpha }{\mathit{\lambda }}_e^T{\mathit{F}}_e^{int}$$

(24)

where ${\mathit{\lambda }}_e^T$ and ${\mathit{F}}_e^{int}$ represents the adjoint vector and internal force vector of the eth element. The formulation of the sensitivity number is simple, which also makes it easy to implement.

Checkerboard and mesh-dependency problems are typical problems in topology optimization applications. An effective solution is to introduce a low pass filter to smooth the elemental sensitivity numbers. The weighted average scheme is employed as a sensitivity number filter, which is shown as follows [2,10]:

$$\widehat{\frac{\partial \mathit{C}}{\partial x_e}}=\frac{{\sum }_{m=1}^N{H}_{e,m}^S\frac{\partial \mathit{C}}{\partial x_m}}{{\sum }_{m=1}^N{H}_{e,m}^S}$$

(25)

Here, $H_{e,m}^S$ is a linear distance factor defined as:

$$H_{e,m}^S=max(0,r_s-D(e,m))$$

(26)

where $r_s$ is the sensitivity filter radius. $D(e,m)$ is the center-to-center distance between the eth and mth element.

The filtered elemental sensitivity numbers are averaged with the historical information to stabilize the evolutionary results, which is:

$${\left(\widehat{\frac{\partial \mathit{C}}{\partial x_e}}\right)}^k\leftarrow \frac{{\left(\widehat{\frac{\partial \mathit{C}}{\partial x_e}}\right)}^k+{\left(\widehat{\frac{\partial \mathit{C}}{\partial x_e}}\right)}^{k-1}}{2},\mathrm{for}k>2$$

(27)

where k denotes the number of iteration steps.

6. Numerical Examples

Numerical examples are provided for the purpose of validating the effectiveness of the proposed method in the optimization of geometrically nonlinear structures. Two of numerical examples are compliance minimization problems, which aim to improve the stiffness of structures. The last numerical example solves displacement maximization problem, which designs a compliant mechanism to achieve anticipated motion. In all these examples, the solid material Young’s modulus is $E_0=3.0$ GPa and the Poisson’s ratio is $\nu =0.4$. The optimization processes will be terminated when the changes in objective function during a series of iterations are sufficiently small, and the optimization stopping criterion is shown below:

$$\beta =\frac{\left|{\displaystyle \sum _{k-4}^k}{J}_i-{\displaystyle \sum _{k-9}^{k-5}}{J}_i\right|}{{\displaystyle \sum _{k-4}^k}{J}_i}<{\beta }^{*}$$

(28)

where $\beta$ represents the variation in the objective function and $J_i$ represents the value of the objective function in the ith iteration. ${\beta }^{*}$ denotes the allowable tolerance with typical values ranging between 0.001 and 0.05.

6.1. Topology Optimization of a Cantilever Structure

The problem of optimizing the classical cantilever structure is presented in the first numerical example. The design domain of cantilever is a rectangular with a size of 0.12 m × 0.03 m, which is shown as Figure 5. The fixed support of the rectangular design domain is on its left side, and the external force, F, is exerted on the middle part of its right side [21]. To relax the stress concentration, the external load is evenly applied on several nodes. The design domain is meshed into 3600 (120 × 30) finite elements and the target volume fraction of the available material is set to be $V^{*}=0.5$. At the start of the optimization, the design domain is filled with elements whose density are set to be 0.5. The filter radius for the element density and the sensitivity numbers are set as $r=0.0015$ m.

When the external force is small, the optimization can be accomplished by the original interpolation scheme shown in Equation (2). By setting $\alpha =0$, the proposed interpolation scheme is equivalent to the interpolation scheme shown in Equation (2). The optimized structures derived by different interpolation schemes are presented in Figure 6.

As shown in Figure 6, both the derived topology optimization results have up–down asymmetric structures, which indicates that the geometrical nonlinearity are considered. The optimized structure in Figure 6a is are derived with no means to help convergence. When the proposed method using the interpolation scheme in Equation (3), the obtained optimization result is presented in Figure 6b, which differs from the result shown in Figure 6a. The compliance of structures shown in Figure 6 are $4.139\times {10}^{-3} \mathrm{J}$ and $4.058\times {10}^{-3} \mathrm{J}$, respectively. In this example, it can be seen that, by using the proposed self-adaptive material interpolation model, a new structure can be derived and its compliance is smaller.

However, when the external force increase to 5N, the optimization process would fail due to the numerical instability if the original interpolation scheme in Equation (2) is applied. The finite element analysis could not produce a result because some elements were highly distorted. As a comparison, the geometrically nonlinear topology optimization based on the proposed material interpolation scheme can successfully produce an optimization result even with larger input forces. Figure 7 shows the optimization results for input forces of 5 N, 10 N, 15 N and 20 N, respectively.

As presented in Figure 7, the external force can be increased to 20 N while the numerical instability would not happen in the optimization. Also, a typical bar structure appears on the right part of the cantilever and will become larger as the input force increases, which is similar to the findings of other researchers [21,22,23].

The structural compliance convergence histories of the optimization results in Figure 7 are shown in Figure 8.

As can be seen from Figure 8, the compliance values of all the cantilever optimization examples converge stably.

In order to visualize how the self-adaptive material interpolation scheme works during procedure of topology optimization, the variation histories of parameter $\alpha$ are recorded and shown in Figure 9. As shown in Figure 9, when the input force is small (e.g., $F=5$ N and $F=10$ N), the deformation of all the elements is small, so the value of $\alpha$ shows a monotonically decreasing trend and eventually stabilizes near a value close to 0. When the input force is larger (e.g., $F=15$ N and $F=20$ N), the deformation of some elements may exceed the given maximum strain value ${\epsilon }^{*}$ during the optimization. Therefore, from the change curve of $\alpha$, it can be seen that there is a phenomenon that its value increases first and then decreases. But, at the end of optimization, the final values of $\alpha$ are reduced to a smaller value than the initial value. Therefore, from Figure 9, it is proven that the adaptive adjustment strategy shown in Equation (4) can stabilize the optimization convergence.

Table 1. Comparison of optimization results of cantilever beams under different loads.

Input Force (N)	Initial Value of $\mathit{\alpha }$	Final Value of $\mathit{\alpha }$	Displacement at Loading Point ($\mathsf{\mu }$m)	Compliance (N $·$ m)	Displacement by Element Removal ($\mathsf{\mu }$m)	Relative Error
5	0.02	$1.00\times {10}^{-6}$	−20.31	0.102	−20.31	0.01%
10	0.02	$4.07\times {10}^{-4}$	−35.98	0.360	−35.87	0.31%
15	0.02	$1.80\times {10}^{-3}$	−47.78	0.717	−47.39	0.81%
20	0.02	$2.20\times {10}^{-3}$	−55.95	1.119	−55.22	1.32%

shows the detailed optimization results. In each numerical examples, the final values of parameter $\alpha$ are different. The initial value of $\alpha$ of all the numerical examples are set to 0.02 and the final value increases with the value of the input force, which indicates that the larger the input force is, the harder the optimization converges.

As discussed earlier, the value of $\alpha$ has an impact on the precision of the finite element analysis and, therefore, the optimization results may be affected. To quantify the effect of $\alpha$, elements with density less than 0.1 were removed from the design domain. To exterminate the true deformation of the structure, an additional finite element analysis was performed on the optimization results. By comparing the displacement at the loading point before and after element removal, the relative error caused by the parameter $\alpha$ can be obtained. As can be seen from Table 1, the relative error also increases with input force. In the numerical examples discussed, the relative error is greatest when the input force is 20 N, with a value of $1.32\%$. From an engineering point of view, the error caused by introducing $\alpha$ is acceptable.

In summary, the proposed self-adaptive interpolation scheme can effectively address the numerical instability problem without decreasing the precision of finite element analysis. Moreover, the proposed method only causes very little computational stress and has no significant impact on the efficiency of the optimization. To illustrate this conclusion, the time consumption for the optimization iterations of the numerical examples shown in Figure 6 is analyzed. In each iteration, the optimization code can be partitioned into two primary parts, which are the finite element part and the numerical optimization part. Taking the first 100 iterations as an example, the average time spent in each part is listed in

Table 2. Comparison of time consumption.

	With Proposed Method	Finite Element Analysis (s)	Optimization (s)	Total Time (s)
Figure 6a	No	5.99	1.52	7.51
Figure 6b	Yes	5.98	1.85	7.83

.

Table 2 shows that, the time spent in the finite element analysis is almost the same, which takes up most of the time in each iteration. The introduction of the parameter $\alpha$ only affects the efficiency of the optimization part, since there is additional code for determining its value. In the case shown in Table 2, the code elapsed time for the optimization part increases from 1.52 s to 1.85 s. The increase of 0.33 s takes only $4.39\%$ of the total time; therefore, by the proposed method, the nonlinear topology optimization is able to be applied to the case of larger input force without sacrificing optimization efficiency. Meanwhile, the geometrically nonlinearity is considered during the topology optimization. Therefore, it is verified that the proposed self-adaptive interpolation scheme is capable of addressing the geometrically nonlinear topology optimization problem.

6.2. Topology Optimization of a Double Clamped Structure

The next numerical example is to optimize the double clamped beam, and its design domain is illustrated in Figure 10. The design domain of double clamped beam is a rectangular with a size of 0.16 m × 0.04 m. The input force F is exerted on the middle of the upper side of the design domain, while the left side and the right side are defined as fixed supports. The design domain is meshed into 160 × 40 four-node square elements. Because of the symmetric design, the optimization procedure considers only half of the domain, i.e., 3200 elements.

Also, $4\times 4$ solid elements are arranged on the location where the external loading is applied. By referring to the work of other researchers, the target volume, $V^{*}$, was set to 0.25. At the start of the optimization, the design domain is filled with elements whose density are set to be $x_e=0.25\left(e=1,2,...,3200\right)$. In this optimization example, the compliance of the structure is to be minimized. The radius of the filter is 0.0018 m, while the initial value of parameter $\alpha$ is 0.2. Four different optimization cases with different values of force are studied in this section. The force values are specified as $F=1.5$ N, $F=5$ N, $F=12.5$ N and $F=15$ N, and Figure 11a–d shows the optimization results.

As can be seen from Figure 11, the optimized topologies of the double clamped structures are similar. All optimized topologies are composed of two bars connected to the fixed support and a triangular structure directly subjected to the load. As the load F increases, the optimized topology only has small changes, such as the top angle of the triangle structure becomes smaller. Triangular structures with smaller top angles have greater stiffness, which helps to resist larger input forces.

The detailed optimized results are shown in

Table 3. Optimization results for the double clamped structure example by the self-adaptive interpolation scheme.

Input Force (N)	Initial Value of $\mathit{\alpha }$	Final Value of $\mathit{\alpha }$	Displacement at Loading Point ($\mathsf{\mu }$m)	Compliance (N $\mathit{·}$ m)	Displacement by Element Removal ($\mathsf{\mu }$m)	Relative Error
1.5	0.03	$1.00\times {10}^{-6}$	−1.141	0.102	−1.144	0.27%
5	0.03	$1.00\times {10}^{-6}$	−3.858	0.360	−3.870	0.29%
12.5	0.03	$1.00\times {10}^{-6}$	−8.452	0.717	−8.472	0.24%
15	0.03	$1.00\times {10}^{-6}$	−10.19	1.119	−10.22	0.33%

. The true finite element analysis of the optimized topologies are also calculated by removing elements with density less than 0.3 and the true displacements are presented in Table 3. Table 3 shows that the error caused by the parameter $\alpha$ is still very small and can be neglected.

6.3. Topology Optimization of a Displacement Inverter

One of a typical application of topology optimization is designing compliant mechanism. The third numerical example is to optimize a compliant displacement inverter, aiming at maximizing its output. Since the displacement inverter has a symmetric structure, the topology optimization is conducted on half of the displacement inverter to save computational resources. Then, the design domain can be simplified by eliminating the symmetric part and adding a symmetrical constraint. Figure 12 shows the sketch of the design domain as well as its boundary condition. The external load is applied on the upper left corner of the design domain, which is distributed uniformly over $4\times 4$ elements at the input port to avoid the excessive distortion. The density of elements in the loaded region is set to be 1 and will not alter during the whole optimization procedure. Moreover, $10\times 4$ solid elements are arranged on the corner where the fixed support is applied. The density of the other elements is set to 0.25 at the start of the optimization. The stiffness of input and output spring are defined as $K_{in}=500\mathrm{N}/\mathrm{m}$ and $K_{in}=100\mathrm{N}/\mathrm{m}$, respectively. The dimension of simplified design domain is 0.1 m × 0.05 m and is meshed into 5000 $\left(100\times 50\right)$ finite elements. The target material volume fraction is set to 0.25.

Figure 13 shows the optimization results of compliant displacement inverter under different loads with the $\alpha$ initialized as 0.03. The input loads corresponding to optimized results in Figure 13a–d are F = 1 N, F = 2.5 N, F = 5 N and F = 7.5 N, respectively. As can be seen in Figure 13, the bar structure connected to the input region gradually evolves into a hinge structure as the input load increases. The emergence of the hinge structure is reasonable since the inverter is undergoing a large deformation and the hinge structure can provide better load transfer by reducing the input stiffness.

The convergence histories are shown in Figure 14a. As can be seen in Figure 14a, the output displacement values converge stably for all numerical examples. The evolution process of the parameter, $\alpha$, for all the examples are shown in Figure 14b. With the proposed regulation scheme, the parameter $\alpha$ decreases rapidly and finally remains at a very low level $\left(\le 1\times {10}^{-4}\right)$, which ensures that the behavior of the optimization results are close to its actual behavior. The optimization results, including the output displacement before and after void element removal, and iteration numbers, are summarized in

Table 4. Optimization results for the displacement inverter example with self-adaptive material interpolation scheme.

Input Force (N)	Output Displacement ($\mathsf{\mu }$m)	Output Displacement after Element Removal ($\mathsf{\mu }$m)	Relative Error
1	−1.751	−1.752	0.07%
2.5	−4.281	−4.284	0.08%
5	−8.318	−8.394	0.90%
7.5	−12.29	−12.44	1.18%

. Since an increase in the input force tends to increase the deformation of structure, the value of $\alpha$ has a tendency to increase with the input force to prevent excessive deformation. As discussed previously, the value of $\alpha$ would affect the accuracy of the finite element analysis. In all cases above, the maximum value of $\alpha$ is $6\times {10}^{-5}$ and the relative error does not exceed $1.18\%$. This discrepancy is deemed satisfactory when considering the perspective of engineering.

6.4. Comparison

There were currently only three available educational codes of topology optimization considering geometrically nonlinearity. Chen et al. [21] introduced a 213-line code using the SIMP method, which was based on MATLAB R2016a and ANSYS 17.1. Zhu et al. [31] proposed a 89-line code, which was established on the open-source program platform FreeFEM. A 137-line code was built by Han et al. [32] based on the ESO method. The latter two codes do not take into account the effects of numerical instability, and therefore optimizing the numerical examples in Section 6 using these two codes, a failure to converge will be found. Based on the additive hyperelasticity technique [22] and SIMP method, the 213-line Matlab code has a better convergence performance. In this section, the optimization results of the 213-line code and the self-adaptive material interpolation scheme are compared. Taking the cantilever in Section 6.1 as an example, both the method are used to optimize the cantilever with loads ranging from 5 N to 25 N. The results of the comparison are presented in

Table 5. The results of the comparison between the self-adaptive material interpolation scheme and 213-line Matlab code.

Method	Compliance ($\mathbf{N}·\mathbf{m}$)					Time (s)
	5N	10 N	15 N	20 N	25 N	5 N	10 N	15 N	20 N	25 N
proposed method	0.102	0.360	0.717	1.119	1.568	24.679	32.946	33.404	36.150	40.370
213-line Matlab code	0.109	0.369	0.730	1.129	1.591	9.469	13.301	13.220	14.322	16.238

.

Compared to 213-lines code, the compliance of the optimization results derived using the self-adaptive material interpolation scheme is smaller. In addition, the optimization takes much less time by using the proposed method. The main reason is that the additive hyperelasticity technique used in 213-lines code requires adding another layer of elements to the design domain. More elements reduce the efficiency of finite element analysis as well as optimization. From the comparison, it can be seen that by using an adaptive material interpolation scheme, the proposed method can produces better optimization result efficiently.

7. Conclusions

A self-adaptive material interpolation scheme for topology optimization of geometrically nonlinear structures is proposed in this study, which is aimed at avoiding the numerical instability resulting from the excessive distortion of low-density elements. The proposed scheme is based on the famous SIMP approach, and a parameter, $\alpha$, is added to the material interpolation function to adjust the material properties of minimum-density and intermediate-density elements, which can change the resistance to deformation of these elements. The main idea behind this approach is to increase the stiffness of the low-density elements so that they do not deform excessively, thus ensuring that the finite element analysis can converge. However, the accuracy of the finite element analysis decreases. Considering the contradictory relationship between the approximation accuracy and convergence performance of the finite element analysis, a self-adaptive strategy is suggested to regulate the parameters.

Based on the proposed scheme, three numerical examples, including two minimum compliance and one compliant mechanism deign problems, are then formulated to validate the proposed scheme. The optimization results show that it is possible to obtain a well-defined solution with the proposed scheme, which is capable of withstanding greater loads and deformations. The obtained solutions are reasonable from an engineering point of view, and the error caused by introduction of parameter $\alpha$ is negligible, since the value of parameter $\alpha$ is decreased to a very small value by the adaptive adjustment strategy. From the analysis of the optimization results, the self-adaptive material interpolation scheme is proved to be effective. Moreover, by comparing the time consumption before and after applying the proposed scheme, it is found that the efficiency of optimization is not significantly affected.