一种基于边界自适应分段方案的嵌入式等几何分析方法

《Computer-Aided Design》：An embedded isogeometric analysis method based on boundary adaptive segmentation scheme

　　

•开发了一种用于结果可视化的 CAD 模型预处理方案。
•开发了一种边界自适应分割方法以提高分析精度。
•研究了网格与边界不匹配对分析结果的影响。
•通过几个数值算例证明了该方法的有效性。

1. 引言
随着计算机技术的快速发展，计算机辅助工程 (CAE) 取得了显著进步，商业 CAE 软件目前广泛应用于航空航天、造船、建筑、汽车工程和电子产品设计领域。目前，市场上大多数主流 CAE 软件都是基于有限元分析 (FEA) 方法开发的 [1]。典型的有限元分析工作流程包括 CAD 模型导入、域分解、网格生成、材料属性和边界条件分配、数值求解以及分析结果的后处理。在这些步骤中，域分解和网格生成是最耗时的，因为它们需要创建独立于底层几何结构的分析模型。此外，网格离散化过程中引入的几何离散化误差以及网格单元之间的低阶连续性严重影响了有限元解的精度。因此，通常需要局部网格细化策略 [2] 来确保精度。这些问题从根本上源于有限元方法中几何模型与分析模型的独立性。

2000 年代初，Hughes、Cottrell 及其同事 [3,4] 提出了等几何分析 (IGA) 方法，为这些局限性提供了解决方案。IGA 方法采用 CAD 建模中常用的非均匀有理 B 样条 (NURBS) 基函数来构建分析模型，从而统一了几何模型和分析模型。在此框架下，基于 NURBS 的几何模型可以直接用于分析，只需简单的网格细化，无需传统网格划分。尽管经过近二十年的发展，但 IGA 尚未在工程应用中广泛采用。主要障碍在于两个方面：(i) NURBS 模型由控制点和 NURBS 基函数的张量积定义，这使得建模过程复杂、耗时且不直观；(ii) 工程结构通常几何高度复杂，在构建适合 IGA 的模型之前需要进行域分解。然后每个子域必须用 NURBS 表示，这实际上将预处理负担从分析阶段转移到了 CAD 建模阶段，而未解决根本问题。

除了 IGA 之外，嵌入域方法也具有高效的预处理工作流程。这些方法通过将复杂几何结构嵌入到简单的背景域中来避免传统 FEA 所需的域分解和网格生成。代表性方法包括扩展有限元法 (XFEM) [[5], [6], [7]]、切分有限元法 (CutFEM) [8]、聚合有限元法 (AgFEM) [9,10] 和有限胞元法 (FCM) [[11], [12], [13]]。文献 [14] 回顾并总结了上述各类嵌入域方法，并指出了嵌入域方法目前面临的三个主要挑战：

1. 狄利克雷边界条件的施加。在嵌入域方法中，背景网格通常不符合嵌入几何模型的边界。因此，无法像在传统有限元方法中那样直接在网格节点或控制点上施加狄利克雷边界条件。相反，需要采用“弱”施加技术，通过额外的边界积分项将狄利克雷边界条件纳入边界值问题的等效变分形式中。常用的“弱”施加方法包括罚函数法 [[15], [16], [17], [18]]、Nitsche 方法 [[19], [20], [21]] 和拉格朗日乘子技术 [[22], [23], [24], [25]]。罚函数法和 Nitsche 方法都需要仔细选择罚参数：值太小无法有效实施边界条件，而值太大则会加剧条件数问题。虽然拉格朗日乘子技术不涉及罚参数，但它们会导致鞍点问题，这要求适当选择离散拉格朗日乘子场以确保系统稳定。在本工作中，采用 Nitsche 方法来施加狄利克雷边界条件。

2. 小切分单元问题。在传统有限元分析中，网格单元的大小在网格生成过程中受到严格控制，以防止系统矩阵的条件数过大，从而导致平衡方程难以求解。然而，在嵌入域方法中，无法控制修剪单元内材料部分的形状和大小。当修剪单元的材料体积分数变得非常小时，所得系统矩阵的条件数往往非常大，导致数值不稳定。为解决此问题，嵌入域方法通常采用预处理策略来降低系统矩阵的条件数 [[26], [27], [28], [29]]。由于该问题不是本研究的主要关注点，我们转而采用一种受基于密度的拓扑优化启发的策略来处理数值奇异性 [30,31]。具体而言，为每个修剪单元的孔洞区域分配一个比固体材料小几个数量级的虚拟杨氏模量。这种虚拟刚度对最终分析结果的影响很小。

3. 修剪单元的积分。由于修剪单元内材料区域的几何构型多样，传统有限元方法中使用的高斯积分不能直接应用。修剪单元上的积分可分为两个方面：(i) 修剪单元材料部分的体积积分，以及 (ii) 修剪单元内包含的狄利克雷和诺伊曼边界上的面积分。目前，针对修剪单元的体积积分方法已进行了广泛的研究，包括八叉树细分方法 [[32], [33], [34], [35]]、切分单元重参数化方法 [[36], [37], [38]] 和多式积分方法 [[39], [40], [41]]。在本研究中，我们采用 Manuel Me?mer 等人于 2024 年开发的多式积分方法 [42]。该方法的核心思想是通过求解矩拟合方程 [43] 来构建每个修剪单元的积分规则（即积分点的位置和权重）。具体而言，首先采用交集算法获得包围每个修剪单元内材料区域的封闭曲面。基于该封闭曲面，利用散度定理 [44] 计算矩拟合方程所需的矩。然后将这些矩代入矩拟合方程，并应用点消元算法 [[45], [46], [47]] 来确定最终的积分点及其相应权重。由于通过交集算法对修剪单元进行了预处理，该方法在矩评估方面具有高精度。本工作主要关注修剪单元上的面积分。在大多数现有的基于 STL 的嵌入域方法中，狄利克雷和诺伊曼边界的积分通常直接使用 STL 模型的三角面片作为积分单位并应用高斯积分。然而，STL 三角面片与背景网格之间通常会发生不匹配，例如单个三角面片可能跨越多个背景单元。这种不匹配会不利影响边界积分的精度，如 4.1 节所示。文献 [42] 通过采用 CAD 算法对跨越多个背景单元的三角面片进行细分，然后在生成的子面片上进行边界积分，从而从根本上解决了不匹配问题。受此方法启发，我们提出了一种专门用于边界积分的边界自适应分割方案。该方案不仅消除了 STL 模型与背景网格之间的不匹配，还最小化了参与边界积分的三角面片数量。此外，所提出的方案可直接用于模拟结果的后处理和可视化。

