Topological design of continuum structures with global stress constraints considering self-weight loads

1.   Introduction

Topology optimization is a design technique in which the best material distribution in the design domain can be obtained subject to given constraints and prescribed load conditions ^[1,2,3]. Since the pioneering work of Bendsøe and Kikuchi ^[4], great progress has been made on topology optimization. A variety of topology optimization methods have been developed, such as the homogenization method ^[5], variable density method ^[6,7], level set method ^[8,9], evolutionary structural optimization method ^[10,11], moving morphable component approach ^[12] and isogeometric method ^[13].

Many studies on topology optimization are performed only considering applied load, and self-weight loads are rarely considered. In engineering practice, self-weight loads are the primary factor for large-scale structures. Therefore, it is necessary that topology optimization of large-scale structures is carried out considering self-weight loads. Topology optimization of continuum structures with self-weight loads was first proposed by Turteltaub et al. ^[14]. Bruneel et al. ^[15] noted some challenges such as the nonmonotonic behavior of the compliance, the possible unconstrained character of the optimum and the parasitic effect for low densities when the solid isotropic material with penalization (SIMP) model is used to carry out topology optimization of continuum structures including self-weight loads. When the ratio between the self-weight load and the stiffness in low densities is too large, the structure itself cannot support the self weight, resulting in large deformation of the region. Erratic intermediate density patterns appear in the optimal topology, and the phenomenon is called the parasitic effect. Ansola et al. ^[16] adopted the improved evolutionary structural optimization method and the modified sensitivity calculation strategy to perform topology optimization of continuum structures considering self-weight loads. Holmberg et al. ^[17] put forward a method for worst-case topology optimization of self-weight loaded structures with semidefinite programming. Xu et al. ^[18] proposed a method for the topological design of continuum structures, including body forces, using the guide-weight method. Jain et al. ^[19] studied the effect of self-weight on optimal configurations obtained by topology optimization when the self-weight and applied load are applied. Zhang et al. ^[20] put forward a method for improved multiobjective topology optimization of continuum structures including self-weight. Han et al. ^[21] and Yang et al. ^[22] adopted the evolutionary structural optimization method to carry out topology optimization of continuum structures considering self-weight loads. Novotny et al. ^[23] proposed a regularized formulation for topology optimization of structures subject to self-weight loading using the topological derivative method. In the aforementioned studies, compliance optimization was adopted to design continuum structures considering self-weight. The structures obtained by topology optimization with minimum compliance usually cannot meet the strength requirement. Therefore, stress constraints should be considered during the topological design of continuum structures considering self-weight so that the maximum stress is less than the allowable stress.

However, most works in the literature focus on the topological design of continuum structures with stress constraints without considering self-weight loads. Duysinx et al. ^[24] first proposed a topology optimization method of continuum structures with local stress constraints. Compared with the stiffness topology optimization, some difficulties occur in the stress-constrained topology optimization problem ^[25,26]. The first is the singularity problem appearing in the variable density method. There are high stress values in low-density elements, which make the optimization algorithm incapable of eliminating the elements. The singularity problem can be remedied by the ε-relaxation ^[27] method and the p-q approach ^[28]. The bi-directional evolutionary structural optimization (BESO) method can be adopted to avoid the singularity problem ^[29,30]. The second difficulty is the local nature of stress constraints, which lead to a large number of constrains and computationally demanding sensitivity evaluation. Global stress strategy such as the Kreisselmeier–Steinhauser (KS) ^[31] functions and the P method ^[32] are employed to reduce the computational cost. Long et al. ^[33] proposed a method for stress-constrained topology optimization of continuum structures under harmonic excitation using the P method. Zhan et al. ^[34] performed topological design of multi-material structures with a global stress constraint. Meng et al. ^[35] put forward the stabilizing control schemes for topology optimization of thermo-elastic structures with stress constraints. Han et al. ^[36] proposed a novel topology optimization method to design geometrical and material nonlinear structures with maximum von Mises stress constraints. The third difficulty is the highly nonlinear behavior when stress levels are affected by changing topologies, particularly in critical regions such as reentrant corners. The density filtering method ^[37,38] was adopted to alleviate the difficulty. Currently, self-weight loads have been seldom considered to perform topology optimization of continuum structures with stress constraints. Recently, dos Santos et al. ^[39] proposed an approach for topology optimization of structures with local stress constraints and self-weight loads based on the topological derivative method. The number of local stress constraints is large, which may lead to low computational efficiency.

