利用均质化方法对周期性梁状结构进行微观结构拓扑优化

《CMES - Computer Modeling in Engineering and Sciences》：Microstructural Topology Optimization for Periodic Beam-Like Structures Using Homogenization Method

基于拓扑优化与均质化理论，设计梁微构型以实现轻量化并满足多变形需求。实验验证其张力-扭转特性及能量吸收性能。

　　### 1. 引言

在许多工程领域中，梁结构因其承载能力和广泛应用而成为重要的设计对象，包括航空工程和土木工程（如桥梁和建筑结构）。梁的宏观性能在很大程度上取决于其微观结构，因此如何找到适合设计要求的微观结构成为关键问题。与尺寸和形状优化相比，拓扑优化提供了更大的设计自由度，并且是实现所需结构性能的有效工具。在过去的四十年里，各种拓扑优化方法被提出，如均质化方法、固体力学与惩罚（SIMP）方法、双向进化结构优化（BESO）方法以及水平集方法（LSM）。其中，SIMP方法因其概念简单而被广泛使用。

Sigmund利用渐近均质化方法（AHM）来评估周期性单元格的等效性能，推动了极端材料特性的优化研究。这种方法被扩展用于设计高性能功能材料，如最大体积模量、改善导热性和优化复杂模量（粘弹性材料）。Zhai等人提出了一种热传导驱动的方法，用于控制微观结构的边界扩展，以确保功能梯度材料的性能连续性。各向同性多孔材料因其方向独立的特性在工程中具有显著的潜力，但自然界中各向同性多孔材料极为罕见。Xiang等人通过拓扑调整策略，从各向异性配置中推导出各向同性多孔结构。基于AHM的拓扑优化方法可以得到自然界难以获得的材料特性，包括零或负的热膨胀系数、螺旋特性以及负的泊松比。Long等人发现，将负泊松比材料与正泊松比材料结合可以增强微观结构的刚度。Wang等人进一步推动了机械隐身的设计，展示了通过战略组装各种微观结构，可以精确操控弹性响应，使被保护物体在周围环境中变得不可探测。

过去十年，两尺度结构优化一直是研究的热点，相较于单一尺度微观结构设计方法，它提供了更高的设计灵活性。Liu等人构建的两尺度优化框架中，AHM建立了宏观结构性能与微观单元格配置之间的联系。这种方法大大减少了计算负担，因为它将原本对大型异质结构的优化分析转换为对等效均质结构的分析。在该优化框架基础上，后续研究扩展到热弹性问题和动态问题，展示了其在多物理领域的广泛适用性。为了充分挖掘两尺度拓扑优化的能力，研究人员进行了系统性的研究，以设计适用于异质宏观结构的功能梯度微观结构。随着增材制造技术的加速发展，对拓扑优化解决方案的可制造性提出了更高的要求。为了与计算机辅助设计集成，学者们发展了结合移动可变形组件/杆方法的两尺度优化框架。在航空航天工程中，加筋壳体有广泛的应用。Zhou等人应用两尺度技术，以最大化网格加筋圆柱壳的临界屈曲载荷。通过实施p范数约束，Wu等人成功设计了骨启发的多孔结构。随着对非线性优化的兴趣增加，Jia等人建立了用于多孔架构的碰撞优化框架，该框架消除了对灵敏度分析的需求。通过使用去均质化技术，Groen等人在桌面计算平台上实现了高保真度的多尺度结构优化。

Giavotto梁理论，考虑到扭曲和剪切变形，已被广泛用于梁截面优化。考虑到梁柱连接，Grubits等人优化了钢I梁，考虑了几何和材料非线性影响。为了减少非均匀梁优化的计算成本，Liu等人开发了一种映射策略，通过从单一参考剖面推导出所有几何相似截面的截面特性，实现属性的外推。在结合AHM的两尺度优化中，宏观结构分析通常依赖于三维弹性理论。然而，使用固体单元进行结构分析时，计算成本大大增加，而精度并没有显著提高。此外，AHM受限于全方向周期性，不适用于宽度和厚度方向上尺度分离有限的梁结构。Yi等人建立了用于周期性梁的AHM，实现了具有指定弯曲刚度和极端扭转刚度的微观结构设计。然而，这种方法忽略了剪切效应。Xu和Qian提出了一种基于放松的Saint-Venant解的梁微观结构拓扑优化方法，通过与两种应变能比相关的约束来解决材料分离问题。然而，通过Saint-Venant解考虑剪切变形，而不是通过剪切刚度，使优化过程变得不直观，分析也更加复杂。