总之，本文提出了一种基于边界自适应分割方案的嵌入等几何分析方法。该方法可以直接基于 STL 模型进行模拟，STL 模型是现有 CAD 系统中的常见数据格式。因此，它可以显著降低 CAD 和 CAE 系统之间数据交换的难度，对于实现 CAD/CAE 集成具有重要意义。在该方法中，使用 NURBS 背景网格来构建结构分析模型，保留了 IGA 固有的相邻单元间高阶连续性优势。狄利克雷边界条件使用 Nitsche 方法施加，而背景网格内切分单元的积分规则（即积分点坐标和权重）通过矩拟合策略和点消元算法获得。最后，应用边界自适应分割方案将几何边界表面离散化为完全符合背景网格的三角面片，从而减少边界积分误差，并实现基于离散几何边界的模拟结果可视化。

本文其余部分组织如下：第 2 节介绍所提出的嵌入等几何分析方法的基本原理。第 3 节提出了边界自适应分割方案。第 4 节通过几个三维算例展示了该技术的性能。最后，第 5 节得出结论。

2. 嵌入等几何分析
2.1. 基于 NURBS 的背景网格
在嵌入等几何分析中，背景域使用 NURBS 基函数构建。对于三维结构问题，NURBS 背景域 ?? 可定义为公式 (1)。(1)???(??,??,??)=∑????=1∑????=1∑????=1???????????????????(??,??,??)??????????????????????????????(??,??,??)=????,???(??)?????,???(??)?????,???(??)???????????∑?????=1∑?????=1∑?????=1?????,???(??)??????,???(??)??????,???(??)??????????????其中 ???????????????????(??,??,??) 表示 NURBS 基函数，由三个参数方向上的 B 样条基函数 ????,???(??),????,???(??),????,???(??) 及其相关权重系数 ?????????? 的乘积构成。由于用于构建背景网格的 NURBS 基函数对应的权重均等于 1，且参数空间中任意点处 B 样条基函数之和等于 1，上述表达式可等效重写为：(2)???(??,??,??)=??∑??=1??∑??=1??∑??=1????,???(??)?????,???(??)?????,???(??)???????????此处，??,??,?? 表示参数空间三个方向上的参数坐标；?????????? 表示第 ijl 个控制点坐标，对于三维问题 ??????????=[??????????,??????????,??????????]T，且 ??????????,??????????,?????????? 为三个空间方向上的物理坐标。整数 ??,??,?? 分别表示三个参数方向上单变量 NURBS 基函数的数量；????,???(??),????,???(??),????,???(??) 表示三个参数方向上的 B 样条基函数，且 ??,??,?? 为其相应阶数。为了简洁，可对三个参数方向上 B 样条基函数的乘积重新编号，允许将公式 (2) 重写为简化形式。(3)???(??,??,??)=??∑??=1?????(??,??,??)?????由于背景域不需要符合物理结构的几何形状，可以方便地将其选为矩形盒，如图 1 所示。背景域的长度、宽度和高度可根据嵌入结构的尺寸进行调整。然后使用 NURBS 特有的标准 k 细化策略对背景域进行细化，即首先升阶 NURBS 基函数的多项式阶数，然后插入新的节点。此过程生成 NURBS 背景网格。网格密度可根据结构的几何复杂度进行手动调整。值得注意的是，在本工作中，X、Y 和 Z 方向上的 NURBS 基函数阶数均设为 2。下载：下载高分辨率图像 (299KB) 下载：下载全尺寸图像图 1. 通过 k 细化生成的 NURBS 背景网格（红点为控制点）。

2.2. 边界值问题
图 2 展示了嵌入域分析的概念。图 2(a) 显示了结构静力分析的物理模型，据此可制定相应的边界值问题。(3a)Ω:????+??????=0Γ??:??=???????Γ??:??·??=???下载：下载高分辨率图像 (200KB) 下载：下载全尺寸图像图 2. 嵌入域分析：(a) 物理模型；(b) 背景域；(c) 嵌入物理模型的背景网格。根据边界值问题的数学描述，在结构域 Ω 内，应力场 ???? 的散度加上外加体力 ?????? 为零。在狄利克雷边界 Γ?? 上，满足狄利克雷条件，即边界上的位移等于给定位移 ???????。在诺伊曼边界 Γ?? 上，牵引矢量（定义为应力张量 ?? 在外法线方向上的投影）等于施加的边界载荷 ???。为求解此边界值问题，可采用变分法，得到以下等效积分形式。(4)∫Ω????·??dΩ=∫Ω??·??????dΩ+∫Γ????·???dΓ此处，?? 表示权函数，其在狄利克雷边界上为零。值得注意的是，公式 (4) 中的 ?? 代表公式 (3) 中出现的应力张量的向量形式。在传统有限元分析中，可以直接在有限元网格的节点上强制执行狄利克雷和诺伊曼边界条件。

然而，在嵌入域框架中，结构几何体完全嵌入在背景域 $\tilde{\Omega}$ 内部，如图2(a)所示。由于结构边界通常与背景网格节点不对齐，边界条件无法直接在节点层次上施加。为解决这一问题，我们采用Nitsche方法施加Dirichlet边界条件，该方法将边界条件转化为边界积分项。具体而言，通过对约束边界上三角面片进行场函数积分求值，将Dirichlet边界条件以边界积分项的形式纳入等价的变分格式中，并装配至整体刚度矩阵，从而在嵌入框架内实现边界条件的施加。相比之下，Neumann边界条件通过对边界三角面片上的外部载荷函数进行积分，并将其等价地分配给包含这些面片的背景网格单元所关联的控制点，从而形成式(7)中的外部载荷向量 $\mathbf{F}$。

