基于SIMP法的板壳结构支撑位置拓扑优化设计

Topology Optimization Design of Plate and Shell Structure Support Position Based on SIMP Method

摘要: 拓扑优化技术通常在具有预定支撑条件的设计域上执行，以生成最优的结构。这种预定支撑条件的设计限制了优化设计空间，提高了使用成本。然而，现有的支撑位置优化技术是有限的，大多数方法需要繁琐的程序来预先定义支持条件。本研究提出了一种基于SIMP法的板壳结构支撑位置拓扑优化设计方法，通过在允许支撑的边界处引入一层单元，该方法可以简单地在有限元模型中实现。文中给出了各种实例来验证新方法的有效性。研究表明，将基于单元的支撑位置作为设计变量可以有效地获得高效和创新的结构设计。

Abstract: Topology optimization techniques are usually performed on a design domain with predetermined support conditions to generate the optimal structure. The design of this predetermined support condition limits the optimization design space and increases the use cost. However, the existing support position optimization techniques are limited, and most methods require cumbersome pro-cedures to pre-define support conditions. In this study, a topology optimization design method for the support position of plate and shell structures based on SIMP method is proposed. By introduc-ing a layer of elements at the boundary of the allowable support, this method can be simply imple-mented in the finite element model. Various examples are given to verify the effectiveness of the new method. This study shows that using the element-based support position as a design variable can effectively obtain efficient and innovative structural design.

1. 引言

结构拓扑优化是通过确定设计域中的材料分布来生成轻质、高性能和经济高效结构的有效方法之一。近年来发展出多种拓扑优化设计方法，包括均匀化方法 [1] 、密度法 [2] 和渐进结构优化方法 [3] 等。

在许多拓扑优化方法中，SIMP法(Solid Isotropic Material with Penalization，固体各向同性材料惩罚法)的简单性使得其在学术界和工程界广泛使用。SIMP法利用最优性准则在0和1之间连续改变单元伪密度并通过惩罚使单元伪密度趋近于0/1，其中1代表实体材料，0代表空隙。但是，SIMP法在求解过程中可能会出现数值不稳定现象，如出现棋盘效应、网格依赖。这时候可以通过使用过滤技术 [4] [5] 或者添加约束 [6] 来处理这些不利影响。在通常的拓扑优化问题中，结构的支撑位置是预先确定的，限制了拓扑优化设计空间。因此，支撑位置的优化可以作为一项单独的任务来实现，以提高结构的性能。

Jang等 [7] 通过将弹簧连接在结构单元所有节点上建立支撑单元，研究了大型板件在自重下的最优支撑位置，使得结构不均匀变形最小化。Buhlt [8] 以最大刚度为目标通过在结构单元节点上设置弹簧模拟支撑件建立了支撑优化区域并在对结构优化的同时获得了最优的支撑位置。Zhu [9] 等采用伪密度变量技术，将离散变量问题转化为连续变量问题，研究了与边界条件优化相对应的支撑布局设计，以最大化结构的固有频率。然而，这些方法对最终用户来说可能并不容易。它们要求明确定义系统刚度矩阵和支撑刚度矩阵，以确定结构系统的平衡并进行灵敏度分析。

本文提出了一种新的基于单元的优化方法来寻找支撑位置。新方法将一层具有边界条件的支撑单元引入到初始设计域，这可以在有限元模型中简单实现，并通过伪密度变量技术将支撑单元设计变量0/1分布，能够简单高效地得到具有最优支撑位置的板壳结构。

2. 设计方法

2.1. 问题描述

如图1所示为悬臂梁示意图，其中图1(a)为几何模型，图1(b)为有限元模型。通过4节点四边形单元离散设计模型，有限元模型由支撑单元和结构单元组成。在本研究中，引入了一层具有边界条件的单元来表示支撑的位置，支撑单元为设计域用于优化结构的支撑位置，结构单元为非设计域。

(a) 几何模型 (b) 有限元模型

Figure 1. Cantilever beam iagram

图1. 悬臂梁示意图

