Topological Optimization of Bi-Directional Progressive Structures with Dynamic Stress Constraints under Aperiodic Load

Abstract

The topology optimization of dynamic stress constraints is highly nonlinear and singular and has been little studied. Dynamic stress based on progressive structural optimization is only available by applying the modal iteration method, but due to the nonlinear limitations of the modal superposition method, there is an urgent need to develop a progressive structural optimization method based on dynamic stress sensitivity under direct integration. This method is for the dynamic stresses under non-periodic loading with iterative cycle updating variations. This article proposes a topological optimization method of continuum structures with stress constraints under an aperiodic load based on the Bi-directional Evolutionary Structural Optimization Method (BESO). First, the P-norm condensation function was used to obtain the global stress to approximate maximum stress. By introducing the Lagrange multiplier, the design goal was to increase the P-norm stress on the basis of the smallest volume. After that, based on the dynamic finite element theory, the sensitivity of each cell formula of the objective function and the constraint conditions of the design variables were strictly derived. Then, the performance evaluation index was put forward based on volume and stress, and the convergence criterion based on the performance evaluation index was defined. This method solves the topology optimization problem of stress constraints under a non-periodic load and the topology optimization problem of stress constraints under a periodic load, such as a simple harmonic load.

1. Introduction

With the increasing demand for lightweight structure and performance improvement in automotive, marine, aerospace, and other fields, structural optimization technology has been widely used in the past two decades. Topology optimization can lead to better performance designs in the conceptual design stage; so, many scholars and designers use this technique widely. In the traditional topology optimization field, the research mainly focuses on the minimization of structural flexibility [1,2,3]. Structural stiffness is increased by optimizing the material distribution. The topology optimization design technology based on structural flexibility minimization is mature, but topology optimization based on dynamic stress correlation is still in the development stage. At the same time, the stress of the structure under a static load and the stress under a dynamic load are very different, and the important criteria for measuring the strength failure of the structure are also different. Therefore, stress-constrained topology optimization is a hot spot in current research, but most of the current research focuses on static stress, and this paper mainly focuses on topology optimization under dynamic stress.

Compared with the minimum flexibility design, there are three challenging problems in the stress-constrained topological optimization problem [4], namely the singular phenomenon of stress, the localization of the stress constraint, and the highly nonlinear nature of stress behavior. Firstly, the stress singular phenomenon mainly appears in the method based on variable density, because the stress value of the low-density element always behaves as a high stress value; so, it cannot be completely eliminated by the optimization algorithm. Rozvany [5] proposed a smooth boundary function, namely the K-S (Kresselmeier Steinhauser) function. The smooth boundary function is used to relax the stress constraint to solve the stress singular value problem. Cheng Gengdong et al. [6] proposed an ε-relaxation method that can effectively search the global optimal solution in the feasible design domain. Matteo et al. [7,8] proposed a stress relaxation method based on the ε-relaxation method to deal with the stress singular value problem. Duysinx et al. [9] applied the ε-relaxation method to the topology optimization problem with local stress constraints.

Secondly, in the process of topology optimization, the position of maximum stress may shift due to the change in structural configuration in the different iterations. To constrain the maximum stress of the structure, the stress of each element must be constrained, and the optimization process needs to calculate the sensitivity of each constraint with respect to the density of the element, which causes a huge amount of calculation. To solve this problem, we can use a condensation function method, such as the P-norm method and the K-S function method. This method condenses multiple local stress constraints into a global stress function and only needs to constrain the global stress function to constrain the stress of the whole structure. Yang [10,11] used the K-S condensation function method to deal with the topology optimization problem of local stress constraints, condensing multiple local stress constraints into an overall stress function; then, it was only necessary to constrain this overall stress function to constrain the stress of the entire structure. This algorithm solves the topology optimization to minimize the volume of the structure under stress constraints. Rong [12] established a topology optimization method for a continuum structure with stress gradient constraints, solved the problem of stress concentration, and proposed a topology optimization algorithm through P-norm and Lagrange multipliers. Fan Zhao [13] and others proposed the following. Using the BESO method [14,15], they adopted the global stress measure based on P-norm to deal with the flexibility minimization problem under volume and stress constraints.