施加Nitsche方法后边界值问题的等价变分格式如下：

(5) $\int_{\Omega} \boldsymbol{\omega} \cdot \boldsymbol{\sigma} \, d\Omega + \int_{\Gamma_D} \delta(\boldsymbol{\sigma}\mathbf{n})^T (\bar{\mathbf{u}} - \mathbf{u}) \, d\Gamma - \int_{\Gamma_D} \boldsymbol{\omega} \cdot (\boldsymbol{\sigma}\mathbf{n}) \, d\Gamma + \alpha \int_{\Gamma_D} \boldsymbol{\omega} \cdot (\bar{\mathbf{u}} - \mathbf{u}) \, d\Gamma = \int_{\Omega} \boldsymbol{\omega} \cdot \bar{\mathbf{b}} \, d\Omega + \int_{\Gamma_N} \boldsymbol{\omega} \cdot \bar{\mathbf{t}} \, d\Gamma$

(6) $\mathbf{n} = \begin{bmatrix} n_x & 0 & 0 & 0 \\ 0 & n_y & 0 & 0 \\ 0 & 0 & n_z & n_y \\ n_z & n_y & n_x & 0 \\ 0 & 0 & n_z & n_x \\ 0 & 0 & n_y & n_x \end{bmatrix}$

其中 $\mathbf{n}$ 表示Dirichlet边界上的外法向单位向量，在式(5)中，$\mathbf{n}$ 以矩阵形式表示，如式(6)所示。$\alpha$ 是确保Nitsche方法稳定性的稳定化参数。其取值必须恰当：若 $\alpha$ 过小，边界条件无法充分施加；若其过大，则可能导致整体刚度矩阵病态并引起数值不稳定。本研究中，稳定化参数根据工程经验选取。在第4节所示的所有数值算例中，将其设为 $10^{10}$，该值足以施加边界条件而不会在仿真过程中引入数值不稳定。

与式(4)相比，经Nitsche方法转换后的格式引入了三个额外的边界积分项，这使得边界积分的精度对嵌入域分析的整体解精度更为关键。值得注意的是，Nitsche方法属于边界条件"弱施加"技术的一类。除Nitsche方法外，其他常用方法还包括罚函数法和Lagrange乘子法。关于这些不同边界条件施加技术的详细比较，可参见文献[48]。

位移场由NURBS基函数与控制点位移系数的张量积表示。将该表示代入等价变分格式，式(5)即转化为如下平衡方程：

(7) $\left(\int_{\Omega} \mathbf{B}^T \mathbf{D} \mathbf{B} \, d\Omega - \int_{\Gamma_D} \mathbf{N}^T \mathbf{n} \mathbf{D} \mathbf{B} \, d\Gamma\right) \mathbf{U} + \left(\int_{\Gamma_D} \mathbf{B}^T \mathbf{D} \mathbf{n}^T \mathbf{N} \, d\Gamma + \alpha \int_{\Gamma_D} \mathbf{N}^T \mathbf{N} \, d\Gamma\right) (\bar{\mathbf{U}} - \mathbf{U}) = \mathbf{F}$

其中 $\mathbf{B}$ 为应变-位移矩阵，其元素为NURBS基函数对空间坐标的导数；$\mathbf{N}$ 为形状函数矩阵，由定义在背景网格上的NURBS基函数组成；$\mathbf{U}$ 为控制点处未知位移向量；$\bar{\mathbf{U}}$ 表示规定的位移约束；$\mathbf{F}$ 为控制点处通过数值积分获得的外部载荷向量；$\mathbf{D}$ 为材料弹性矩阵，对于三维问题可表示如下，其中 $E$ 和 $\mu$ 分别为杨氏模量和泊松比。

(8) $\mathbf{D} = \frac{E(1-\mu)}{(1+\mu)(1-2\mu)} \begin{bmatrix} 1 & \frac{\mu}{1-\mu} & \frac{\mu}{1-\mu} & 0 & 0 & 0 \\ \frac{\mu}{1-\mu} & 1 & \frac{\mu}{1-\mu} & 0 & 0 & 0 \\ \frac{\mu}{1-\mu} & \frac{\mu}{1-\mu} & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\mu}{2(1-\mu)} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\mu}{2(1-\mu)} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\mu}{2(1-\mu)} \end{bmatrix}$

**2.3 单元刚度矩阵的计算**

将结构几何体嵌入背景域后，背景网格单元可分为三类：(1) 完全包含在几何体内部的实体单元；(2) 完全位于几何体外部的空隙单元；以及(3) 与几何边界相交的裁剪单元，如图2(c)所示。为区分这三类单元，本研究采用基于STL几何模型的两步分类方案。第一步，按照文献[49]提出的策略识别所有裁剪单元，其对应的实现算法可在开源项目CGAL[50,51]中获取。第二步，通过判断质心相对于几何边界的位置，将所有非裁剪单元进一步划分为内部单元或外部单元[52]。

对于空隙单元，不存在材料，因此无需计算其刚度矩阵。对于填充有材料的实体单元，可使用标准Gauss求积进行刚度矩阵的求解。然而，裁剪单元仅部分填充有材料，部分为空隙，因此无法直接适用Gauss求积。为解决这一问题，本研究采用矩拟合方案结合点消除算法，为裁剪单元构建精确的积分规则。

对于定义在三维参数空间 $\hat{\Omega}$ 中的 $p$ 阶多项式函数 $f(\boldsymbol{\xi})$，其在参数域上的积分可通过有限个求积点处的函数值乘以相应权重之和来近似，即：

(9) $\int_{\hat{\Omega}} f(\boldsymbol{\xi}) \, d\hat{\Omega} \approx \sum_{j=1}^{m} f(\boldsymbol{\xi}_j) w_j$