To improve computational efficiency, a new approach for topology optimization of continuum structures with a global stress constraint subject to both self-weight loads and mechanical loads is proposed in this paper. The rational approximation of material properties is employed to describe the material distribution for overcoming the parasitic effect for low densities. The structure volume is developed as the objective function to be minimized. The local stress constraints for all elements are aggregated into a global stress constraint using the improved P-norm method. The element density is employed as an elemental scale factor in the improved P-norm method, which has a positive effect on convergence of the stress constrained topology optimization problem. A model for stress-constrained topology optimization considering self-weight loads is established. The projection filtering method is adopted to avoid the phenomenon of numerical instability, and the optimization problems are solved using the method of moving asymptotes.

The remainder of this paper is organized as follows. In Section 2, the optimization formulation for stress-constrained topology optimization of continuum structures considering self-weight loads is presented. In Section 3, the global stress constraint strategy is introduced. In Section 4, the sensitivity analysis is described. In Section 5, several numerical examples are presented to showcase the validity of the proposed method. Conclusions are stated in Section 6.

2.   Optimization formulation

2.1.   Material interpolation model

For the structural topology optimization problem considering self-weight load, the parasitic effect occurs in the low-density region when the solid isotropic material with penalization (SIMP) model is adopted. The ratio of the element self-weight to the element stiffness in such a low-density region approaches infinity. The low-density configuration is not enough to bear its gravity, resulting in excessive displacement of the region. This may cause the optimization problem to not converge, and many unreasonable gray elements occur in the region. When the rational approximation of material properties (RAMP) model is applied, the ratio of the element self-weight to the element stiffness in the low-density area does not become unbounded. The parasitic effect for low densities can be avoided. Therefore, the RAMP model ^[40] is used to describe the material distribution.

where $E$ is the elastic modulus of element $e$, ${\rho _e}$ is the element density, that is, the design variable of element $e$, ${E_0}$ is the elastic modulus of the solid material, and $k$ is the penalization factor.

2.2.   Projection filtering

The density filtering ^[41,42] method is adopted to avoid numerical instabilities such as checkerboard and mesh dependency. The filtered density ${\tilde \rho _e}$ can be obtained by the weighted average of the densities of the adjacent elements,

where ${{\mathit{\Omega}} _e}$ is the neighboring set of elements within the minimum allowable radius ${r_{\min }}$. ${w_j}$ is the weight function, which can be expressed as

where ${r_j}$ is the distance between the centroids of element $e$ and element $j$.

Assume that the weight factor ${w_j}$ is defined as 0 for other elements beyond the filter radius. The filtered element density can be written as

where $N$ is the total number of finite elements, and ${W_{ej}}$ is the modified weight factor and defined as

To reduce the intermediate densities in the topological configurations, the smoothed Heaviside projection ^[43] is used to modify the filtered density. The Heaviside function is formulated as

where ${\bar \rho _e}$ is the physical density, $\beta$ is the parameter that controls the steepness of the projection, and $\eta$ represents the threshold value.

2.3.   Topology optimization model

Considering the combined action of mechanical loads and self-weight loads, minimization of the volume fraction of the structure is used as the objective function, and the maximum stress is developed as the constraint to meet the strength requirements. The model for stress-constrained topology optimization of continuum structures considering self-weight loads is established as