Thirdly, the highly nonlinear stress behavior and the stress value of any point on the structure are closely related to the topological configuration of the whole structure, especially in some areas with high stress concentration. The high stress distribution is extremely sensitive to subtle topological changes. The addition and deletion of elements will lead to strong stress constraint discontinuity, making the calculation convergence more difficult. There are intermediate density elements that appear in the topological structure, particularly in the framework of the variable density method; so, it is difficult to accurately calculate the stress of the boundary elements. To solve this problem, the current common solution is to use cell density filtering [16] to smooth the density of each cell and its surrounding cells, so as to reduce the instability of the optimization process.

When the structure is subjected to transient dynamic load or periodic dynamic load, the stress of the structure changes with time. Therefore, there is the problem of slow convergence, in addition to the above three problems. Due to these problems, there are few studies on the dynamic stress constraint under dynamic response [17]. Long [18] and others studied the topology optimization of stress constraints under simple harmonic loads based on the variable density method and the P-norm stress aggregation function. Pingyao [19] and others converted the dynamic stress response of the manipulator into a static load based on the equivalent static load method and combined this with the variable density method to carry out topological optimization. Lei Zhao [20,21] and others proposed the following. Based on the BESO method, they proposed a topological optimization method using a continuum structure constrained by a dynamic stress response under a random load and used a P-norm aggregation function instead of a dynamic stress response constraint. Lei Zhao [22] and others used the constraint limitation variational method to obtain a stable and convergent topology and derived the sensitivity of the dynamic fatigue constraint based on design variables to form an approximate function of the dynamic fatigue constraint function and objective function, thus realizing topology optimization.

The rest of the article is arranged as follows. The second section expounds the stress-constrained topology optimization of continuum structures under an aperiodic load. The third section derives the expression of element sensitivity based on the Newmark method. The fourth section defines the performance index by considering the changes in structural stress and overall volume, and it improves the convergence criterion. The fifth section gives several typical numerical examples to verify the effectiveness of this method. The sixth part summarizes the main conclusions of this article.

2. Topological Optimization Problem with Stress Constraints under Aperiodic Load

2.1. Transient Dynamics Analysis

In the actual engineering problems, the load form borne by the structure is often dynamic. The general motion control equation expressed in terms of time is as follows:

$$M\ddot{u}(t)+C\dot{u}(t)+Ku(t)=F(t)$$

(1)

where M is the structural mass matrix, C is the structural damping matrix, K is the structural stiffness matrix, $F(t)$ is the node load vector varying with time, and $\ddot{u}$, $\dot{u}$, and $u$ are the node acceleration, node velocity, and node displacement vector, respectively.

There are two methods to solve the dynamic finite element equations, namely the direct integration method and the modal superposition method. The basic idea of the direct integration method is to divide the time domain into a series of time points, which divides the time domain into several time intervals and then solves the solution that can satisfy Formula (1) at the discrete time points. The modal superposition method is to superimpose the responses of the system under various modes, but the proportion of the weight factors of each order is different. In this article, the birth–death element method is used to realize the addition and deletion of elements in the process of topological optimization. This process is a nonlinear analysis; so, it can only use the direct integration method but not the mode superposition method. The common direct integration methods include the Newmark method, Wilson method, Houbolt method, central difference method, etc. This section focuses on the Newmark method.

The Newmark method assumption:

$${\dot{u}}^{t+\Delta t}={\dot{u}}^t+\left[\left(1-\beta \right)\right]{\ddot{u}}^t+\beta {\ddot{u}}^{t+\Delta t}\Delta t$$

(2)

$$u^{\mathrm{t}+\Delta \mathrm{t}}=u^t+\Delta t{\dot{u}}^t+\left[\left(\frac{1}{2}-\gamma \right){\ddot{u}}^t+\gamma {\ddot{u}}^{t+\Delta t}\right]\Delta t^2$$

(3)

where $\gamma$ and $\beta$ are the Newmark integral parameters. When $\beta \ge 0.5$ and $\gamma \ge 0.25{\left(0.5+\beta \right)}^2$, the Newmark method has unconditional stability. $\ddot{u}$, $\dot{u}$, and $u$ represent, respectively, the node displacement vector, velocity vector, and acceleration vector at time t; $u^{t+\Delta t}$, ${\dot{u}}^{t+\Delta t}$, and ${\ddot{u}}^{t+\Delta t}$ represent, respectively, the node displacement vector, velocity vector, and acceleration vector at time $t+\Delta t$.

To calculate the structural displacement vector $u^{t+\Delta t}$, the equation of motion at time $t+\Delta t$ can be expressed as:

$$M{\ddot{u}}^{t+\Delta t}+C{\dot{u}}^{t+\Delta t}+K{u}^{t+\Delta t}=F^{t+\Delta t}$$

(4)

The effective stiffness matrix and effective load at time $t_{n+1}$ can be expressed as:

$$\overline{K}=K+\frac{1}{{\alpha }_0\Delta t^2}M+\frac{\beta }{{\alpha }_1\Delta t^2}C$$

(5)

$$\begin{array}{cc}\hfill {\overline{F}}^{t+\Delta t}=& \hfill F^{t+\Delta t}+M\left({\alpha }_0{u}^t+{\alpha }_2{\dot{u}}^t+{\alpha }_3{\ddot{u}}^t\right)+\\ \hfill & C\left({\alpha }_1{u}^t+{\alpha }_4{\dot{u}}^t+{\alpha }_5{\ddot{u}}^t\right)\hfill \end{array}$$

(6)

where ${\alpha }_0=\frac{1}{\gamma \Delta t^2}$, ${\alpha }_1=\frac{\beta }{\gamma \Delta t}$, ${\alpha }_2=\frac{1}{\gamma \Delta t}$, ${\alpha }_3=\frac{1}{2\gamma }-1$, ${\alpha }_4=\frac{\beta }{\gamma }-1$, ${\alpha }_5=\frac{\Delta t}{2}\left(\frac{\beta }{\gamma }-2\right)$, ${\alpha }_6=\Delta t\left(1-\beta \right)$, ${\alpha }_7=\beta \Delta t$.

The equation of displacement at time $t+\Delta t$ can be expressed as:

$$\overline{K}{u}^{t+\Delta t}={\overline{F}}^{t+\Delta t}$$

(7)

The velocity and acceleration at time $t+\Delta t$ can be expressed as:

$${\ddot{u}}^{t+\Delta t}={\alpha }_0\left(u^{t+\Delta t}-u^t\right)-{\alpha }_2{\dot{u}}^t-{\alpha }_3{\ddot{u}}^t$$

(8)

$${\dot{u}}^{t+\Delta t}={\dot{u}}^t+{\alpha }_6{\ddot{u}}^t+{\alpha }_7{\ddot{u}}^t$$

(9)

2.2. Topology Optimization Model Design

The stress level of the structure is an important criterion to measure whether the structure will have strength failure, and it is also a design factor that the structural designer must consider at the early stage of the design stage.

Based on the BESO method, the discrete topology optimization design variables can solve the stress singular value problem in essence, which has certain advantages for solving the stress-constrained topology optimization problem. With the minimum volume of the structure as the objective, the maximum stress as the constraint, and the presence or absence of materials as the design variable, the optimization model is as follows:

$$\begin{array}{l}\mathrm{find}:x=[x_1,x_2,\cdots ,x_n]\\ \min :V={\displaystyle \sum _{i=1}^n{V}_i{x}_i}\\ s.t.:M\ddot{U}+C\dot{U}+KU=F\\ {\sigma }_{vm}^{\max }\le {\sigma }_{vm}^{\ast }\\ x_i=0 \mathrm{or} 1\end{array}$$

(10)

Where $x$ is the design variable and $x_i$ represents the design variable value corresponding to the ith element. The design variables can be understood as the relative density of the element and, more directly, as an indicator of the presence or absence of the element material. When $x_i$ equals 1, it indicates that element i is a solid element with material. When $x_i$ equals 0, it indicates that cell i is an empty element without material. n is the total number of finite elements in the design domain, V is the overall volume of the structure, V_i is the volume of ith element, ${\sigma }_{vm}^{\max }$ is the maximum Mises stress of the structure, and ${\sigma }_{vm}^{\ast }$ is the artificial allowable stress.

The above optimization model means that the structure with the smallest volume can be obtained by updating the design variables under the conditions necessary to satisfy the given load, boundary conditions, and stress constraints.

2.3. Topology Optimization Model after Coalescing Stress Replacement

The design domain is divided into n finite element meshes, and each element corresponds to a design variable. The design variable can be described as:

$$x_i\in \{0,1\} ,i=1,2,\cdots n$$

(11)

The design variable can be understood as the relative density of the element, and more directly, it can be considered as an indicator of the presence or absence of the element material. When $x_i$ = 1, it indicates that the ith element is a solid element with material. Conversely, it indicates that the ith element is an empty element without material. The interpolation model between the design variables and the elastic matrix of the element material is established.