2.2. 优化设计数学模型

基于SIMP法的板壳结构支撑位置拓扑优化设计的优化目标是最小化结构的柔度。以支撑单元的伪密度作为设计变量，在一定的体积约束下更新支撑单元的伪密度，获得最优的支撑位置，实现板壳结构最优的静力学性能。其数学模型如下：

$\begin{array}{l}\text{find}X\\ \min C\\ \text{s}\text{.t}.F=KU\\ g\left(X\right)=v-\eta v_0\le 0\\ 0<x_{\min }\le x_i\le 1\end{array}$ (1)

式中，X为支撑单元伪密度集合，C为结构总柔度，F表示外载荷，K表示结构整体刚度矩阵，U表示结构整体位移矩阵。v表示支撑单元总体积， $\eta$ 表示体积约束因子，v₀表示初始支撑单元总体积，x_min表示设计变量下限，x_i表示第i个支撑单元伪密度。上式中的体积约束公式可以理解为：

$N^{*}={\displaystyle \underset{i}{\sum }{x}_i}$ (2)

式中N^*表示预设的支撑单元数量。

2.3. 灵敏度分析

由于SIMP法是基于梯度进行拓扑优化的算法，有必要求解目标函数和约束函数对于设计变量的偏导数即灵敏度，它反映了设计变量的变化对系统的影响，常被用于指导设计和搜索方向。目标函数C对于支撑单元设计变量x_i的灵敏度为：

$\frac{\partial C}{\partial x_i}=U^{\text{T}}\frac{\partial K}{\partial x_i}U+2U^{\text{T}}K\frac{\partial U}{\partial x_i}$ (3)

静力学平衡方程两边对设计变量x_i求偏导得：

$\frac{\partial U}{\partial x_i}=-K^{-1}\frac{\partial K}{\partial x_i}U$ (4)

将式(4)代入式(3)可得目标函数C对于支撑单元设计变量x_i的灵敏度：

$\frac{\partial C}{\partial x_i}=-U^{\text{T}}\frac{\partial K}{\partial x_i}U$ (5)

在本研究中SIMP插值模型为：

$E_i=x_i^p{E}_0$ (6)

上式中，E_i表示第i个支撑单元的弹性模量，E₀表示材料弹性模量。p为惩罚因子，在本文中取3。

式(5)可改为：

$\frac{\partial C}{\partial x_i}=-p{x}_i^{p-1}{u}_i^{\text{T}}{k}_0{u}_i$ (7)

上式中，u_i表示第i个支撑单元的位移矩阵，k₀表示其单元刚度矩阵。

3. 算例

为了验证本文方法的有效性，对不同边界条件和不同载荷条件的悬臂梁结构进行拓扑优化设计，分别是施加集中载荷、均布载荷和线性变化载荷的悬臂梁结构。结构材料的弹性模量为210 GPa，泊松比为0.3。

3.1. 施加集中载荷悬臂梁

设计模型如图2(a)所示，为一自由边顶点施加集中载荷的悬臂梁结构，其厚度/长度等于0.001，N^* = 1。支撑单元优化设计的主要过程如图2(b)~(e)所示。优化迭代31步后，最优支撑位置如图2(e)所示，图中黄色单元为最优支撑单元，右侧色条表示支撑单元的伪密度大小。最终结构的柔度与初始结构的柔度之比C/C₀ = 0.021。图2(f)为迭代历程曲线，可以看出，随着迭代优化的进行，结构的总柔度不断减小并最终收敛。

为了进一步说明本文方法的有效性，将N^*分别设为2、3、4和5时，对板壳结构进行优化设计。具有最优支撑位置的板壳结构如图3(a)所示，可以发现，支撑位置总是分布在施加载荷的一侧，这样的结构拓扑形态可以在加载点和固定端之间形成最短传力路径。随着板壳结构支撑单元的增加，板壳结构总柔度呈现下降趋势。图3(b)是优化迭代历程曲线，当N^*取不同数值时，本文方法均能够快速收敛到最优解。