where ${V_f}$ is the volume fraction of the structure, ${V_0}$ is the volume of the design domain, ${v_0}$ is the volume of elements filled with solid material, ${\sigma _{\max }}$ is the maximum von Mises stress, $\sigma *$ is the allowable stress value, ${\boldsymbol{K}}$ is the global stiffness matrix of the structure, ${\boldsymbol{U}}$ is the vector of the nodal displacement, ${{\boldsymbol{F}}^m}$ is the vector of the mechanical load, ${{\boldsymbol{F}}^g}$ is the vector of the self-weight load, and $N$ is the number of finite elements. ${\bar \rho _{\min }}$ is the minimum value of the element densities, which is set to be 10^-3 to avoid the singularity of the global stiffness matrix.

Square four-node elements are applied for structural finite element analysis, and the self-weight load vector of the element $e$ can be expressed as

where ${\rho _m}$ is the material density, and $g$ represents the gravitational acceleration.

3.   Global stress constraint strategy

For the plane stress problem, the element stress is expressed as

where ${{\boldsymbol{D}}_0}$ is the constitutive matrix, ${\boldsymbol{B}}$ is the strain-displacement matrix, and ${{\boldsymbol{u}}_e}$ represents the vector of the element displacement. ${\sigma _{11}}$, ${\sigma _{22}}$ and ${\sigma _{12}}$ are the stress components of the element $e$.

The von Mises stress of the element $e$ can be calculated as

where ${\boldsymbol{V}}$ is the auxiliary matrix and expressed as

To avoid the singularity phenomenon of the stress-constrained topology optimization problem, the stress relaxation method ^[44] is applied to define the stress of intermediate density elements. The relaxation stress of the $e$th element can be represented as

If the local stress constraints are applied for every element, the total number of local constraints is too large, and the calculation efficiency is very low. To overcome this problem, the improved P-norm method ^[34,45] is adopted to aggregate the local stress constraints for every element into a global stress constraint. In the improved P-norm method, the element density is employed as an elemental scale factor which has the positive effect on convergence of the stress constrained topology optimization problem. The maximum stress is calculated as

where $\gamma$ is constraint scaling coefficient, and $p$ represents the norm parameter. When the iteration step is n, $\gamma$ can be expressed as

where ${\alpha ^{(n)}}$ is the parameter that can adaptively update the coefficient $\gamma$ during the optimization process. ${\alpha ^{(1)}}$ and ${\gamma ^{(1)}}$ are set to 0.5 and 1, respectively.

4.   Sensitivity analysis

The sensitivity of the objective can be obtained by differentiating Eq (2.7) to a design variable

By taking the derivative of Eq (2.4), the term $\partial {\tilde \rho _e}/\partial {\rho _j}$ in Eq (4.1) can be obtained as

The term $\partial {\bar \rho _e}/\partial {\tilde \rho _e}$ in Eq (4.1) is computed by differentiating Equation (2.6).

Substituting Eqs (4.2) and (4.3) into Eq (4.1) yields

The sensitivity of the maximum von Mises stress is expressed as

The term $\partial {\tilde \sigma _{{\text{PN}}}}/\partial {\sigma ^e}$ in Eq (4.5) can be calculated as

The term $\partial {\sigma ^e}/\partial {\bar \rho _e}$ in Eq (4.5) is obtained by differentiating Eq (3.4).

The partial derivative of the element von Mises stress to the element stress can be expressed as

The partial derivative of element stress to physical density can be computed as

Taking the derivative of the finite element equilibrium equation yields

The term $\partial {\boldsymbol{U}}/\partial {\bar \rho _e}$ in Eq (4.10) can be calculated as

Substituting Eq (4.11) into Eq (4.9), the term $\partial {{\boldsymbol{\sigma }}_e}/\partial {\bar \rho _e}$ is rewritten as

where $\partial {{\boldsymbol{u}}_e}/\partial {\boldsymbol{U}}$ is a transformation of local element degrees of freedom into global degrees of freedom.

Substituting Eqs (4.8) and (4.11) into Eq (4.7), the term $\partial {\sigma ^e}/\partial {\bar \rho _e}$ is rewritten as