其中 $\boldsymbol{\xi}_j = [\xi_j, \eta_j, \zeta_j]^T$ 表示第 $j$ 个求积点的参数坐标，$w_j$ 为对应的权重，$m$ 为求积点总数。在嵌入等几何分析中，参数空间中的NURBS函数表示为NURBS基函数与控制点的张量积。由于NURBS基函数为分段多项式函数，将 $S(\boldsymbol{\xi})$ 代入上式可得：

(10) $\int_{\hat{\Omega}} S(\boldsymbol{\xi}) \, d\hat{\Omega} = \sum_{i=1}^{n} P_i \int_{\hat{\Omega}} N_i(\boldsymbol{\xi}) \, d\hat{\Omega} \approx \sum_{i=1}^{n} P_i \sum_{j=1}^{m} N_i(\boldsymbol{\xi}_j) w_j$

将式(10)右端改写为矩阵形式，并令其等于NURBS基函数在裁剪单元实体部分上的精确积分，即构建矩拟合方程：

(11) $\begin{bmatrix} N_1(\boldsymbol{\xi}_1) & \cdots & N_1(\boldsymbol{\xi}_m) \\ \vdots & \ddots & \vdots \\ N_n(\boldsymbol{\xi}_1) & \cdots & N_n(\boldsymbol{\xi}_m) \end{bmatrix} \begin{bmatrix} w_1 \\ \vdots \\ w_m \end{bmatrix} = \begin{bmatrix} \int_{\hat{\Omega}_t} N_1(\boldsymbol{\xi}) \, d\hat{\Omega} \\ \vdots \\ \int_{\hat{\Omega}_t} N_n(\boldsymbol{\xi}) \, d\hat{\Omega} \end{bmatrix}$

其中 $\hat{\Omega}_t$ 表示裁剪单元实体部分对应的参数域。右端积分作为参考值，可利用散度定理将体积分转化为面积分来求解，详细步骤见文献[19]。本质上，矩拟合方法优化预定义的求积点集，使得计算结果逼近参考值。求积点的初始化遵循文献[53]。随后，将矩拟合方程重构为关于求积点及其权重的最小二乘问题：

(12) 求：$(\boldsymbol{\xi}, \mathbf{w}) = ([\boldsymbol{\xi}_1, \ldots, \boldsymbol{\xi}_m]^T, [w_1, \ldots, w_m]^T)$

最小化：$\mathbf{r}^2 = \mathbf{f}_c - \mathbf{f}(\boldsymbol{\xi})\mathbf{w} \leq \delta$

其中 $\|\cdot\|_2$ 表示欧几里得范数，$\mathbf{r}$ 为表示参考解与近似解之间差异的残差向量，$\mathbf{f}_c$ 为参考项（式(11)右端）的简化记号，$\mathbf{f}(\boldsymbol{\xi})\mathbf{w}$ 表示近似项（式(11)左端）。容差 $\delta$ 手动规定，本研究设为 $10^{-8}$。

基于该最小二乘问题，采用点消除算法迭代优化初始求积点集。在每次迭代步骤中，利用非负最小二乘法更新求积权重，并移除权重为零或极小的求积点。迭代持续进行，直至残差向量的欧几里得范数超过规定容差 $\delta$，最终获得求积点和权重。点消除算法的详细信息见文献[54]。需要注意的是，为避免裁剪单元实体占比极小时导致整体刚度矩阵病态，在装配刚度矩阵时，为裁剪单元的空隙区域赋予一个虚构的杨氏模量。该虚构模量设为实体材料杨氏模量的 $10^{-6}$ 倍，其对整体分析结果的影响可忽略不计。

**3. 边界自适应剖分方案**

如2.2节所述，嵌入域分析通过Nitsche方法在等价变分格式中引入额外的边界积分来施加Dirichlet边界条件。因此，边界积分的精度直接影响嵌入域分析的整体精度。为提高边界积分精度，常规做法是在导出三维几何时对STL模型的边界三角面片进行细化，即手动控制曲面三角化的网格尺寸。更细的三角化会增加边界表面上的三角单元数量，如图3所示。然而，该方法无法保证边界面片与背景网格之间的完全一致性。实际上，三角面片经常跨越两个或更多背景单元，导致积分误差。此外，过度细化将显著增加STL模型的存储需求。

> 下载：下载高分辨率图片（329KB）
> 下载：下载全尺寸图片
> 图3. 以STL模型形式导出的几何模型。

为克服上述局限性，本研究基于背景网格为STL模型开发了一种边界自适应剖分策略。该策略增强了嵌入域分析处理STL模型的灵活性，并实现了对几何边界的更精确处理。其优势在第4.1节中具体加以说明。该策略包含三个主要步骤，以图3中的几何模型为例进行说明，其中L形表面代表积分边界。

**步骤1.** 提取对应于积分边界的STL曲面。对于每个STL三角面片，根据其三个顶点计算对应的空间平面函数。然后，通过比较这些平面函数的参数来确定共面性。随后，根据端点坐标识别并标记面片的边：出现两次的边被归类为内部边，仅出现一次的边构成外轮廓。最后，移除外轮廓包围的内部面片边，并将三角面片合并为一个平面多边形，如图4所示。

> 下载：下载高分辨率图片（186KB）
> 下载：下载全尺寸图片
> 图4. 平面多边形的构建。

**步骤2.** 计算多边形轮廓边与背景网格面的交点，以及背景网格边与多边形的交点，如图5所示。

> 下载：下载高分分辨率图片（472KB）
> 下载：下载全尺寸图片
> 图5. 两类交点的计算。

**步骤3.** 利用步骤2中的交点及多边形顶点，对每个背景单元内的多边形区域进行Delaunay三角剖分。计算每个生成三角面片的质心，并丢弃质心位于多边形之外的面片。最终获得三角剖分结果，如图6所示。

> 下载：下载高分辨率图片（227KB）
> 下载：下载全尺寸图片
> 图6. 新三角面片的生成。