我们提出了一种结合完全耦合均质化理论的拓扑优化框架，用于周期性微观结构设计。与传统的基于AHM的方法不同，我们的方法采用Timoshenko梁理论进行宏观分析。这种方法解决了梁结构在三维均质化分析中的等效本构建模误差，同时通过宏观尺度的维度减少提高了计算效率。与其他基于梁理论的研究相比，我们的方法不仅考虑了拉伸、弯曲和扭转，还考虑了剪切变形。优化目标聚焦于极端拉伸（压缩）-扭转特性、最小临界屈曲载荷的最大化以及结构柔顺性的最小化。数值和实验结果验证了所提出的优化框架的有效性。值得注意的是，实验结果揭示了梁长度对扭转行为的影响以及优化的螺旋结构在能量吸收机制中的表现。该方法在精度和效率方面具有显著优势，对具有大量梁组件的实际工程结构的优化设计具有广阔前景。

### 2. 周期性Timoshenko梁的均质化方法

本节简要回顾周期性梁结构的均质化方法，用于预测其有效刚度。等效Timoshenko梁的广义本构关系可以表示为：

$$
[N_1, M_3, M_2, T_1, Q_2, Q_3] = [D_{11}, D_{12}, D_{13}, D_{14}, D_{15}, D_{16}][\varepsilon_1, \kappa_3, \kappa_2, \kappa_1, \gamma_{12}, \gamma_{13}]
$$

其中，$[N_1, M_3, M_2, T_1, Q_2, Q_3]$ 和 $[\varepsilon_1, \kappa_3, \kappa_2, \kappa_1, \gamma_{12}, \gamma_{13}]$ 分别表示梁的广义力和应变，包括拉伸、两个方向的弯曲、扭转以及两个方向的剪切；$D_{\alpha\beta}$ 是有效刚度矩阵 $D$ 的分量，其中 $\alpha, \beta = 1, 2, ..., 6$。为了得到上述关系，均质化方法考虑了以下单元格问题：

$$
\frac{\partial \sigma_{ij}^{\varepsilon[\alpha]}}{\partial x_j} = 0 \quad \text{在} \quad V
$$
$$
\sigma_{ij}^{\varepsilon[\alpha]}n_j = 0 \quad \text{在} \quad S_{np}
$$

其中，$i, j$ 和 $m, n$ 分别等于1、2、3，用于三维问题；$V$ 是单元格域；$S_{np}$ 和 $S_p$ 分别表示非周期面（或自由边界面）和周期面，如图1所示；一个带有上标 $\varepsilon$ 的变量表示具有微观振荡的实际场；$\sigma_{ij}^{\varepsilon[\alpha]}$ 是实际应力，其计算公式为：

$$
\sigma_{ij}^{\varepsilon[\alpha]} = C_{ijmn} \varepsilon_{mn}^{\varepsilon[\alpha]}
$$

通过这些处理，我们可以获得梁结构的有效刚度。

### 3. 拓扑优化公式

#### 3.1 材料插值方案

基于修改后的SIMP插值规则，梁结构中每个单元格 $e$ 的杨氏模量可以表示为：

$$
E_e = E_{min} + \rho_e^{~} (E_0 - E_{min})
$$

其中，$E_{min}$ 和 $E_0$ 分别表示空单元格和实单元格的杨氏模量；$\rho_e^{~}$ 是单元格的物理密度，且 $\rho_e^{~} \in [0, 1]$。在考虑剪切效应的情况下，$E_e$ 可以通过插值方案进行调整，以确保结构的稳定性和优化的可行性。

#### 3.2 优化问题与灵敏度分析

在考虑梁结构的宏观性能时，优化问题可以被定义为：

$$
\text{find } \rho = (\rho_1, \rho_2, ..., \rho_N)
$$
$$
\text{max or min } f(\rho)
$$
$$
\text{s.t. } v_f = \frac{\sum_{e=1}^N \rho_e^{~} v_e}{\sum_{e=1}^N v_e} \leq v_f^*
$$
$$
0 \leq \rho_e \leq 1, \quad e = 1, 2, ..., N
$$

其中，$f(\rho)$ 是优化目标函数，$v_f$ 是当前单元格的体积分数，$v_f^*$ 是体积分数的上限，$v_e$ 是单元格 $e$ 的体积，$N$ 是单元格的数量。在微观结构设计中，有两种主要方法：一种直接使用等效微观结构特性作为目标，另一种则评估包含异质微观结构的梁的结构性能作为设计目标。我们的工作系统地涵盖了这两种方法。对于微观结构的极端刚度优化，我们通过拉伸（压缩）-扭转的超材料设计展示了方法的有效性，其他极端刚度优化问题遵循类似的原则。对于以宏观性能为目标的微观结构设计，我们的框架集成了两个关键目标：刚度最大化，这是结构拓扑优化的主要目标；以及屈曲阻力增强，这是细长梁应用中的关键要求。

### 4. 优化结果与验证

#### 4.1 极端拉伸（压缩）-扭转刚度的优化

