徐衍睿等：面向视网膜脱离手术的硅油填充模拟 1237 2.1水-硅油两相液体耦合及临界面精细化模拟 度h内的所有邻居粒子j的值A近似为： 对眼内的水和硅油两相流体填充过程中不同 A=∑A iW(xjh) (3) 性质流体交互进行分析，需构建基于粒子的离散 化模型，用于模拟多相流体整体运动，综合分析多 A的梯度可直接由光滑核函数梯度近似得 相流物理特性，建立多相流体运动模型. 到，即： 2.1.1纳维-斯托克斯方程与SPH离散化 (4) 纳维-斯托克斯方程(Navier--Stokes Euqation, %=∑Agw- N-S)为描述黏性不可压缩流体运动的物理方程， 由于光滑核函数二阶导数计算会产生高阶数 可以被表示为7： 值误差不稳定现象，故A的拉普拉斯算子近似可 Dy 采用人工形式， PDi =P8-Vp+uvv (1) 2A,=2+2∑m-A)w (5) 式中，P为流体密度，为流体速度场，p为压强，g为 jPk+0.01h2 体外力，如重力加速度，“为动力学黏性系数，根据 式中，x=-,W=Wx-x,h,d表示仿真物理 不同流体物理性质而不同，D/Dt=(8/t+v.V)被称 空间维度，三维空间为d=3.0.01h2作为极小值用于 为物质导数(Material derivative),描述空间中运动 保证分母非0.将SPH方法应用于N-S方程求解， 物质点物理量随时间变化率，72为拉普拉斯算子， 可对流体运动进行数值模拟.流体在空间上被离 表示物理量梯度的散度.式(1)为描述不可压缩流 散化为流体粒子，粒子所处区域的密度可由其邻 体运动的动量方程，等号左边为流场中某一质点 居粒子所具有的质量和密度近似为： 的速度变化率，右边三项表示这一变化率同三个 力有关系，分别是体外力，压强力，以及黏滞力. A=∑，W=∑，m (6) 在基于SPH方法的流体模拟中，模拟区域内 而N-S方程中的梯度、拉普拉斯算子可由式 任意物理场A在坐标x处的取值可以根据狄拉克 (4)和式(5)进行计算，本文基于流体不可压缩特 函数性质，被近似为P: 性，通过SPH近似密度与流体静态密度的压缩 A()=∫A(r)Wx-x,h)dV (2) 比，利用压强力抵消数值计算过程流体压缩性， 其中，W为一归一化类高斯核函数，称为光滑核函 使用隐式不可压缩SPH算法ISPH进行流体 数(Smoothing kernel function),其中h被称为支持 模拟 半径(Supporting radius),本文选取三次样条函数 2.1.2水-硅油两相液体耦合 (Cubic spline kernel)P作为光滑核函数，dV'为对应 为实现多相流体相间交互动态效果，需提取 的体积积分变量.SPH方法将经典物理学中的连 多相流体交互界面位置.本文采用颜色场方法1 续介质离散化为宏观质点，如图4所示 通过设定阈值提取边界，实现对于多相流的界面 效果控制 界面张力的定义如下： Fi-5akn (7) 式中，α为张力系数，表示促使界面收缩程度，k 为交互界面曲率，n为界面法线，6=∑W为数 量密度.该界面力力图平滑高曲率的界面区 域和最小化总的表面区域.n和x可通过颜色场 定义进行计算，然而传统颜色场c=∑(mjlP)W 受粒子分布影响易出现数值不稳定影响，故使用 图4SPH方法下粒子数值近似示意图 Fig.4 Numerical approximation of particles using smoothed particle 数值密度方法进行平滑，平滑后颜色场的表达式 hydrodynamics method 如下2： 对于处于x:位置的粒子i,其所具有的任意物 (8) 理量A的取值A可以根据其本身取值以及光滑长2.1    水‒硅油两相液体耦合及临界面精细化模拟 对眼内的水和硅油两相流体填充过程中不同 性质流体交互进行分析，需构建基于粒子的离散 化模型，用于模拟多相流体整体运动，综合分析多 相流物理特性，建立多相流体运动模型. 2.1.1    纳维‒斯托克斯方程与 SPH 离散化 纳维‒斯托克斯方程（Navier‒Stokes Euqation， N‒S）为描述黏性不可压缩流体运动的物理方程， 可以被表示为[27] ： ρ Dv Dt = ρg− ∇p+µ∇ 2 v （1） ρ v D/Dt = (∂/∂t+v · ∇) ∇ 2 式中， 为流体密度， 为流体速度场，p 为压强，g 为 体外力，如重力加速度，μ 为动力学黏性系数，根据 不同流体物理性质而不同， 被称 为物质导数（Material derivative），描述空间中运动 物质点物理量随时间变化率， 为拉普拉斯算子， 表示物理量梯度的散度. 式（1）为描述不可压缩流 体运动的动量方程，等号左边为流场中某一质点 的速度变化率，右边三项表示这一变化率同三个 力有关系，分别是体外力，压强力，以及黏滞力. 在基于 SPH 方法的流体模拟中，模拟区域内 任意物理场 A 在坐标 x 处的取值可以根据狄拉克 函数性质，被近似为[28] : A(x) = w A(x ′ )W(x− x ′ ,h)dV ′ （2） dV ′ 其中，W 为一归一化类高斯核函数，称为光滑核函 数（Smoothing kernel function），其中 h 被称为支持 半径（Supporting radius），本文选取三次样条函数 （Cubic spline kernel)[28] 作为光滑核函数， 为对应 的体积积分变量. SPH 方法将经典物理学中的连 续介质离散化为宏观质点，如图 4 所示. W S h j i rij 图 4    SPH 方法下粒子数值近似示意图 Fig.4     Numerical  approximation  of  particles  using  smoothed  particle hydrodynamics method xi i Ai 对于处于 位置的粒子 ，其所具有的任意物 理量 A 的取值 可以根据其本身取值以及光滑长 度 h 内的所有邻居粒子 j 的值 Aj 近似为: Ai = ∑ j Aj mj ρj W(xi − xj ,h) （3） Ai 的梯度可直接由光滑核函数梯度近似得 到，即： ∇Ai = ∑ j Aj mj ρj ∇W(xi − xj ,h) （4） Aj 由于光滑核函数二阶导数计算会产生高阶数 值误差不稳定现象，故 的拉普拉斯算子近似可 采用人工形式[28] : ∇ 2Ai = 2(d +2)∑ j mj ρj (Ai−Aj)· xi j   xi j    2 +0.01h 2 ∇Wi j （5） xi j=xi − xj Wi j = W ( xi − xj ,h ) d d=3 0.01h 2 i 式中， ， ， 表示仿真物理 空间维度，三维空间为 . 作为极小值用于 保证分母非 0. 将 SPH 方法应用于 N‒S 方程求解， 可对流体运动进行数值模拟. 流体在空间上被离 散化为流体粒子，粒子 所处区域的密度可由其邻 居粒子所具有的质量和密度近似为: ρi = ∑ j mj ρj ρjWi j = ∑ j mjWi j （6） 而 N‒S 方程中的梯度、拉普拉斯算子可由式 （4）和式（5）进行计算，本文基于流体不可压缩特 性，通过 SPH 近似密度与流体静态密度的压缩 比，利用压强力抵消数值计算过程流体压缩性， 使用隐式不可压缩 SPH 算法 IISPH[16] 进行流体 模拟. 2.1.2    水‒硅油两相液体耦合 为实现多相流体相间交互动态效果，需提取 多相流体交互界面位置. 本文采用颜色场方法[29] 通过设定阈值提取边界，实现对于多相流的界面 效果控制. 界面张力的定义如下： F s i = 1 δi ακn （7） α κ n δi = ∑ j Wi j n κ ci = ∑ j ( mj/ρj ) Wi j 式中， 为张力系数，表示促使界面收缩程度， 为交互界面曲率 ， 为界面法线 ， 为数 量密度 . 该界面力力图平滑高曲率的界面区 域和最小化总的表面区域. 和 可通过颜色场 定义进行计算，然而传统颜色场 受粒子分布影响易出现数值不稳定影响，故使用 数值密度方法进行平滑，平滑后颜色场的表达式 如下[29] ： c ′ i = ∑ j c j δj Wi j/∑ j 1 δj Wi j （8） 徐衍睿等： 面向视网膜脱离手术的硅油填充模拟 · 1237 ·