所提出的边界自适应剖分策略确保积分边界与背景网格完全一致，每个三角面片完全包含在单个背景单元内，并可最小化参与边界积分的三角面片数量。此外，该三角剖分可自动适应背景网格尺寸的变化以及积分边界与背景网格之间空间关系的变化，无需人工干预。图7(a1)-(c1)展示了三种不同情况：图7(a1)对应 $3\times3\times10$ 的背景网格；图7(b1)对应 $5\times5\times10$ 的背景网格；在图7(c1)中，L形边界嵌入在 $5\times5\times10$ 的网格中，沿 $Z$ 轴负方向平移70 mm，然后绕其底边逆时针旋转30°。图7(a2)-(c2)展示了相应的自适应剖分结果。结果表明，无论网格密度或相对位置如何，该策略始终能生成与背景网格完全一致的三角面片。

> 下载：下载高分辨率图片（781KB）
> 下载：下载全尺寸图片
> 图7. 不同背景网格及相对位置的自适应剖分结果。

除平面积分边界外，所提出的策略还可以处理曲面，如图8所示。值得强调的是，与后续的嵌入域分析相比，利用自适应剖分策略生成新积分边界的计算代价微乎其微。

> 下载：下载高分辨率图片（833KB）
> 下载：下载全尺寸图片
> 图8. 曲面的自适应剖分结果。

边界自适应剖分策略不仅可作为嵌入等几何分析中评估边界积分的预处理步骤，还支持后处理阶段。基于边界自适应剖分策略处理的几何边界，嵌入等几何分析所得仿真结果以位移云图、应力云图及其他物理场的形式进行可视化展示。

**4. 数值算例**

为验证所提出的具有边界自适应剖分方案的嵌入等几何分析方法的有效性，本节考察了四个三维结构算例。

所有数值算例均在 MATLAB R2022b 中进行，运行在一台配备 64 位 3.7 GHz CPU 和 64 GB 内存的台式计算机上。第一个算例分析了一个圆柱悬臂梁，以评估所提方法的准确性并检验自适应边界分割方案对数值结果的影响。第二个算例考虑了一个四分之一环形结构，将所提方法的性能与传统等几何分析 (IGA) 方法进行了比较。第三个算例分析了一个飞机发动机支架结构，并将使用所提方法获得的模拟结果与商业有限元分析软件的结果进行了比较。最后，第四个算例分析了一个转向节，以展示所提方法在处理复杂工程结构模拟问题的能力。

4.1. 圆柱悬臂梁
如图 9 所示，第一个基准问题是圆柱悬臂梁。圆柱长度为 8 米，端面半径为 1 米。一端完全固定，另一端施加垂直向下的均布载荷 F，大小为 0.01N/m2。圆柱材料的杨氏模量为 10 Pa，泊松比为 0。图 10 显示了导出为 STL 模型后的相应几何形状，其中两个端面均被视为积分边界。根据参考文献 [55] 中的 Timoshenko 梁理论，悬臂梁中心线上的挠度解析解 ???? 可以表示为：
(13) ????=FA?[??2?(3?L???)6?EI???GA??]
其中 F 是均布载荷的大小，A 是圆柱的横截面积，?? 表示中心线上某点到固定端面的距离，I 是截面二次矩，I=?????4/4，G 是剪切模量，其与杨氏模量和泊松比的关系为 G=E/(2+2???)，?? 是剪切系数，给定由：??=(6+12???+6???2)/(7+12???+4???2)。
下载：下载高分辨率图像 (92KB) 下载：下载全尺寸图像 图 9. 圆柱悬臂梁。
下载：下载高分辨率图像 (309KB) 下载：下载全尺寸图像 图 10. STL 模型及其积分边界。

采用嵌入域方法对圆柱悬臂梁进行模拟分析。首先，将图 10 所示的 STL 模型嵌入到使用 NURBS 基函数构造的 3 × 3 × 5 背景网格中。图 11(a1) 和 (b1) 分别展示了使用本文提出的边界自适应分割策略和参考文献 [42] 中采用的边界细分策略处理后的积分边界，用于计算矩拟合方程中的矩。两种分割策略的比较表明，与参考文献 [42] 相比，所提方法将三角面片的数量减少了 60% 以上，并且随着构成积分边界的共面且连接的三角面片数量增加，这一优势变得越来越明显。此外，使用参考文献 [42] 中策略获得的分割结果取决于 STL 模型的初始三角面片分布。如图 11(b2) 所示，当圆柱模型的 STL 导出参数更改时，积分边界上的三角面片分布随之改变，导致使用参考文献 [42] 中的方法获得不同的分割结果。相比之下，如图 11(a2) 所示，使用所提策略获得的分割结果与图 11(a1) 中的结果几乎相同，表明所提方法具有优越的鲁棒性。
下载：下载高分辨率图像 (772KB) 下载：下载全尺寸图像 图 11. 不同方案下的离散几何边界：(a) 边界自适应分割方案；(b) 参考文献 [42] 中的边界细化方案。

图 12(a) 和 (b) 分别展示了与图 11(a1) 和 (b1) 中所示的两种分割策略相对应的嵌入域分析结果。在图 12 中，左列显示总位移云图，右列显示冯·米塞斯应力云图。可以观察到，尽管两种分割策略之间的三角面片分布不同，但最终的分析结果相同。这表明所提边界自适应分割策略可以在保持分析精度的同时最小化三角面片的数量，从而降低求解过程中与边界积分相关的计算成本，如表 1 所示。
下载：下载高分辨率图像 (465KB) 下载：下载全尺寸图像 图 12. 圆柱悬臂梁的模拟分析结果：(a) 边界自适应分割方案；(b) 参考文献 [42] 中的边界细化方案。

表 1. 求解消耗。
处理的边界 | 边界积分时间 (s) | 组装时间 (s) | 求解时间 (s)
--- | --- | --- | ---
图 11(a1) | 0.266 | 0.375 | 0.008
图 11(b1) | 0.757 | 0.377 | 0.008