$$D_i=x_i{D}_0$$

(12)

where $D_i$ represents the material elastic matrix corresponding to the ith element.$D_0$ represents the elastic matrix of the solid materials, which is only related to the material properties and is a constant matrix.

According to the elastic mechanics theory, the stress at a point $a$ of the solid structure is as follows:

$${\sigma }_a={[{\sigma }_{ax},{\sigma }_{ay},{\sigma }_{az},{\tau }_{axy},{\tau }_{ayz},{\tau }_{azx}]}^T$$

(13)

In the finite element analysis, the following formula can be described as:

$${\sigma }_a=D{B}_{a}u$$

(14)

where $B_a$ is the strain matrix of point $a$, and $u$ is the displacement vector of point $a$.

Therefore, the equivalent stress model of the element can be described as:

$${\sigma }_i=\frac{D_i{B}_{i}u(x)}{x_i}=D_0{B}_i{u}_i$$

(15)

where ${\sigma }_i$ is the stress vector of the ith element, $B_i$ is the strain matrix of the centroid of the ith element, and $u_i(x)$ is the node displacement matrix of the ith element. The von Mises stress can be described as:

$${\sigma }_{vm}={({\sigma }^{T}V\sigma )}^{{}_{{}^{1/2}}}$$

(16)

where ${\sigma }_{vm}$ is the von Mises stress, $\sigma$ is the element stress vector, and V is the stress coefficient matrix.

Because the constraint condition in this optimization problem model is the maximum Mises stress, which has the local attribute of the stress constraint, its size cannot be described by a clear mathematical analytical formula, and the stress sensitivity of each element needs to be calculated in the optimization process. In order to reduce the calculation cost caused by the local stress constraint, the aggregation function is usually used to obtain a global stress measure to approximate the maximum Mises stress of the structure. In this article, the P-norm condensation function is adopted, and it can be described specifically as:

$${\sigma }_{PN}={({\displaystyle \sum _{j=1}^n{\sigma }_{vm,j}^p})}^{\frac{1}{p}}$$

(17)

where ${\sigma }_{PN}$ is the cohesion stress value of the structure and p is the stress norm value considered to be set; ${\sigma }_{vm,i}$ is the Mises stress value at the center of the ith element. From the expression form of this expression, when p = 1, ${\sigma }_{PN}$ is the sum of the Mises stress values of all the elements. When p approaches infinity, ${\sigma }_{PN}$ is the exact Mises stress value of the structure.

After the maximum Mises stress is replaced by the condensed stress, the new topological optimization mathematical model can be described as:

$$\begin{array}{l}\mathrm{find}:x_i=[x_1,x_2,\cdots ,x_n]\\ \min :V={\displaystyle \sum _{i=1}^n{V}_i{x}_i}\\ s.t.:M\ddot{U}+C\dot{U}+KU=F\\ {\sigma }_{PN}\le {\sigma }_{vm}^{\ast }\\ x_i=0 \mathrm{or} 1\end{array}$$

(18)

3. Cell Sensitivity Analysis

In order to solve the problem that the local properties of the constraints that appear in the optimization problem model cannot be described by a clear mathematical analytic formula, this chapter derives the element sensitivity and analyzes and calculates BESO as the basic method.

In order to solve this problem using BESO method, the Lagrange multiplier needs to be introduced, namely λ. The stress constraint and volume constraint are added to the objective function, while the BESO method takes the structural volume as the evolution parameter. So, the volume constraint is easy to meet. So, the objective function justly introduces a Lagrange multiplier, and the objective function can be described as:

$$f_1=V+\lambda ({\sigma }_{PN}-{\sigma }_{vm}^{\ast })$$

(19)

${\sigma }_{vm}^{\ast }$ is considered to be a given stress constraint value and is a constant value. Therefore, the derivative of this term is 0 during the sensitivity analysis. In sensitivity analysis, the objective function $f_1$ is equivalent to the following formula:

$$f=V+\lambda {\sigma }_{PN}$$

(20)

The element sensitivity calculation formula is derived. The partial derivative of the objective function $f$ to the design variable $x_i$ is obtained as follows:

$$\frac{\partial f}{\partial x_i}=\frac{\partial C}{\partial x_i}+\lambda \frac{\partial {\sigma }_{PN}}{\partial x_i}$$

