Level-set topology optimisation with unfitted finite elements and automatic shape differentiation

技术摘要：基于非匹配有限元与自动形状微分的水平集拓扑优化

问题陈述
拓扑优化是一种关键的计算技术，在工业领域有着广泛的应用，传统上依赖于密度法或水平集法。传统的水平实现方法通常采用“平滑（smudging）”技术，即通过 Heaviside 函数在界面处对边界进行平滑处理，从而允许在静态背景域上进行积分。虽然这种方法计算效率高，但在非线性或流固耦合问题中可能会引入不理想的伪影，并且无法充分捕捉边界扰动下的泛函形状导数。

非匹配有限元方法（如 CutFEM）提供了一种替代方案，通过将复杂域嵌入到简单的包围域中，且不增加额外的自由度，从而提供对边界的精确描述。然而，这些方法面临着挑战：

1. 解析推导： 对非匹配离散化进行一阶和二阶形状导数的推导在数学上非常困难且易出错。
2. 边界交截： 现有的解析形状微积分结果（例如 Berggren [16]）通常假设设计域边界完全包含在背景域内，排除了设计域边界与背景边界相交的情况。
3. 孤立体积： 在非匹配方法中，不连通的设计域区域（孤立体积）会导致系统矩阵出现奇异性，因为它们不像密度方法中那样通过弱“伪（ersatz）”材料与 Dirichlet 边界相连。
4. 精度与通用性： 目前尚不清楚自动微分技术能否在非匹配离散化的背景下，将自动微分恢复到与解析形状导数相同的机器精度，特别是在处理复杂的形状海森矩阵（shape Hessians）时。

方法论
作者开发了一个结合非匹配有限元方法与自动形状微分的框架，该框架在 Julia 包 GridapTopOpt 及更广泛的 Gridap 生态系统中实现。

1. 扩展解析形状微积分：
   论文将现有的解析结果（Berggren [16]）扩展到了设计域边界 $\partial\Omega$ 与背景域边界 $\partial D$ 相交的情况。作者推导了在这些交截条件下，对于一般体积分、面积分以及通量积分的定向形状导数。他们证明，只要满足特定的网格交截假设，半导数（semiderivatives）即使在设计域边界接触背景边界时仍保持一致性。

2. 自动形状微分 (ASD)：
   其核心方法论贡献是将前向模式自动微分应用于非匹配离散化。

• 对偶数 (Dual Numbers)： 该方法利用对偶数 ($a + b\epsilon$) 在切割单元（cut cells）的数值积分规则中传播导数。
• 实现： 将水平集函数值在节点处替换为对偶数。切割单元的几何形状（由水平集零等值面与背景网格的交集定义）以及相关的雅可比矩阵均使用这些对偶数进行计算。
• 验证： 作者证明，无论网格尺寸如何，ASD 都能以机器精度恢复出解析定向导数。这消除了手动推导复杂网格相关量的需求。此外，该技术允许计算形状海森矩阵以及对于那些在 Berggren 框架内尚未建立严谨解析对应物的泛函导数。

3. 水平集演化与重初始化：
   水平集通过基于灵敏度的 Hilbertian 延展正则化（Hilbertian extension-regularisation）产生的速度场进行传输方程演化。作者提出了一个改进的惩罚项，其中包含了速度幅值 $|n_F \cdot \beta|$，以防止非设计区域被过度稳定。重初始化为符号距离函数采用带有人工粘性或连续内部惩罚项的稳态方法。

4. 孤立体积检测：
   提出了一种基于图的新颖方法来检测孤立体积。不同于以往依赖热对流问题的方法（这类方法无法区分仅由单层单元分隔的体积），该方法构建了一个符合性的切割网格连通图（使用多面体切割以确保符合性），并应用广度优先搜索（BFS）着色算法。这可以识别不连通的区域，从而允许应用惩罚项来约束这些体积内的平均位移或压力，进而使系统矩阵正则化。

5. 计算框架：
   该实现支持通过 GridapDistributed 进行串行和内存分布式计算。它利用 GridapEmbedded 在子三角剖分上进行计算，并利用 GridapTopOpt 进行拓扑优化循环，该循环结合了 Cea 方法（伴随方法）与 ASD 进行 PDE 受约束优化。

主要结果
该框架通过两个三维示例进行了验证：

1. 线性弹性轮的最小柔度： 一个涉及具有体积约束的三维弹性轮的拓扑优化问题。该问题在包含约 132 万个单元的非结构化背景网格上，利用 48 个处理器在分布式内存框架下求解。结果表明，该方法能够成功收敛，并能在无需解析推导形状导数的情况下处理复杂几何形状。
2. 流固耦合 (FSI)： 一个涉及与 Stokes 流相互作用的线性弹性结构的拓扑优化问题。该问题需要处理流体相与固体相之间的单向耦合。作者在目标函数中引入了一个惩罚项，以防止形成孤立的固体体积，否则这些体积会通过阻挡流体流动来人为地提高柔度目标。该问题通过 Stokes 方程和弹性方程的交替（顺序）求解来解决。

在这两个案例中，自动微分技术都成功计算了所需的梯度，且基于图的孤立体积检测有效地管理了不连通区域。

意义与主张
本文声称为非结构化网格上的水平集拓扑优化提供了一个通用且鲁棒的实现，避免了解析形状导数推导的负担。具体而言：

• 精度： 作者断言，其自动微分方法能以机器精度恢复解析导数，验证了 ASD 与 Berggren [16] 的“先离散后微分”方法之间的理论联系。
• 通用性： 通过将解析框架扩展到包括设计域与背景边界的交截情况，该方法填补了现有文献中的空白，这对于通用的拓扑优化问题至关重要。
• 灵活性： 使用 ASD 使得优化涉及复杂物理现象（例如 FSI、非标准函数空间如 $H(\text{div})$ 和 $H(\text{curl})$）的问题成为可能，而无需为每一个新泛函进行新的数学推导。
• 可扩展性： GridapTopOpt 的实现支持分布式计算，能够解决大规模非结构化三角剖分上的问题。

作者将这项工作定位为使非匹配拓扑优化对于更广泛的工程问题更加易于获取和应用的重要一步，通过自动化处理其中最具数学挑战性的环节。该软件已作为开源软件发布。