Let $\eta = \frac{{\partial {{\tilde \sigma }_{{\text{PN}}}}}}{{\partial {\sigma ^e}}} = \gamma {\left({\sum\limits_{e = 1}^N {{{\bar \rho }_e}{{\left({{\sigma ^e}} \right)}^p}} } \right)^{\frac{1}{p} - 1}}{\bar \rho _e}p{\left({{\sigma ^e}} \right)^{p - 1}}$. Then, the term $\partial {\sigma _{\max }}/\partial {\bar \rho _e}$ is simplified as

where ${\boldsymbol{\lambda }}$ is the adjoint vector and obtained as

The term $\partial {\boldsymbol{F}}_e^g/\partial {\bar \rho _e}$ is obtained by differentiating Equation (2.8).

The method of moving asymptotes (MMA) ^[46] is suitable for solving complex topology optimization problems, and the method has good robustness. In this paper, the MMA method is employed to solve the optimization problem of stress-constrained topology optimization considering self-weight loads with high nonlinearity.

5.   Numerical examples

In this section, the method for stress-based topology optimization considering self-weight load is demonstrated by two numerical examples. The material properties used are givens as follows: Young's modulus ${E_0}$ is 68.9 GPa, Poisson's ratio $\mu$ is 0.33, and the material density ${\rho _m}$ is 2.7 × 10³ kg/m³. The gravitational acceleration is 9.8 m/s². The minimum allowable radius ${r_{\min }}$ is set to 3.5 times the element size, and the parameter $p$ of the P-norm is 8. The initial value of the parameter $\beta$ is 1, which will double up to 16 every 10 every iterations. The initial value of the design variables is set to 0.85.

5.1.   L-shaped beam

The design domain of the L-shaped beam is shown in Figure 1. The upper edge of the design domain is fixed, and the input load ${{\boldsymbol{F}}^m}$ = 7 kN is applied at the upper corner of the right edge. The load is evenly distributed on six nodes, as shown in Figure 1, to avoid stress concentration in the applied load region. The allowable stress constraint limit is set to 50 MPa. The design domain of the L-shaped beam has the dimensions of L = 1 m and H = 1 m, and it is discretized into 100 × 100 square four-node elements.

First, topology optimization of the L-shaped beams without and with stress constraints considering the self-weight load are performed. In the optimization problem without stress constraints, minimization of the compliance is used as the objective subject to the volume constraint, and the volume allowable limit value is obtained by the results of stress-constrained topology optimization. The results of topology optimization of the L-shaped beams without and with stress constraints considering the self-weight load are shown in Figures 2 and 3, respectively. The performances of the obtained L-shaped beams are listed in Table 1.

For the L-shaped beam in Figure 2, one can find that the reentrant corner occurs in the optimal layout, and the phenomenon of stress concentration at the corner is obvious, as shown in Figure 2(b). The maximum stress value of the L-shaped beam is 73.803 MPa, which exceeds the allowable stress limit value. Compared with the results obtained by topology optimization without stress constraints considering self-weight loads, the layout of the L-shaped beam obtained by topology optimization with stress constraints considering self-weight loads is different, and the reentrant corner is removed. The maximum stress value of the L-shaped beam obtained by the proposed method is 50.033 MPa, and the maximum stress constraint is satisfied. The stress in the L-shaped beam obtained by the proposed method is more uniformly distributed. However, the compliance of the L-shaped beam obtained by the proposed method is larger than that of the L-shaped beam obtained by topology optimization without stress constraints considering self-weight loads, as shown in Table 1. This shows that the L-shaped beam obtained by the proposed method has higher strength but smaller stiffness. The same material volume is used in both cases. The compliance of the L-shaped beam obtained by compliance topology optimization is minimized to maximize structural stiffness. However, the stress levels of the L-shaped structure obtained by stress-constrained topology optimization considering the self-weight load can be controlled by eliminating the reentrant corner.