(21)

According to the chain rule, it is easy to obtain ${\sigma }_{PN}$ with respect to the design variables in Formula (21). The sensitivity can be described as:

$$\frac{\partial {\sigma }_{PN}}{\partial x_i}=\frac{\partial {\sigma }_{PN}}{\partial {\sigma }_{vm,i}}\frac{\partial {\sigma }_{vm,i}}{\partial {\sigma }_i}\frac{\partial {\sigma }_i}{\partial x_i}$$

(22)

It can be seen from the above formula that to calculate the sensitivity of the global stress cohesion function, the derivative of the stress cohesion function to the Mises stress, the derivative of the Mises stress to the stress component, and the derivative of the stress component to the design variable must be calculated first. The following three items are calculated separately.

(1)

The derivative of the stress condensation function to the Mises stress

For the Mises stress of the P-norm in Formula (22) with respect to each element, we should strive to derive:

$$\frac{\partial {\sigma }_{PN}}{\partial {\sigma }_{vm,i}}={({\displaystyle \sum _{j=1}^n{\sigma }_{vm,j}^p})}^{(\frac{1}{p}-1)}·{({\sigma }_{vm,i})}^{p-1}$$

(23)

(2)

The derivative of the Mises stress to the stress component

The Mises stress in Formula (22) is derived from the stress vector of each element:

$$\frac{\partial {\sigma }_{vm,i}}{\partial {\sigma }_i}={\sigma }_{vm,i}^{-1}{\sigma }_i^{T}V$$

(24)

(3)

The derivative of the stress vector to the design variable

The stress vector in Formula (22) is derived from the design variable $x_i$, as follows:

$$\frac{\partial {\sigma }_i}{\partial x_i}=D_0{B}_i\frac{\partial u_i}{\partial x_i}$$

(25)

In the structural dynamics analysis, the stress produced varies with time. The dynamic stress is calculated according to the element, which is the same as the stress in the static finite element analysis. That is, the stress at the node of the element or the stress at any point on the element can be calculated through the element stress calculation formula. The dynamic sensitivity can also be calculated according to the element; that is, the element stress sensitivity can be obtained by calculating the partial derivative of the design variable with the element stress calculation formula. As the dynamic stress changes with time, the dynamic stress sensitivity also changes with time. If the solution time step is n in the dynamic finite element analysis, the stress at n time points can be obtained, and the stress sensitivity at n time points can also be obtained.

It can be seen from Formula (15) that the dynamic stress sensitivity calculation formula at time $t+\Delta t$ is:

$$\frac{\partial {\sigma }_i^{\mathrm{t}+\Delta t}}{\partial x_i}=D_{0}B\frac{\partial u_i^{t+\Delta t}}{\partial x_i}$$

(26)

where ${\sigma }_i^{\mathrm{t}+\Delta t}$ is element stress. It can be seen from the above formula that as long as the dynamic displacement sensitivity is obtained, the dynamic stress sensitivity can be obtained according to the above formula.

Based on the finite analysis, the displacement at each time can be obtained according to Equation (7); then, the velocity and acceleration at each time can be obtained according to Equations (8) and (9). The dynamic displacement sensitivity analysis is carried out according to the Newmark method.

The partial derivative of the design variable $x_i$ on both sides of Equation (7) is obtained as follows:

$$\frac{\partial \overline{K}}{\partial x_i}{u}^{t+\Delta t}+\overline{K}\frac{\partial u^{t+\Delta t}}{\partial x_i}=\frac{\partial {\overline{F}}^{t+\Delta t}}{\partial x_i}$$

(27)

So,

$$\overline{K}\frac{\partial u^{t+\Delta t}}{\partial x_i}=\frac{\partial {\overline{F}}^{t+\Delta t}}{\partial x_i}-\frac{\partial \overline{K}}{\partial x_i}{u}^{t+\Delta t}$$

(28)

wherein

$$\overline{K}=K+{\alpha }_{0}M+{\alpha }_{1}C$$

(29)

$$\frac{\partial \overline{K}}{\partial x_i}=\frac{\partial K}{\partial x_i}+{\alpha }_0\frac{\partial M}{\partial x_i}+{\alpha }_1\frac{\partial C}{\partial x_i}$$

(30)

