作者： Huu Phuoc Bui、Satyendra Tomar、Stéphane P.A. Bordas
发表期刊： Computer Methods in Applied Mechanics and Engineering
发表时间： 2018年10月接收，2018年11月在线发表，正式刊出为2019年第345卷
本研究属于计算力学与生物医学工程交叉领域，核心目标是解决实时手术模拟中网格生成与几何描述之间的耦合难题。在传统有限元方法（Finite Element Method， FEM）中，计算网格必须与模拟对象的几何边界保持一致，即所谓“共形网格”（conforming mesh）。然而，医学模拟所涉及的解剖结构（如肝脏、大脑、前列腺等）通常具有极其复杂的几何形态，且不同患者之间的个体化几何差异显著。当针对患者个体化几何进行模拟时，每次都需要重新生成高质量的共形网格，这一预处理过程的计算成本极高，严重制约了实时仿真的可行性。
为应对这一挑战，作者提出了共旋Cut有限元方法（Corotational Cut Finite Element Method， Corotational CutFEM）。该方法仅需一个不必然贴合模拟对象边界或界面的背景网格（background mesh），即可自动捕捉模拟对象的表面几何细节。研究旨在实现三个核心目标：（1）使离散化过程独立于几何描述；（2）避免复杂几何的网格生成困难；（3）在保持标准FEM精度的同时达到实时仿真所需的计算速度。
研究首先建立软组织内针插入问题的数学描述。组织域Ω通过边界条件分为Dirichlet边界γ_u和Neumann边界γ_t，控制方程包含位移场、Cauchy应力张量、体力项以及针与组织间的交互力。材料本构关系采用基于共旋公式的线弹性模型，以正确处理软组织在手术仿真中常见的大旋转变形问题。针的建模采用Euler–Bernoulli梁理论，使用一维单元离散，每个节点含有6个自由度（3个平移和3个转动）。通过引入有限维空间v_h，连续变分问题被转化为离散问题，最终得到质量矩阵、阻尼矩阵和刚度矩阵组装而成的系统方程。时间离散采用隐式向后Euler格式，以数值稳定性换取较大的时间步长。
针与组织的交互通过三种约束来建模：组织表面的穿刺约束（满足Kuhn–Tucker条件）、针尖处切割约束以及沿针杆的摩擦约束（遵循Coulomb摩擦定律）。所有交互力均通过Lagrange乘子法求解。
对于浸入边界或非共形界面的几何离散，研究提出了系统的处理方案。首先，利用水平集方法（level-set method）识别哪些背景网格单元被表面切割：若一条边的两个端点相对表面的有符号距离函数符号相反，则该边被界面切割，交点的重心坐标由距离函数值确定。切割单元与表面的交面可以是三角形或四边形两种情况。
由于背景网格可能较粗而界面表面为曲面，可能出现“无效切割单元”情况——例如一条边被切割两次，或四面体被切割的边数不在3至4之间。为此，作者采用递归细分策略：将无效切割单元递归嵌入一组子单元，直至所有子单元均成为有效切割单元。此过程不引入新的自由度，仅用于辅助积分。
此外，切割单元可能被表面切出极小体积的部分，这会严重影响系统矩阵的条件数和稳定性。为此，作者提出移动背景节点的稳定化策略：当切割后较小部分的体积小于单元总体积的5%时，将该单元最靠近交点的节点沿表面外法线方向移动一个与单元尺寸成比例的随机距离，从而消除极小切割体积。
实现层面，作者将包围盒细分为若干子立方体，建立每个表面三角形与相邻子立方体、以及每个子立方体与相邻四面体的双向索引关系。通过这种空间哈希结构，可以高效地从表面三角形出发找到潜在切割的四面体，以及从域边界出发向外传播以标记外部单元。切割单元被嵌入由8个模板四面体组成的集合，模板节点位于单元边的中点，若边被实际切割则中点移至实际交点处。最终使用树状数据结构管理背景单元与其嵌入子单元之间的层级关系。
数值积分方面，完全位于域内或域外的单元按经典FEM方式积分；切割单元的积分被分割为内、外两个子域，分别累加各子四面体的贡献。对于线性四面体单元，由于应变-位移矩阵在单元内为常数，每个子四面体的刚度贡献可简化为B矩阵转置、材料张量和子体积的乘积。
共旋公式通过极分解提取单元总位移中的刚体运动部分，从旋转后的当前构型计算单元变形，有效地消除了大旋转变形下线性弹性模型产生的伪应力。对于被切割的单元，由于自由度仍仅定义在父单元节点上，旋转矩阵的计算方式不变。
Neumann边界条件的处理采用主-从（master–slave）方案：作用在隐式表面上的外力通过重心坐标映射到包含该作用点的背景网格单元节点上。Dirichlet边界条件则通过Lagrange乘子法隐式施加，形成带约束的系统方程。整个带约束系统采用三步求解策略：先对无约束系统进行Cholesky分解并求解自由解，然后利用柔度矩阵求解Lagrange乘子，最后求解修正解并更新总解。该策略避免了增广矩阵的不定性问题，且适合GPU并行化。
研究首先验证了几何逼近的收敛性。以直径为1.4的球面为测试对象，将其隐式嵌入背景网格，计算CutFEM离散化后的体积与离散表面精确体积（1.40005）之间的误差。结果表明体积误差随网格细化以最优速率收敛，其收敛阶与位移误差的L2范数收敛阶近似相同。
力学求解的收敛性通过拉伸试验和弯曲试验进行。试验对象为6×2×2 mm的梁，内部嵌入一个直径1.4的球面（材料属性在球内外一致）。采用四组不同密度的四面体网格（节点数从7×3×3到49×17×17），以97×33×33节点细网格的经典FEM解作为参考。结果显示，L2范数的位移误差收敛阶为O(n^{-2⁄3})，能量范数收敛阶为O(n^{-1⁄3})，均与理论最优势一致。
在梁几何模型上，CutFEM与经典FEM在不同网格密度下的位移结果几乎完全一致。在复杂的肝脏几何模型上，虽然两方法使用的单元平均尺寸略有不同（CutFEM为0.0059953，FEM为0.00635283），位移结果仍表现出良好的一致性。关键的计算效率比较显示：在相近精度下，FEM需要983.41毫秒求解系统方程，而CutFEM仅需429.85毫秒，加速约2.3倍。
研究者将半径为0.7的球形包涵体（模拟肿瘤）隐式嵌入尺寸为6×2×2的软组织模型中，针以3.5度倾角插入。包涵体与基体组织的杨氏模量比被设为1、2、4和8。力-位移曲线显示，当针尖到达组织表面时交互力开始出现并持续上升至穿刺强度阈值，随后针穿刺进入组织。随着针尖接近包涵体，不同刚度比下的力曲线差异显著增大——刚度比越大，交互力越高。当针尖位移达到3时开始退针，交互力变为负值，直至针完全退出后力归零。网格细化研究（从498到3183个自由度）表明Lagrange乘子法给出了收敛的结果。
在肝脏针插入仿真的对比中，CutFEM隐式施加的Dirichlet约束和针尖-隐式表面间的接触交互均与标准FEM的显式处理结果吻合良好。针插入和退针过程中的位移测量曲线在两方法间高度一致。
该方法被应用于深部脑刺激（Deep Brain Stimulation， DBS）手术中电极导线的植入模拟，并考虑了开颅后脑脊液泄漏引起的脑移位（brain shift）现象。模拟结果表明，脑移位导致丘脑底核（subthalamic nucleus， STN）靶点产生约1.1 cm的稳定位移；随后插入导管（cannula）时，由于摩擦相互作用，STN靶点位移进一步增加；当导管开始退管后，STN位移逐渐回降并稳定在脑移位后的位置附近。这一结果验证了CutFEM在非共形网格上处理脑组织与颅骨、导管及电极之间复杂交互的能力。
本研究提出的共旋CutFEM在保持标准FEM精度的前提下，实现了离散化与几何描述的彻底解耦。该方法在科学层面证明了非共形背景网格与隐式几何描述相结合时，最优收敛率仍然可达；在应用层面，它极大简化了患者个体化手术模拟的预处理流程，将预处理成本从繁琐的网格重建降低为仅需将几何表面嵌入背景网格。对于实时手术训练与术前规划系统而言，这意味着可以在使用较粗计算网格的同时保持几何细节的准确性，从而同时满足精度与速度要求。
本研究的方法论贡献体现在多个层面：多层嵌入算法有效解决了复杂曲率表面与粗网格相交时的无效切割问题；移动节点稳定化策略简单高效地消除了极小切割体积对系统稳定性的损害；模板化子四面体嵌入方案将切割单元的积分统一为固定模板的旋转与映射操作，高度利于代码实现。此外，在GPU并行化的框架下，所提出的约束求解三步策略保证了实时性。这些方法创新共同构成了一套从几何处理到力学求解的完整流程，为CutFEM在生物医学模拟领域的实用化铺平了道路。
作者明确指出当前实现仍以线性单元为主，且嵌入子单元也为线性，因此在精确捕捉曲面几何方面仍存在局限。未来工作将探索使用高阶单元或等几何分析（isogeometric analysis）以进一步提高几何逼近精度。此外，文中也提到了针对该方法严谨数学分析的后续研究方向。