(a) 设计模型 (b) k = 1 C = 23.41 J (c) k = 10 C = 0.72 J (d) k = 20 C = 0.64 J (e) k = 31 C = 0.49 J (f) 迭代历程曲线

Figure 2. Design results of cantilever beam support position under concentrated load

图2. 施加集中载荷悬臂梁支撑位置设计结果

3.2. 施加线性变化载荷悬臂梁

厚度/长度等于0.001的悬臂梁结构如图4(a)所示，其自由端施加线性变化载荷，N^* = 1。支撑单元优化设计的主要过程如图4(b)~(e)所示，优化迭代34步收敛，具有最有支撑位置的板壳结构如图4(e)所示。图4(f)是支撑单元优化设计的迭代历程曲线，可以看出，随着迭代步的增加，板壳结构总柔度不断下降并最终收敛。最终结构的总柔度与初始结构的总柔度之比C/C₀ = 0.017。

N^* = 2 C = 0.31 J N^* = 3 C = 0.23 J

N^* = 4 C = 0.19 J N^* = 5 C = 0.17 J

(a) 支撑位置优化设计结果

N^* = 2 N^* = 3 N^* = 4 N^* = 5 (b) 迭代历程曲线

Figure 3. Different number of support positions corresponding to the design results

图3. 不同数量支撑位置对应设计结果

(a) 设计模型 (b) k = 1 C = 41.48 J (c) k = 10 C = 1.03 J (d) k = 19 C = 0.80 J (e) k = 34 C = 0.71 J (f) 迭代历程曲线

Figure 4. Design results of cantilever beam support position under linear variable load

图4. 施加线性变化载荷悬臂梁支撑位置设计结果

当N^*分别设为2、3、4和5时，对板壳结构进行优化设计，具有最优支撑位置的板壳结构如图5所示。施加线性变化载荷后，支撑位置会随着悬臂梁自由端载荷的分布出现非对称式排列，并且随着支撑单元数量的增加，具有最优支撑位置的板壳结构总柔度不断减小。

N^* = 2 C = 0.45 J N^* = 3 C = 0.35 J

N^* = 4 C = 0.29 J N^* = 5 C = 0.25 J

图5. 不同数量支撑位置对应设计结果

3.3. 施加均布载荷悬臂梁

自由端施加均布载荷的悬臂梁结构如图6(a)所示，其厚度/长度等于0.001，N^* = 2。支撑单元的主要优化设计过程如图6(b)~(e)所示，迭代至26步收敛，最优支撑位置的板壳结构如图6(e)所示。最终结构的柔度与初始结构的柔度之比C/C₀ = 0.081。图6(f)为迭代历程曲线，可以看出，随着迭代优化的进行，结构的总柔度不断减小并最终收敛。

(a) 设计模型 (b) k = 1 C = 2.98 J (c) k = 10 C = 0.94 J (d) k = 19 C = 0.26 J (e) k = 27 C = 0.24 J (f) 迭代历程曲线

Figure 6. Design results of cantilever beam support position under uniform load

图6. 施加均布载荷悬臂梁支撑位置设计结果

当N^*分别设为3、4、5和6时，具有最有支撑位置的板壳结构如图7所示，随着支撑单元数量的增加，板壳结构总柔度呈下降趋势。由于施加了对称的均布载荷，支撑位置的分布也是对称的。

N^* = 3 C = 0.21 J N^* = 4 C = 0.15 J

N^* = 5 C = 0.14 J N^* = 6 C = 0.12 J

图7. 不同数量支撑位置对应设计结果

4. 结论

本文提出一种基于SIMP法的板壳结构支撑位置拓扑优化设计方法，通过在板壳结构中引入一层支撑单元，采用伪密度变量技术将支撑单元设计变量0/1分布，能够简单高效地得到具有最优支撑位置的板壳结构，可以显著降低结构总柔度。本研究实现了不同载荷条件和不同边界条件下，板壳结构支撑位置的拓扑优化设计，优化设计可获得具有最优支撑位置和最大刚度的板壳结构。

参考文献