一种硅油填充术眼内硅油乳化过程模拟可视化方法

班晓娟 王佳敏王笑琨* 张雅斓徐衍睿宋重明黄厚斌 朱志鸿

①(北京材料基因工程高精尖创新中心(北京科技大学) 北京 100083)

②(计算机与通信工程学院、人工智能研究院(北京科技大学) 北京 100083)

③(材料领域知识工程北京市重点实验室(北京科技大学) 北京 100083)

④(中国人民解放军总医院海南医院 三亚 572013)

1 引言

孔源性视网膜脱离(Rhegmatogenous Retinal Detachment， RRD)以视网膜裂孔形成为特征，是视网膜神经上皮层和视网膜色素上皮层分离的一类视网膜病变。玻璃体切割(Pars Plana Vitrectomy，PPV)联合硅油填充手术是孔源性视网膜脱离的主要治疗方式。硅油是一种安全有效的眼内填充剂[1]，其理化性质稳定，没有生物毒性相对安全，良好的透光性和屈光性不影响术后对眼底的观察与治疗，有利于视力恢复，强表面张力可以有效封闭视网膜裂孔，并防止眼内低张力液体穿透裂孔造成视网膜2次脱离。但是，硅油长期存在于眼内会发生乳化，乳化部分表面张力降低。填充过少，可能导致剩余硅油不足以完整覆盖视网膜裂孔，降低手术成功率；填充过多，乳化硅油进入前房，会导致一系列眼部并发症，如增生性玻璃体视网膜病变、复发性视网膜脱离、眼内炎症、继发性青光眼等[2]。因此，考虑硅油乳化情况并适量填充硅油，预测硅油乳化状态对手术预后的影响，是提高手术成功率、减少眼部并发症发生的关键。

随着医学与计算科学的交叉学科发展，医学可视化研究与计算机技术辅助医疗诊断已取得了长足进展。蔡轶珩等人[3]利用机器学习方法，将血管像素与非血管像素看作二分类，实现从背景中迅速分割出视网膜血管。陈强等人[4]利用随机森林分类器，通过频谱域光学相干层析技术对视网膜神经纤维层进行分割，以辅助诊断青光眼等疾病。最近徐衍睿等人[5]提出一种针对硅油填充手术的3维建模与可视化方案辅助医生进行手术流程规划，但是该方法的相间压力计算并不精确，并且没有考虑硅油乳化的扩散现象。本文利用医学影像眼球3维模型，提出一种计算机图形学基于物理的流体模拟方法，对眼内硅油填充过程以及乳化过程进行可视化仿真建模，从而辅助医生判断硅油填充量和取出时间，以减少眼部并发症的发生，提升手术预后。

基于物理的流体模拟可视化方法可分为欧拉网格法、拉格朗日粒子法以及混合方法[6]。欧拉网格法将空间离散化为网格，可简化多相流场景相间动量交换的计算过程，便于耦合数据驱动方法以加快计算速度[7]，但容易产生数值耗散；拉格朗日粒子法将流体和刚体离散化为粒子，更适合模拟具有复杂流变性的液体[8]以及流体自由表面的精细动态[9]。其中光滑粒子流体动力学(Smoothed Particle Hydrodynamics， SPH)方法采用粒子积分近似[10]，能够捕捉流体各方面动态细节。Becker等人[11]提出基于Tait方程的弱可压缩SPH(Weakly Compressible SPH， WCSPH)方法以逼近流体不可压缩特性。随后Solenthaler等人[12]提出预测-矫正SPH(Predictive-Corrective Incompressible SPH， PCISPH)方法，通过预测-矫正隐式迭代方案确定流场压强，相较于WCSPH有了较大的效率提升。Ihmsen等人[13]提出隐式不可压缩SPH(Implicit Incompressible SPH， IISPH)方法，进一步提高了求解器的收敛速度与鲁棒性。在此基础上，Bender等人[14]组合无散度求解器和恒定密度求解器，提出了无散度SPH(Divergence-Free SPH， DFSPH)方法，进一步提升模拟性能。对于SPH方法在预测-矫正过程中存在的数值耗散问题，Wang等人[15]提出一种基于朗肯涡模型的湍流精细化方法，通过恢复由于粒子缺乏旋转而损失的能量来增强表面细节。进一步，Liu等人[16]提出了一种基于流函数的涡度恢复方法，将涡量场和速度场联系起来，该方法可以增强现有的涡旋并且捕获额外的湍流。