Second, to analyze the effect of the self-weight loads on the results of stress-constrained topology optimization, stress-constrained topology optimization of L-shaped beam without considering self-weight loads is carried out. The topology optimization results are shown in Figure 4. The reentrant corner is also removed, and the maximum stress constraint can be satisfied. Compared with the results of stress-constrained topology optimization without considering self-weight loads, the layout of the L-shaped beam obtained by the proposed method is different, and more materials are distributed in the region near the boundary constraint. More materials are used in the optimal configuration obtained by the proposed method, as shown in Table 1. Compared with that of the case without considering self-weight loads, the total load applied on the structure in the case considering self-weight loads is greater. Therefore, more materials are applied in the structure to satisfy the same stress constraints. The maximum stresses of the two L-shaped beams obtained by stress-constrained topology optimization without and with considering self-weight loads are 49.998 MPa and 50.033 MPa, respectively. The strengths of the of two L-shaped beams are nearly equal. The compliance of the L-shaped beam obtained by stress-constrained topology optimization considering self-weight loads is smaller than the other. This shows that the L-shaped beam obtained by the proposed method has better stiffness. Compared with the result obtained by topology optimization without considering self-weight loads, more materials are used in the L-shaped beam obtained by the proposed method to satisfy the same stress constraint limit subject to the mechanical load and self-weight loads. The L-shaped beam obtained by the proposed method has better stiffness because of more material usage. Moreover, the three iteration processes are stable, and the constraints can be satisfied, as shown in Figures 2(c), 3(c) and 4(c).

In order to study the effects of different norm parameter p values on the obtained L-shaped beam, stress-constrained topology optimization of L-shaped beam considering self-weight loads wase carried out with different norm parameter p values. A comparison of optimal results obtained by the proposed method with different norm parameter p values is given in Figure 5 and Table 2. Optimal layouts of L-shaped beams obtained by the proposed method with different norm parameter p are different, and the layout of L-shaped beam obtained by using a low p value is similar to that obtained by compliance topology optimization. As the norm parameter p is increased, the stress distribution of L-shaped beam is more uniform, and fewer materials are used to satisfy the stress constraint.

Finally, the effects of different mesh densities on the obtained L-shaped beam are investigated. Stress-based topology optimization of L-shaped beams considering self-weight loads are performed under different mesh densities. The optimal layout and stress distribution of L-shaped beams obtained by the proposed method under different mesh densities are given in Figure 6. One can find that the optimal layouts of L-shaped beams under different mesh densities are different, and the stress constraints can be satisfied in all the cases. As the mesh density is increased, more materials are also used in the obtained L-shaped beam, as shown in Table 3. As mentioned in ^[30], this is because the element stress is closer to point stress due to denser mesh.

5.2.   Cantilever beam

The design domain of the cantilever beam is shown in Figure 7. The design domain has a size of L = 2 m and H = 1 m. The left edge of the design domain is fixed, and the input load ${{\boldsymbol{F}}^m}$ = 50 kN is applied at the middle of the right edge. To avoid stress concentration in the applied load region, the load is evenly distributed on five nodes. The allowable stress constraint limit is set to 50 MPa. The design domain is discretized into 100 × 50 square four-node elements.

The optimal configurations, stress distribution and iteration histories for topology optimization of the cantilever beams with and without stress constraints considering self-weight loads are shown in Figures 8 and 9, respectively. The maximum stress of the cantilever beam obtained by topology optimization without stress constraints considering self-weight loads is 73.320 MPa, which exceeds the allowable stress limit value. The region with high stress is mainly concentrated in the left fixed boundary, and the stress distribution of the structure is extremely uneven, as shown in Figure 8.

Compared with the results of topology optimization without stress constraint considering self-weight loads, the optimal configuration of the cantilever beam obtained by the proposed method is different. The maximum stress of the cantilever beam obtained by the proposed method is 49.999 MPa, and stress constraints can be satisfied. The stress in the cantilever beam is more uniformly distributed. However, the compliance of the cantilever beam obtained by the proposed method is larger than the other, as shown in Table 4. It shows that the cantilever beam obtained by the proposed method has higher strength but smaller stiffness. The same material volume is used in both cases. The compliance of the cantilever beam obtained by compliance topology optimization is minimized to maximize structural stiffness. The stress levels of the cantilever structure obtained by the proposed method can be controlled by changing the layout.