$$\begin{array}{l}\frac{\partial {\overline{F}}^{t+\Delta t}}{\partial x_i}=\frac{\partial F^{t+\Delta t}}{\partial x_i}+\frac{\partial M}{\partial x_i}({\alpha }_0{u}^t+{\alpha }_2{\dot{u}}^t+{\alpha }_3{\ddot{u}}^t)+\frac{\partial C}{\partial x_i}({\alpha }_1{u}^t+{\alpha }_4{\dot{u}}^t+{\alpha }_5{\ddot{u}}^t)a\\ +M({\alpha }_0\frac{\partial u^t}{\partial x_i}+{\alpha }_2\frac{\partial {\dot{u}}^t}{\partial x_i}+{\alpha }_3\frac{\partial {\ddot{u}}^t}{\partial x_i})+C({\alpha }_1\frac{\partial u^t}{\partial x_i}+{\alpha }_4\frac{\partial {\dot{u}}^t}{\partial x_i}+{\alpha }_5\frac{\partial {\ddot{u}}^t}{\partial x_i})\end{array}$$

(31)

Equations (28)–(31) constitute the calculation formula of dynamic displacement sensitivity. The calculation formula of dynamic displacement sensitivity can be obtained by synthesizing (28)–(31):

$$\frac{\partial u^{t+\Delta t}}{\partial x_i}=(K+{\alpha }_{0}M+{\alpha }_{1}C{)}^{-1}\left[\begin{array}{l}\frac{\partial F^{t+\Delta t}}{\partial x_i}+\frac{\partial M}{\partial x_i}({\alpha }_0{u}^t+{\alpha }_2{\dot{u}}^t+{\alpha }_3{\ddot{u}}^t)\\ +\frac{\partial C}{\partial x_i}({\alpha }_1{u}^t+{\alpha }_4{\dot{u}}^t+{\alpha }_5{\ddot{u}}^t)\\ +M({\alpha }_0\frac{\partial u^t}{\partial x_i}+{\alpha }_2\frac{\partial {\dot{u}}^t}{\partial x_i}+{\alpha }_3\frac{\partial {\ddot{u}}^t}{\partial x_i})\\ +C({\alpha }_1\frac{\partial u^t}{\partial x_i}+{\alpha }_4\frac{\partial {\dot{u}}^t}{\partial x_i}+{\alpha }_5\frac{\partial {\ddot{u}}^t}{\partial x_i})\\ -(\frac{\partial K}{\partial x_i}+{\alpha }_0\frac{\partial M}{\partial x_i}+{\alpha }_1\frac{\partial C}{\partial x_i})u^{t+\Delta t}\end{array}\right]$$

(32)

When calculating the displacement sensitivity at time $t+\Delta t$, in addition to the displacement at time $t+\Delta t$, the displacement, velocity, and acceleration sensitivity at time t are also used. It can be seen from the formula that the dynamic displacement sensitivity analysis based on finite element analysis is used to solve the displacement sensitivity at the current time by solving the dynamic displacement response at the current time and combining the displacement, velocity, and acceleration sensitivity at the previous time. The displacement sensitivity of the whole time domain can be obtained by iteration. According to the initial conditions, it is easy to determine that the sensitivity of displacement, velocity, and acceleration at time t = 0 is 0; F is the external load vector, independent of the design variable $x_i$; so, the partial derivative of this term is 0. Therefore, the displacement sensitivity at time $\Delta t$ is:

$$\frac{\partial u^{\Delta t}}{\partial x_i}=-(K+{\alpha }_{0}M+{\alpha }_{1}C{)}^{-1}(\frac{\partial K}{\partial x_i}+{\alpha }_0\frac{\partial M}{\partial x_i}+{\alpha }_1\frac{\partial C}{\partial x_i})u_{{}^i}^{\Delta t}$$

(33)

Finally, the expression of dynamic stress sensitivity at time t can be obtained by introducing Equations (23), (24), (26) and (33) into Equation (22), as shown in Equation (34):

$$\frac{\partial f}{\partial x_i}=V_i-\lambda {{\displaystyle \sum _{j=1}^n\left[{({\sigma }_{vm,j}^{\Delta \mathrm{t}})}^p\right]}}^{(\frac{1}{p}-1)}{({\sigma }_{{}_{vm,i}}^{\Delta t})}^{p-2}{({\sigma }_i^{\Delta t})}^{T}V{D}_0{B}_i{L}_i{\overline{K}}^{-1}(\frac{\partial K}{\partial x_i}+{\alpha }_0\frac{\partial M}{\partial x_i}+{\alpha }_1\frac{\partial C}{\partial x_i})u_{{}^i}^{\Delta t}$$