物理学研究中，“相”指不同物态或同一物态的不同物理性质或力学状态，多相流即系统中同时存在两个或两个以上的流体相态。硅油填充手术过程中，涉及硅油与水两相液体的交互运动，需要考虑不同液体的物理性质，以实现混溶、不混溶等复杂现象的模拟。与欧拉网格法相比，粒子表示方法能更加清晰地表示流体界面，更适合应用于多相流模拟。Müller等人[17]采用扩散方程，使用粒子方法实现了可混溶液体的扩散模拟。为了进一步丰富多相流交互模拟效果，Liu等人[18]引入体积分数方案与SPH方法相结合。这些方案仅利用浓度差驱动不同相之间的扩散效应，忽略了流体运动和作用力分布导致的多相流体的混合与分离。Ren等人[19]将漂移速度引入SPH公式，成功模拟了由于各相之间的相对运动而引起的混合和非混合现象。然而，该方法将质心速度作为混合粒子的速度，这会造成非无散度速度场和粒子内部膨胀，并且该方法与不可压缩求解器不兼容。Yang等人[20]基于亥姆霍兹自由能，扩展了多相流模拟方法，以灵活地捕捉各种混合与分离过程。Yan等人[21]进一步提出多相流与刚体耦合方法，实现可变形体、颗粒材料、多种流体之间的相互作用效果模拟。Jiang等人[22]在体积加权混合速度的基础上，建立了一种新的多相流无散度混合模型，通过自适应单相不可压缩求解器求解混合速度，并确保质量守恒，大大提高了模拟的精度。但是该方法数密度的简单实现限制了可以处理的多相流密度比，并且引入了非守恒的粒子内界面动量源，使得模拟过程不稳定。

本文针对治疗孔源性视网膜脱离的玻璃体切割联合硅油填充术，提出一种眼内硅油乳化过程模拟可视化方法，对硅油填充过程及乳化过程的多相流体动态交互与互溶扩散过程进行模拟。为处理眼内硅油-水两种不同密度流体交互，提出体积不可压缩多相流体计算框架进行耦合模拟；结合内聚力与曲率力描述宏观两相流体耦合中的表面张力作用；引入力平衡散体动力学模型，描述硅油乳化过程中的互溶扩散效应。

2 基于SPH方法的流体模拟基础

2.1 流体力学控制方程

物理上用于描述不可压缩流体动量守恒的纳维-斯托克斯(Navier-Stocks， N-S)方程为

2.2 SPH离散化

3 可混溶多相流体耦合求解器

图1 SPH数值近似示意图

面向硅油填充术中的硅油-水多相流体交互问题，表面张力所导致的交互界面曲率效果，以及硅油乳化效应所产生的多相流体互溶扩散运动，本节基于传统用于单相流体模拟过程的SPH不可压缩流体模拟方法，构建了硅油-水耦合乳化可视化模拟框架，如图2所示。首先，针对传统SPH模拟方法在处理非均一静态密度流场的数值计算误差问题，使用基于体积不可压缩的多相流体压强计算模型，增强数值计算准确性；其次，采用内聚力与曲率力模拟两相流体间表面张力作用；再次，引入力平衡散体动力学模型，模拟硅油-水两相间因乳化产生的互溶扩散效应。

图2 可混溶硅油-水耦合乳化模拟框架

3.1 体积不可压缩的多相流体模拟

时下最先进的无散度隐式不可压缩方法DFSPH[14]通过流体压缩状态描述压强数值。即根据流体密度压缩状态判定压强大小，抵消其他项导致的压缩，使流体整体不可压缩。根据式(2)，有

3.2 相间表面张力建模

为模拟硅油液滴的强表面张力作用，本文引入单相流体间的内聚力(cohesion)以及交互过程中令表面积最小化的曲率力(curvature)[24]。

对于内聚力，可直接表现为同种粒子间的相互吸引作用，在考虑牛顿第二定律的基础上，引入符合高斯曲线形状的吸引规律

即为A=1情况下的SPH方法数值计算过程。位于流体表面的流体粒子由于缺乏邻居，故会产生颜色场偏小问题，通过式(5)可计算颜色场梯度，从而获取界面法线值

3.3 硅油乳化力平衡散体动力学建模

对于硅油和水之间表现为互溶扩散现象的乳化效果，本文采取体积分数方案，令同一离散粒子内部同时存在多相流体。相较于Jiang等人[22]提出的无散度混合模型，本文做了两个改进，分别是粒子内相交换过程和粒子间预测-矫正过程的解耦及隐式计算漂移速度，使得本文方法在满足质量守恒的同时满足粒子内动量守恒，提高了可混溶多相流体交互的稳定性，并且该方法可以便利地集成到各种不可压缩流体求解器中。

对于多相流场景，在每个粒子内部各相混合的状态下，粒子的速度可以描述为

对于动力学过程，本文引入力平衡方程模型[25]，将粒子间相互作用与粒子内物相交换两个过程解耦，显著提高了隐式SPH方法的稳定性，在粒子i中相k的速度变化率可写成

从而可获取粒子间与粒子内部多相流运动交互过程联系。

图3 具有混合相的粒子

4 实验