接下来，采用如图 13(a)-(e) 所示的五种不同单元数量的背景网格进行模拟，以研究沿圆柱轴线的网格细化对所提方法准确性的影响。
下载：下载高分分辨率图像 (671KB) 下载：下载全尺寸图像 图 13. 沿圆柱轴线具有不同单元数量的五种背景网格：(a) 网格 3 × 3 × 3；(b) 网格 3 × 3 × 4；(c) 网格 3 × 3 × 5；(d) 网格 3 × 3 × 6；(e) 网格 3 × 3 × 7；
图 14 绘制了嵌入等几何分析获得的数值挠度 ???? 与圆柱中心线上解析解 ???? 之间的相对误差。横轴表示每个轴点到固定端面的距离。图 13(a)-(e) 中五种背景网格对应的最大相对误差分别为 1.1%、0.46%、0.18%、0.068% 和 0.030%。显然，随着背景网格沿轴线的细化，嵌入等几何分析的准确性提高并收敛于解析解。
下载：下载高分辨率图像 (280KB) 下载：下载全尺寸图像 图 14. 五种背景网格的挠度相对误差。

接下来，保持梁轴线上的单元数量不变，改变圆柱横截面上的单元数量，以检验边界自适应分割方案的效果。如图 15(a)-(c) 所示，采用了三种不同横截面分辨率的背景网格，在 x、y、z 方向上的单元数量分别为 3 × 3 × 10、5 × 5 × 10 和 8 × 8 × 10。首先在不应用边界自适应分割方案的情况下进行模拟。所有模拟结果均基于使用 8 × 8 × 10 背景网格通过边界自适应分割策略处理的几何边界进行可视化，如图 15(d) 所示。
下载：下载高分辨率图像 (644KB) 下载：下载全尺寸图像 图 15. 具有不同横截面分辨率的三种背景网格：(a) 网格 3 × 3 × 10；(b) 网格 5 × 5 × 10；(c) 网格 8 × 8 × 10；(d) 用于可视化的离散几何边界。

图 16(a)-(c) 分别展示了未应用边界自适应分割方案时，对应图 15(a)-(c) 的三种嵌入域背景网格的圆柱悬臂梁模拟分析结果。对于每种情况，左侧显示总位移云图，右侧显示冯·米塞斯应力云图。如图 16(b) 和 (c) 所示，固定端附近的冯·米塞斯应力幅值和分布均出现显著误差。这是因为背景单元的尺寸小于集成边界上的 STL 三角面片的尺寸，导致背景网格与积分边界不匹配。因此，变分方程中边界积分的评估变得不准确，导致结果错误。
下载：下载高分辨率图像 (500KB) 下载：下载全尺寸图像 图 16. 未应用边界自适应分割方案时圆柱悬臂梁的模拟分析结果：(a) 网格 3 × 3 × 10；(b) 网格 5 × 5 × 10；(c) 网格 8 × 8 × 10。

随后，应用自适应分割方案。图 17(a)-(c) 显示了为三种不同背景网格生成的新三角面片分布。细分后的三角面片数量从左到右分别为 91、115 和 176。
下载：下载高分辨率图像 (391KB) 下载：下载全尺寸图像 图 17. 三种不同背景网格的自适应分割边界：(a) 网格 3 × 3 × 10；(b) 网格 5 × 5 × 10；(c) 网格 8 × 8 × 10。

图 18(a)-(c) 显示了应用自适应分割方案后悬臂梁的相应模拟结果。与图 16 的比较清楚地表明，带有边界自适应分割的嵌入等几何分析有效地消除了由背景网格和边界不匹配引起的数值误差。此外，该方案增强了嵌入等几何分析方法的鲁棒性，同时避免了每当背景网格改变时重新生成 STL 模型的额外预处理成本。
下载：下载高分辨率图像 (475KB) 下载：下载全尺寸图像 图 18. 应用边界自适应分割方案后的模拟分析结果：(a) 网格 3 × 3 × 10；(b) 网格 5 × 5 × 10；(c) 网格 8 × 8 × 10。

4.2. 四分之一环形
如图 19 所示，第二个数值算例是一个四分之一环形结构，其几何参数如下：内半径为 1 米，外半径为 2 米，厚度为 1 米。结构的一端完全固定，另一端承受垂直向上的均布载荷 F，大小为 0.01N/m2。结构材料的杨氏模量为 10 Pa，泊松比为 0.3。使用传统等几何分析 (IGA) 方法和所提出的带有边界自适应分割方案的嵌入等几何分析方法进行模拟。
下载：下载高分辨率图像 (102KB) 下载：下载全尺寸图像 图 19. 四分之一环形及其工况。

对于 IGA 方法，必须构建结构的 NURBS 模型，如图 20 所示。左侧面板显示了初始 NURBS 模型，其中 ??,??,?? 参数方向上的节点向量均为 [000111]。NURBS 基函数在所有三个参数方向上的阶数均设为 2，图中的红点表示控制点。然后使用 h-细化策略对初始模型进行细化，得到图 20 右侧所示的细分模型，总共包含 9 × 8 × 16=1152 个跨度单元。
下载：下载高分辨率图像 (379KB) 下载：下载全尺寸图像 图 20. 四分之一环形的 NURBS 模型。

图 21 说明了四分之一环形 STL 模型嵌入到背景网格中的情况。在这种情况下，使用了 15×15×8 的背景网格。总共激活了 1152 个单元，包括 528 个实体单元和 624 个修剪单元。经过边界自适应细分后，Dirichlet 和 Neumann 边界表面均由 128 个新的三角面片表示，用于可视化的离散几何模型如图 21 所示。基于图 20 和图 21 中的预处理步骤，我们随后进行了 IGA 和嵌入 IGA 分析。相应的模拟结果如图 22 所示，其中上方面板显示总位移云图，下方面板显示冯·米塞斯应力云图。
下载：下载高分辨率图像 (732KB) 下载：下载全尺寸图像 图 21. 四分之一环形 STL 模型嵌入到背景网格中。
下载：下载高分辨率图像 (577KB) 下载：下载全尺寸图像 图 22. 模拟结果：(a) 传统 IGA 方法；(b) 嵌入 IGA 方法。