为了获得具有极端拉伸（压缩）-扭转刚度的微观结构，目标函数被选为 $f(\rho) = D_{14}$，在体积约束 $v_f^* = 25\%$ 的条件下进行优化。如图3所示，目标值在迭代初期迅速增加，然后在每次 $\xi$ 变化时保持收敛，而体积分数的变化相对稳定，最终满足目标值 $25\%$。优化后的微观结构由倾斜杆和方形环组成，表现出轴向拉伸（压缩）-扭转特性，其刚度达到 $D_{14} = 373.6$。此外，其他刚度值如 $D_{11} = 77$、$D_{22} = 5834.6$、$D_{44} = 2898.6$、$D_{55} = D_{66} = 11.9$ 也得到满足。

为了比较螺旋结构的性能，还测试了具有等质量和长度的方形管样品（与3个单元格的螺旋结构相对应）。通过力-位移曲线计算能量吸收（EA），比较结果如图6所示，螺旋结构的峰值力远低于具有类似能量吸收能力的方形管，表明这种具有拉伸（压缩）-扭转特性的螺旋结构具有强大的缓冲能力，适用于结构保护。

#### 4.2 宏观梁性能的优化

在本部分中，分别讨论了不同长度（由5个和20个单元格组成）和不同截面（矩形24×20和正方形24×24）的悬臂梁的临界屈曲载荷最大化问题。初始设计域如图7所示，其中外框厚度为2，是不可设计区域。对于矩形截面（24×20）的单元格，其不可设计区域的外框厚度同样为2。体积约束设置为 $v_f^* = 40\%$。

从表1可以看出，优化后的最小临界屈曲载荷随着结构长度的增加而减少，表明长梁更容易发生屈曲。此外，材料倾向于从表面中心移动到八个角，这是由于长梁中的剪切效应减少所致。对于不同截面的单元格，优化后的微观结构也有所不同。对于24×20的矩形截面单元格，上-下和前-后表面的优化拓扑不同，这是由于两个主平面的刚度不同所致。对于24×24的正方形截面单元格，由于两个主平面的刚度一致，因此四个表面的优化拓扑相同。

以20个单元格的正方形梁为例，通过均质化方法预测的最小弯曲刚度为 $D_{22} = 1.75 \times 10^4$，通过公式 $\frac{\pi^2 E I}{(2l)^2}$ 计算的悬臂梁的最小临界屈曲载荷为0.19，与表1中的值一致，这验证了计算方法在预测屈曲行为方面的准确性。

#### 4.3 柔顺性的最小化

在此案例中，我们考虑了20个单元格组成的悬臂梁结构，在右端施加剪切力 $F = 0.2$。体积约束设置为 $v_f^* = 40\%$。

如图8所示，当仅考虑体积分数约束时，材料主要分布在上表面和下表面，而中间没有材料。这种优化结果对于这种细长梁是合理的，因为其变形主要由弯曲引起，而弯曲应力在顶部和底部最大。然而，从工程应用的角度来看，这种材料分离是不合理的。因此，除了体积约束外，还引入了剪切刚度约束 $D_{66} \geq D_{66}^*$。当剪切刚度从10.1增加到35时，柔顺性从35.3、36.3增加到39.2，这表明优化后的结构虽然在刚度上有所损失，但材料连接得到了保障，进一步证明了剪切约束的必要性和有效性。

### 5. 结论

本研究提出了一种用于周期性梁微观结构设计的拓扑优化框架，该框架同时考虑了极端微观结构刚度和增强宏观结构性能，并引入了剪切刚度约束。通过该方法，系统地研究了拉伸（压缩）-扭转刚度、临界屈曲载荷和结构柔顺性的优化问题。实验结果揭示了梁长度对扭转行为和能量吸收特性的影响。值得注意的是，优化后的螺旋结构在抗冲击应用中表现出巨大的潜力。通过讨论临界屈曲载荷问题，揭示了梁长度和截面形状对优化结果的影响。剪切刚度主要控制材料在侧面上的分布，有效消除了细长梁在柔顺性优化中的材料分离问题。

### 6. 资助声明

本研究得到了国家自然科学基金（项目编号11902015）、航天科技集团第五研究院减速与着陆实验室开放基金（项目编号EDL19092138）以及教育部春晖计划（HZKY20220014）的支持。

### 7. 作者贡献

作者们确认了论文的贡献如下：概念化，贾娇；方法论，贾娇；软件，贾娇和贺欣；验证，贺欣和刘振晨；数据管理，贺欣；撰写初稿，贾娇；撰写和编辑，贾娇；可视化，贺欣；监督，贾娇；项目管理，贾娇和吴世清；资金获取，贾娇、刘振晨和吴世清。所有作者都审查了结果并批准了论文的最终版本。

### 8. 数据和材料的可获得性

支持本研究发现的数据可以从通讯作者贾娇处获得，经合理请求后提供。

### 9. 伦理批准

本研究不涉及伦理问题。

### 10. 利益冲突

作者们声明在本研究中没有需要报告的利益冲突。