Laplacian regularized eikonal equation with Soner boundary condition on polyhedral meshes

技术摘要：基于多面体网格的拉普拉斯正则化 Eikonal 方程与 Soner 边界条件

问题陈述
本文研究了在三维计算域 $\Omega$ 内，计算给定物体（由曲线、曲面或边界部分 $\Gamma$ 表示）的距离函数的数值计算问题。该域使用多面体单元进行离散化，这是处理复杂边界形状的工业问题的常见要求。目标是求解 Eikonal 方程的粘性解：
$$|\nabla u(x)| = 1, \quad x \in \Omega \setminus \Gamma$$
且满足 $u(x) = 0$ 在 $\Gamma$ 上。一个关键挑战在于必须在定义域边界 $\partial \Omega \setminus \Gamma$ 上应用 Soner 边界条件（或状态约束条件）：
$$\nu(x) \cdot \nabla u(x) \geq 0$$
其中 $\nu$ 为外法向量。若不满足此条件，特别是在非凸定义域中，数值解可能无法收敛到正确的粘性解。

现有方法面临以下局限性：

• 快速行进法 (FMM)： 虽然高效，但其堆数据结构在利用领域分解进行并行计算时难以高效维护。
• 快速扫描法 (FSM)： 对简单定义域有效，但依赖于在通用多面体网格上并不容易实现的特定排序。
• 时间松弛形式 (Time-Relaxed Formulations)： 求解随时间变化的 Eikonal 方程（随时间传播前沿）虽然稳健，但计算成本高昂，尤其是在感兴趣区域远离 $\Gamma$ 或使用具有大量单元的细密网格时。

方法论
作者提出了一种基于单元中心有限体积法 (FVM) 的数值算法，用于求解拉普拉斯正则化 Eikonal 方程。受消失粘性法（vanishing viscosity method）的启发，通过引入拉普拉斯项对控制方程进行修改：
$$-\epsilon \Delta u_\epsilon(x) + |\nabla u_\epsilon(x)| = 1, \quad x \in \Omega \setminus \Gamma$$
其中 $u_\epsilon = 0$ 在 $\Gamma$ 上，并在 $\partial \Omega \setminus \Gamma$ 上应用 Soner 边界条件。这里 $\epsilon > 0$ 是一个正则化参数。

关键方法组成部分包括：

1. 线性化： 利用已知的近似值 $u_{\epsilon'}$（来自之前具有较大 $\epsilon$ 的迭代）将非线性项 $|\nabla u_\epsilon|$ 线性化，从而形成一个线性对流-扩散方程。对流速度定义为 $v = \nabla u_{\epsilon'} / |\nabla u_{\epsilon'}|_\sigma$。
2. 多面体网格上的离散化： 作者改编了适用于多面体单元的通量平衡近似法。通过将面（faces）细分为三角形来处理非平面面。该离散化分别处理内部面和边界面。
3. Soner 边界条件的实现： 该数值方案通过将边界三角形根据法向通量的符号进行分组，显式地强制执行 Soner 条件：将对应于边界流出（可能违反条件的项）的项设为零，而流入项则通过改进的基于流入的梯度进行处理。
4. 顺序正则化策略： 为确保收敛至粘性解，算法求解一系列具有递减正则化参数 $\epsilon_n$ 的问题。参数通过经验公式 $\epsilon_n = (h_L)^{1/2n}$ 选择，其中 $h_L$ 是特征网格长度。步骤 $n-1$ 的解作为步骤 $n$ 的初始猜测值。
5. 并行化： 该方法利用标准的单元中心 FVM 结构和 1-ring 面邻域 进行领域分解。这使得无需像堆（heap）那样复杂的全局数据结构即可实现简单的并行化实施。

主要贡献

• 针对多面体网格的算法： 一种稳健的单元中心 FVM 方案，能够在通用多面体网格上求解带有 Soner 边界条件的 Eikonal 方程，克服了对锐角三角形或四面体有几何限制的现有方法。
• 处理非线性与正则化： 提出了一种特定的线性化 Eikonal 方程及选择递减正则化参数的策略，以平衡数值稳定性（避免由于 $\epsilon$ 过小导致的振荡）与准确性（避免由于 $\epsilon$ 过大导致的过度平滑）。
• 计算效率： 证明了拉普拉斯正则化方法相比于时间松弛形式具有显著的计算成本降低，特别是在处理大规模网格或远离 $\Gamma$ 的区域时，这归功于椭圆正则化带来的无限传播速度。

数值结果
作者通过十个测试用例（EX1–EX10）验证了该方法，涉及球体、平面、圆盘以及具有不同 $\Gamma$ 配置的复杂定义域。

• 收敛性： 对于光滑解（例如距离球体的距离），该方法在 $L^1$ 范数下实现了二阶实验收敛阶 (EOC)。对于非光滑解（例如靠近角点或奇异点处），EOC 趋于 1，符合理论预期。
• 参数敏感性： 实验确认了选择 $\epsilon$ 至关重要。$\epsilon$ 过小会导致不稳定，而 $\epsilon$ 过大则会降低精度。研究发现，采用 $\epsilon_n = (h_L)^{1/2n}$ 且 $n=5$ 的序列可获得最优结果。
• 性能比较： 与在同一代码库中实现的随时间演化的双向 Eikonal 方程求解器相比，所提方法展示了巨大的加速效果。随着单元数量的增加以及感兴趣区域距离 $\Gamma$ 的增加，效率提升愈发显著。在某些情况下，在相同的误差容限下，所提方法的运行速度几乎快了 70 倍。

意义与主张
本文声称，所提出的算法为在复杂、工业级多面体网格上计算距离函数提供了一种实用且高效的解决方案。其主要意义在于：

1. 并行可扩展性： 能够使用具有 1-ring 邻域的标准领域分解来实现求解器，使其适用于内存消耗成为瓶颈的高性能计算环境。
2. 成本降低： 能够比时间推进法大幅降低计算成本，特别是在计算大型定义域或远离源物体时。
3. 稳健性： 在有限体积框架内成功应用 Soner 边界条件，以确保在非凸定义域上获得正确的粘性解。

作者得出结论，该方法已准备好集成到标准的单元中心 FVM 代码中，并且对于需要细微尺度 3D 模拟的应用（如燃烧建模、湍流建模（壁面距离）以及心脏电生理学）特别具有优势。