(34)

When a uniform grid is used (i.e., the cells have the same volume), the relative ranking of the sensitivities of each cell can be defined by the following sensitivities:

$${\alpha }_i={{\displaystyle \sum _{j=1}^n\left[{({\sigma }_{vm,j}^{\Delta \mathrm{t}})}^p\right]}}^{(\frac{1}{p}-1)}{({\sigma }_{{}_{vm,i}}^{\Delta t})}^{p-2}{({\sigma }_i^{\Delta t})}^{T}V{D}_0{B}_i{L}_i{\overline{K}}^{-1}(\frac{\partial K}{\partial x_i}+{\alpha }_0\frac{\partial M}{\partial x_i}+{\alpha }_1\frac{\partial C}{\partial x_i})u_{{}^i}^{\Delta t}$$

(35)

In the whole process of structural optimization, a series of displacement values in the time domain will be output in detail after each finite element analysis. If these displacement values are processed in parallel, the whole calculation will be more complicated and difficult, and the time cost will increase. Moreover, in the optimization iteration step under a fixed mode determined by the structure, even if the displacement is always changing, the displacement is only related to the time history of the load; so, the series of displacement values is the same in the physical sense. Therefore, the element sensitivity at only one time point is sufficient to meet the requirements.

4. Convergence Criterion

The BESO method uses volume evolution to add and delete elements. It is obviously inappropriate to use the change in the objective function value as the convergence criterion. Therefore, based on the dimensionless concept, this section defines a performance index, which can comprehensively consider the changes in structural stress and overall volume at the same time. With this index, the changes in structural performance with the number of iterations can be reflected. The performance index is shown in Formula (36):

$$PI=\frac{{\sigma }_{PN}^0{V}_0}{{\sigma }_{PN}^t{V}_t}$$

(36)

where PI is the performance index; ${\sigma }_{PN}^0$ and ${\sigma }_{PN}^t$ represent, respectively, the initial and tth iteration stress cohesion function values; and $V_0$ and $V_t$ represent, respectively, the initial and tth iteration structure volumes.

Therefore, the termination criterion of the whole optimization process is as shown in Formula (37):

$$error=\frac{\left|{\displaystyle \sum _{t=1}^N\left(P{I}_{k-t+1}-P{I}_{k-N-t+1}\right)}\right|}{{\displaystyle \sum _{t=1}^{N}P{I}_{k-N-t+1}}}\le \tau$$

(37)

Generally, N is set to 5, which means that the change range of the performance index value in the past 10 iterations is small enough to reach convergence when it reaches less than the convergence factor.

5. Numerical Example

5.1. L-Shaped Two-Dimensional Plate

The zone size and boundary conditions of the L-shaped plate are shown in Figure 1. It is an L-shaped area with fixed constraints on the upper boundary. The design area is meshed by a four-node plane element with a side length of 2 mm. In order to avoid stress concentration at the load position, the load shown in Figure 2 is applied on the right boundary of the L-shaped beam. The stress constraint threshold is set to 4.5 MPa. In this article, the calculation example uses the same material parameters. Where Young’s modulus, Poisson’s ratio, and the material density are set as E = 71 Gpa, $v$= 0.33, and $\rho$= 2793 Kg/m³. With the BESO method, the material evolution rate er = 0.2%, the sensitivity filter radius is set to two times that of the cell grid, and the stress norm parameter p is set to 6.

According to the developed topology optimization program, the change curve of the structure volume fraction and the performance index is shown in Figure 3. It can be seen from Figure 4 that the structural volume decreases continuously with the optimization. When the iteration proceeds to step 90, the change tends to be gentle. After step 45, the performance index fluctuates violently; this is caused by the highly nonlinear stress. When the performance index reaches around 1.7, it tends to converge. At this time, the objective function value is 0.51.

The final topological structure and its stress cloud diagram are shown in Figure 4 and Figure 5. Due to the high nonlinearity of the dynamic stress, there are fluctuations in the whole optimization process. The cohesion stress reaches convergence after 105 iterations, and the final value is 4.33 Mpa. Xia developed the static stress topology optimization program based on BESO by using P-norm stress condensation and carried out the topology optimization of an L-shaped plate. Figure 6 shows the topology configuration diagram when P = 6. The same stress condensation parameters are used in this article. It can be seen from the comparison that the support bifurcation on the right side in the topology configuration in the literature makes the support structure on the right side smaller. In addition, the configuration is basically the same.