图 22 中结果的比较表明，两种方法产生的位移和应力场基本相同，仅在局部区域有轻微差异。IGA 获得的最大位移和最大冯·米塞斯应力分别为 0.04821 m 和 0.1549 Pa，而嵌入 IGA 产生的结果为 0.04826 m 和 0.1566 Pa。两种方法之间最大位移和最大冯·米塞斯应力的相对误差分别为 0.103% 和 1.09%。应力相对误差比位移相对误差高一个数量级，这可以归因于应力取决于位移的梯度，因此对变化更敏感。

表 2 展示了该算例中 IGA 方法和嵌入 IGA 方法的求解消耗。如表 2 所示，两种方法的激活单元数量相同。

然而，嵌入式IGA方法涉及稍多的活跃控制点，导致更大的全局刚度矩阵，从而使得平衡方程的求解时间更长。尽管如此，两种方法在计算成本上最显著的差异在于预处理阶段。传统的等几何分析需要将图19所示的B-rep模型转换为NURBS表示。目前，这一过程严重依赖人工干预：首先必须将B-rep模型转化为节点向量、控制点坐标和权重系数，然后基于这些数据构建NURBS模型（如图20所示）。在本研究中，该转换过程耗时203秒，且随模型复杂度增加，时间成本进一步增加。相比之下，嵌入式等几何分析的预处理阶段仅包括将B-rep模型导出为STL格式以及进行边界积分。这两个步骤均可通过算法过程自动执行，所需计算时间可忽略不计。具体而言，在单元刚度矩阵的计算和组装方面，嵌入式框架中的实体单元相同，这意味着其刚度矩阵只需计算一次，从而降低了计算成本。相比之下，传统的等几何分析需要对每个单元分别进行刚度矩阵计算，导致更高的计算时间。总体而言，对于复杂工程结构的分析，嵌入式等几何分析方法相比传统IGA方法具有明显优势。

表2. 求解消耗。
| 方法 | 单元数 | 全局刚度矩阵 | 预处理时间（s） | 组装时间（s） | 求解时间（s） |
|------|--------|--------------|----------------|--------------|--------------|
| IGA | 1152 | 5940×5940 | 203 | 9.42 | 10.160 |
| 嵌入式IGA | 1152 | 6240×6240 | 0.494 | 5.625 | 0.198 |

4.3 航空发动机支架

如图23所示，第三个数值算例是一个喷气发动机支架结构模型（下载链接：https://www.3dwhere.com/models/m4b65km2yzl6vvcq）。支架材料赋予杨氏模量为10 Pa，泊松比为0.3。在支架顶部红色高亮的铰接孔上施加大小为0.01N/m2的垂直向上的均布载荷F，同时在黄色高亮的四个螺栓孔上施加完全固定位移约束。仿真使用所提出的带边界自适应划分方案的嵌入式等几何分析方法和常规有限元分析（FEA）进行。FEA使用商业软件Ansys 2022R2进行。

图23. 喷气发动机支架结构。

对于该算例，采用七种具有不同单元数的背景网格进行嵌入式等几何分析。沿x、y、z方向的单元数以及这七种网格的活跃单元总数、实体单元数和切割单元数列于表3中。随着背景网格的逐步细化，实体单元与切割单元的比值相应增加。图24展示了支架模型转换为STL表示后嵌入60×40×20背景网格的过程。STL模型包含1646个三角面片。与黄色螺栓孔相关的Dirichlet边界曲面和与红色铰接孔相关的Neumann边界曲面进行自适应细分，所得三角面片分布如图24中的绿色圆柱面所示。自适应细分后，两个铰接孔和四个螺栓孔分别由1960个和1092个新生成的三角面片表示。用于可视化的离散几何模型如图24底部所示。

表3. 七种背景网格。
| 背景网格 | 总活跃单元数 | 实体单元数 | 切割单元数 |
|----------|--------------|------------|------------|
| 42×28×14 | 6594 | 3336 | 3258 |
| 45×30×15 | 7807 | 4103 | 3704 |
| 48×32×16 | 9920 | 5520 | 4400 |
| 51×34×17 | 11440 | 6528 | 4912 |
| 54×36×18 | 13101 | 7672 | 5429 |
| 57×38×19 | 16036 | 9817 | 6219 |
| 60×40×20 | 18214 | 11358 | 6856 |

图24. 支架模型嵌入背景网格的过程。

使用60×40×20背景网格的嵌入式IGA和FEA得到的总位移云图分别如图25(a)和(b)所示，相应的冯·米塞斯应力云图分别如图26(a)和(b)所示。在FEA中，支架几何采用四节点四面体单元离散，共有202316个单元。对比图25(a)和(b)以及图26(a)和(b)可知，两种方法获得的总位移和冯·米塞斯应力分布高度吻合，仅在量值上存在微小差异。具体而言，最大总位移差异小于1.0%，最大冯·米塞斯应力差异小于5.0%。

(14) $er = \frac{|u_{EIGA}^{max} - u_{FEA}^{max}|}{u_{FEA}^{max}}$

图25. 支架的总位移云图：(a) 嵌入式等几何分析方法；(b) 有限元分析方法。

图26. 支架的冯·米塞斯应力云图：(a) 嵌入式等几何分析方法；(b) 有限元分析方法。

为了评估嵌入式等几何分析的收敛性，将商业有限元分析得到的最大总位移记为 $u_{FEA}^{max}$，作为参考解。使用公式(14)评估表3中七种背景网格下嵌入式等几何分析预测的最大总位移与参考解之间的相对误差，所得值分别为10.14%、7.49%、4.89%、3.48%、2.50%、1.20%和0.62%。图27展示了相对误差与背景网格中活跃单元数之间的关系。显然，随着背景网格的逐步细化，嵌入式等几何分析结果稳步收敛于参考解。