Stress-constrained topology optimization of the cantilever beam without considering self-weight loads was performed. The results obtained by topology optimization without considering the self-weight load are shown in Figure 10. Compared with the results of stress-constrained topology optimization without considering self-weight loads, the cantilever beam by the proposed method is different, and more materials are distributed in the region near the boundary constraint. More materials are used in the optimal configuration obtained by the proposed method, as shown in Table 3. Compared to that of the case without considering self-weight loads, the total load applied on the structure considering the self-weight load is larger. Therefore, more materials are applied in the structure to satisfy the same stress constraints. The maximum stresses of the cantilever beams obtained by the proposed method and the other are 49.999 MPa and 50.026 MPa, respectively. This shows that the strength of the cantilever beam obtained by the proposed method is nearly equal to that obtained by the other method. However, the compliance of the cantilever beam obtained by the proposed method is smaller than that obtained by the other method. This shows that the cantilever beam obtained by the proposed method has better stiffness. Compared with the results obtained by topology optimization without considering self-weight loads, more materials are used in the cantilever beam obtained by the proposed method to satisfy the same stress constraint limit under the mechanical load and self-weight loads. The cantilever beam obtained by the proposed method has better stiffness because of more material usage. Moreover, the three iteration processes are stable, and the constraints can be satisfied, as shown in Figures 8(c), 9(c) and 10(c)).

To investigate the effects of different stress constraints on the obtained cantilever beam, stress-constrained topology optimization of the cantilever beam considering self-weight loads was performed under different stress constraint limits. Optimal results obtained by the proposed method under different constraint limits are given in Figure 11 and Table 5. One can find that the optimal layout of cantilever beams obtained by the proposed method under different constraint limits is similar, and the stress constraints can be satisfied in all the cases. As the stress limit value is decreased, more materials are used in the obtained cantilever beam to meet stress constraints.

6.   Conclusions

A method for topology optimization of continuum structures with a global stress constraint considering self-weight loads was proposed. The rational approximation of material properties was applied to describe the material distribution. The minimization of the structure volume was used as the objective function. The improved P-norm method was applied to aggregate the local stress constraints for all elements into a global stress constraint. The projection filtering method was adopted to modify the sensitivities of the objective function and constraints. Stress-constrained topology optimization of continuum structure with mechanical and self-weight loads was performed. The topological structures obtained by the proposed method can meet the requirement of static strength.

Compared with the results of compliance topology optimization considering self-weight loads, the optimal configurations of the structures obtained by the proposed method are quite different, and the stress constraints are satisfied. The structures obtained by the proposed method have better strength, and the stress is more uniformly distributed. Compared with the results of stress-constrained topology optimization without self-weight loads, the optimal configurations of the structures obtained by the proposed method are different, and more materials are used in the topological structures to satisfy stress constraints. The strengths of the two structures are nearly equal, but the structures obtained by the proposed method have better stiffness. This shows that the structures obtained by the proposed method have better comprehensive performance. As the norm parameter is increased, the stress distribution of the structures obtained by the proposed method is more uniform, and fewer materials are used. As the stress limit value is decreased and the mesh density is increased, more materials are used in the obtained structures to meet stress constraints. The proposed method is performed based on linear elastic finite element theory. How to model stress-constrained topology optimization with geometrical nonlinearities considering self-weight loads will be investigated in the near future.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was funded by the National Natural Science Foundation of China (Grant Nos. 52065019, 52165002), Natural Science Foundation of Jiangxi Province, China (Grant Nos. 20224BAB204041, 20202BAB204015), and Project funded by China Postdoctoral Science Foundation (2022M723189). We thank Professor Krister Svanberg for providing the source code of MMA.

Conflict of interest

The authors declare there is no conflict of interest.