5.2. L-Plate with Notch

The initial design area and boundary conditions of the L-plate with a notch are shown in Figure 7. The design area is meshed by a four-node plane element with a side length of 2 mm. Fixed constraints are imposed on the left boundary of the L-plate. In order to avoid stress concentration at the load position, five finite element nodes at the lower vertex of the right boundary of the L-plate act as shown in Figure 2. The stress constraint value is set to 4.5 MPa.

Figure 8 shows the change curve of the structural volume fraction and performance index in the process of the structural optimization of the notched L-shaped beam. It can be seen from the figure that the structural volume decreases continuously with the optimization. When the iteration proceeds to step 60, the change tends to be gentle. The performance index fluctuates sharply after step 45, which is caused by the high nonlinearity of the stress. When the performance index reaches around 1.80, it tends to converge, and the objective function value is 0.44. Figure 9 shows the maximum stress value of the structure and the P-norm condensed stress value change curve in the optimization process. Due to the high nonlinearity of the stress, there are fluctuations in the whole optimization process. A total of 83 iterations were performed to achieve convergence.

The final topological structure and its stress cloud diagram are shown in Figure 10 and Figure 11. The maximum stress of the structure is 3.14 Mpa, and the P-norm cohesion stress is 4.00 Mpa. Long realized the topology optimization of the dynamic stress constraint under harmonic excitation by sequential quadratic programming. Figure 12 shows the topological configuration results of the same calculation example. Due to the different load time histories, the left support leg in the configuration of this article is obviously wider, and in addition, a fine structure is generated in the middle. The other topological configurations are basically the same.

5.3. Horizontal Z-Shaped Pipe

The initial design area and boundary conditions are shown in Figure 13. The design area is divided by a four-node plane element, finite element mesh with a side length of 2 mm. Fixed constraints are imposed on the left boundary of the Z-shaped beam. In order to avoid stress concentration at the load position, the loads shown in Figure 3 are applied to five finite element nodes at the right boundary vertex of the Z-shaped beam. The stress constraint value is set to 4.5 MPa.

Figure 14 shows the change curve of the structural volume fraction and performance index during the optimization of the Z-shaped beam structure. It can be seen from the figure that the structural volume decreases continuously with the optimization. When the iteration proceeds to step 155, the change tends to be gentle. The performance index fluctuates sharply after 50 steps; this is caused by the high nonlinearity of the stress. In Figure 14, when the performance index reaches around 1.68, it tends to converge, and the objective function value is 0.52. Figure 15 shows the change curve of the maximum stress value and P-norm condensed stress value of the structure during the optimization process. Due to the large size change in the structure, a certain degree of stress concentration will occur at the corner, and the dynamic stress will fluctuate greatly during the iteration, which makes it very difficult to converge. A total of 175 iterations were carried out. The final topological structure and its stress cloud diagram are shown in Figure 16 and Figure 17. The maximum stress of the structure is 4.07 Mpa, and the P-norm condensed stress is 6.29 Mpa.

6. Conclusions

There are three difficult problems in the topology optimization of stress constraints, and the dynamic stress constraints are highly nonlinear and difficult to converge. In this article, the topology optimization of dynamic stress constraints is attempted, and a topology optimization method using a continuum structure based on dynamic stress response constraints is proposed. A topological optimization model of structural volume minimization under dynamic stress response constraints is established. At the same time, the P-norm aggregation function is used to reduce the number of dynamic stress response constraints, and the stress sensitivity is derived based on the Newmark method. In addition, the P-normcondensation stress is used to improve the performance evaluation index based on volume and stress, and the convergence criterion based on the performance evaluation index is defined. Finally, through three examples of topology optimization constrained by minimum structural volume and dynamic stress response, the optimal topology obtained is reasonable and has good material distribution, which verifies the applicability and rationality of the dynamic stress response-constrained topology optimization method proposed in this article.

The stability and optimization efficiency of the method need to be further improved. The structures studied in this paper are all isotropic materials, all of which are in the linear elastic range. The application of the fatigue-constrained topology optimization method based on the use of the BESO method in anisotropic materials, functionally graded materials, and lattice structures needs to be further studied, so as to further expand its application scope.