图27. 相对误差与活跃单元数的关系。

值得注意的是，由于FEA采用一阶线性单元，位移和应力场在单元边界上仅具有 $C_0$ 连续性。相比之下，本算例中嵌入式等几何分析使用二次NURBS基函数，在单元边界上实现了 $C_1$ 连续性，从而能够实现更平滑的过渡。

4.4 转向节

为了进一步展示所提出方法在处理复杂工程结构仿真分析方面的能力，最后一个算例考虑一个具有复杂几何和小几何特征的转向节（下载链接：https://grabcad.com/library/steering-knuckle-rh-1），如图28所示。这里，$\Gamma_N^1$、$\Gamma_N^2$ 和 $\Gamma_N^3$ 表示Neumann边界，在其上施加垂直于表面且指向内部的均布载荷，大小分别为2 MPa、1 MPa和0.3 MPa。$\Gamma_D$ 表示Dirichlet边界，在此处结构完全约束。转向节的材料特征为杨氏模量210000 MPa，泊松比0.3。与4.3节中的算例类似，仿真使用所提出的嵌入式等几何分析方法和常规有限元分析（FEA）进行。FEA使用商业软件Ansys 2022R2进行。

图28. 转向节。

图29展示了转向节转换为STL模型后嵌入背景网格的过程。STL模型由6174个三角面片组成。在本算例中，使用50×50×50的背景网格，其中总共激活11324个单元，包括3774个实体单元和7550个切割单元。与蓝色圆形平面相关的Dirichlet边界曲面以及与红色圆形平面相关的三个Neumann边界曲面进行自适应细分。所得三角面片分布如图29中的绿色圆形平面所示。自适应细分后，$\Gamma_N^1$、$\Gamma_N^2$、$\Gamma_N^3$ 和 $\Gamma_D$ 分别由155、232、327和311个三角面片表示。用于可视化的离散几何模型如图29底部所示。

图29. 转向节嵌入背景网格的过程。

使用50×50×50背景网格的嵌入式IGA和FEA得到的仿真结果分别如图30(a)和(b)所示，其中左半部分显示总位移云图，右半部分显示冯·米塞斯应力云图。在FEA中，转向节几何采用四节点四面体单元离散，共有307260个单元。对比嵌入式等几何分析和有限元分析的结果表明，总位移分布高度吻合，最大总位移差异小于1.2%。冯·米塞斯应力分布也基本一致，明显的差异主要集中在图30中红框标注的应力集中区域。在该区域内，最大冯·米塞斯应力差异为5.51%。

图30. 仿真结果：(a) 嵌入式等几何分析方法；(b) 有限元分析方法。

造成该差异的主要原因是应力集中区域的应力状态高度依赖于局部几何特征。然而，在嵌入式等几何分析方法中，几何边界主要由切割单元表示，这些单元在准确捕捉边界几何方面的能力有限。因此，在这些集中区域容易出现应力偏差。细化背景网格会增加应力集中区域中切割单元的数量，从而改善这些区域边界的几何表示精度。如图31所示，图31(a)至(c)分别展示了使用40×40×40、45×45×45和50×50×50三种背景网格的嵌入式等几何分析得到的冯·米塞斯应力分布。可以观察到，随着网格的细化，应力分布逐渐收敛于图30中显示的有限元分析结果。然而，网格细化也会导致计算成本的增加。

图31. 嵌入式等几何分析结果：(a) 网格40×40×40；(b) 网格45×45×45；(c) 网格50×50×50。

5. 结论与展望

本文提出了一种基于背景网格边界自适应划分方案的嵌入式等几何分析方法。该方法中采用的边界自适应划分方案能够自动生成与背景网格完全匹配的积分边界，并且能够基于背景网格对整个几何模型边界进行重新三角化，从而生成可用于可视化的模型。此外，它还允许几何边界上三角面片的分布适应模型与背景网格之间空间关系以及背景网格单元数量的变化。基于第4节中展示的基准研究，可得出以下四个结论：

(1) 网格与边界不匹配的影响。背景网格单元与STL积分边界的三角面片之间的不匹配会导致嵌入式域分析中的误差，特别是当背景网格单元尺寸远小于STL面片尺寸时。所提出的边界自适应划分方案完全消除了此类不匹配，从而提高了嵌入式等几何分析的边界积分精度。

(2) 边界积分的三角面片数量。与文献[42]中的边界细分策略相比，所提出的边界自适应划分策略不仅能够生成与背景网格完全匹配的积分边界，还能通过包含共面单元的合并过程，将参与边界积分的三角面片数量降至最低。当积分边界中连通且共面的三角面片数量增加时，这一优势日益显著。

(3) 与传统等几何分析的对比。不同于传统的等几何分析，所提出的方法不需要构建NURBS模型。相反，可以直接使用由商业CAD软件生成的STL模型进行仿真。这显著提高了预处理效率，并便于对复杂工程组件进行分析。

(4) 与有限元分析的对比。嵌入式等几何分析和有限元分析所得结果总体一致。然而，在应力集中区域可能会出现应力场差异。这主要是因为背景网格不符合几何边界，切割单元无法准确表示局部几何特征，从而导致应力偏差。虽然细化背景网格可以改善这些区域的精度，但也会增加计算成本。

在未来的研究中，我们计划研究一种结合截断层次B样条（THB样条）技术的嵌入式等几何分析方法，以实现背景网格的局部细化。该方法旨在以计算成本的小幅增加来提高局部几何边界附近的分析精度，从而解决嵌入式域方法的自适应分析难题。

CRediT作者贡献声明

赵凡：撰写—初稿、可视化、验证、方法论。高亮：撰写—审阅与编辑、监督。王泽民：撰写—审阅与编辑。兰新强：撰写—审阅与编辑、监督、资金获取。李浩：撰写—审阅与编辑、监督、方法论、资金获取。

赵凡 | 高亮 | 王泽民 | 兰新强 | 李浩