为验证本文方法的有效性，本节首先进行硅油-水两相流体耦合实验，并验证两相流体耦合中的表面张力作用；然后开展墨滴扩散实验，模拟可混溶多相流相间交互效果；最后对眼球硅油填充进行仿真，进行术后硅油乳化的模拟。本文仿真算法基于C++编写，利用Eigen作为数学计算工具，并使用文献[26]的表面重构方法由粒子生成3D网格。可视化方面，使用OpenGL作为流体粒子可视化工具，显示实时仿真效果，并使用Blender进行离线动画渲染。本文开展实验所用的计算配置是一台128 GB内存工作站，处理器为2个16核CPU，主频为2.30 GHz。

4.1 瑞利泰勒不稳定性实验

本节通过实验验证体积不可压缩SPH方法模拟两相流体耦合的优越性，以及表面张力在流体耦合中的作用。瑞利泰勒不稳定性实验如图4所示。图4(a)—图4(d)使用DFSPH方法进行模拟，图4(e)—图4(h)使用本文方法，图4(i)—图4(l)使用本文方法结合表面张力作用，黄色相为硅油，蓝色相为水，两相流体不可混溶，密度比为0.7:1，硅油相因为密度较小而上浮，从而形成对流效果。

图4 瑞利泰勒不稳定性实验

第1行是使用DFSPH方法进行模拟得到的结果，第2行使用本文所提体积不可压缩多相流体压强计算模型，第3行在体积不可压缩模型计算压强的同时处理了表面张力的作用。对比3种方法不同时刻模拟结果可知，DFSPH密度不可压缩计算方法由于相间压强计算误差产生了严重的渗透现象，本文体积不可压缩计算方法相间压强计算更加稳定，可以保留更多的耦合细节；考虑表面张力作用，硅油相产生表面积最小化趋势，可以有效缓解低密度比流体对流效果的紊乱状态，避免产生硅油小滴。

图5对比了DFSPH和本文方法在该场景下每帧压强求解器的迭代次数。DFSPH方法由于没有表面张力，两相流体对流效果紊乱，易产生硅油小滴，相间压强计算不稳定，迭代次数周期性出现峰值；本文方法考虑了硅油相的表面张力作用，可以稳定模拟两相流体对流效果，在流体块碰撞容器壁后迭代次数趋于稳定。

图5 DFSPH和本文方法迭代次数对比

该实验证明了本文方法可有效模拟眼球内硅油-水两相耦合，并表现了宏观两相流体耦合过程中表面张力的作用。

4.2 墨滴实验

本实验模拟红色墨滴落入水中的扩散过程，展示可混溶多相流体相间交互效果，以辅助硅油乳化过程可视化分析。红色相为墨滴，透明相为水，密度比为5:1，墨滴坠落的过程中会产生相间扩散效果，如图6所示。

该实验表明本文的力平衡散体动力学模型可有效模拟可混溶多相流体相间交互效果，以辅助实现硅油乳化过程可视量化分析。

4.3 硅油注入与乳化实验

本实验模拟硅油填充术的硅油注入过程及术后的硅油乳化过程。该实验通过力平衡散体动力学模型中的扩散系数控制硅油乳化程度，并考虑了表面张力作用，表面张力系数设置为3。图7为术后同一时刻眼内腔硅油乳化不同程度的仿真结果，(a)图扩散系数为0.01，硅油乳化较轻，(b)图扩散系数为0.1，硅油乳化情况较为严重。

该实验验证了本文方法可以对硅油填充术后眼内腔硅油乳化过程进行有效模拟，以辅助术者判定所需合适的硅油注入量及术后硅油取出时间。

5 结束语

玻璃体切割联合硅油填充术在治疗孔源性视网膜脱离手术中日趋流行，然而硅油长期存在于眼内会发生乳化，进而导致角膜带状变性、继发性青光眼、并发性白内障等严重的眼部并发症，所以预测硅油乳化的程度并及时取出硅油，可以在保证最好的治疗效果的同时减少并发症的发生。本文提出一种硅油填充术后眼内硅油乳化过程计算机模拟可视化方法，构建了体积不可压缩多相流体计算框架进行眼内腔硅油-水的耦合模拟，采用内聚力和曲率力相结合处理表面张力作用，并引入力平衡散体动力学模型，实现眼内腔硅油乳化的仿真过程。实验结果表明，本文方法能够较好地实现硅油-水稳定耦合、表面张力、乳化扩散等效果的模拟与可视化。

本文为预测术后眼内腔硅油乳化状态提供了一种有效方法，也给物理模拟辅助眼科学带来了新的研究思路。但是，影响硅油乳化的因素有很多，例如硅油自身的黏度、纯度、挥发性等，以及眼内出血、眼球过度运动等眼内因素，如何考虑实际情况合理设置硅油乳化的难易程度，精确地模拟眼内腔硅油乳化过程，以辅助判定硅油取出时间，是下一步研究的重点与难点。