ABSTRACT

Objective: To reconstruct the three-dimensional model of glaucoma eyeball and explore the influence of different degrees of pupillary block and trabecular meshwork lesion on anterior and posterior chamber aqueous humor flow. Based on ultrasound biomicroscopy (UBM), the three-dimensional modeling of human eyeball was carried out. The iris-lens space and trabecular meshwork structure were modeled by annular channel. The annular channel was divided into 20 and 32 equal parts respectively. The degree of pupillary block and the degree of trabecular meshwork lesions were indicated by closing different numbers of bisections. The finite volume method was used to simulate the aqueous humor flow. The results showed that there was an exponential relationship between the anterior and posterior chamber pressure difference and the opening of iris-lens space; in addition, the relationship between anterior chamber pressure and the opening of trabecular meshwork was also exponential, but the closure of trabecular meshwork had little effect on the pressure difference between anterior and posterior chamber. The results of this study are helpful to predict the degree of pupillary block and trabecular meshwork lesions in glaucoma patients, and it has reference significance for the disease analysis of glaucoma patients.

Keywords:Glaucoma, Numerical Simulation, Pupillary Block, Trabecular Meshwork Lesions, Aqueous Humor Flow Characteristics

青光眼房水流动的数值模拟分析

王苗¹，陈建^1*，雷苑²，潘雄¹

¹上海理工大学，能源与动力工程学院，上海

²复旦大学附属眼耳鼻喉科医院，上海

收稿日期：2020年7月29日；录用日期：2020年8月12日；发布日期：2020年8月19日

摘 要

重建青光眼眼球三维模型，探究不同程度瞳孔阻滞和小梁网病变对前后房房水流动的影响。基于超声生物显微镜影像，进行人眼眼球三维建模，采用环形通道模化虹膜–晶状体间隙以及小梁网结构，将环形通道分别划分为20和32等分，通过闭合不同数量的等分来表明瞳孔阻滞程度以及小梁网病变程度，应用有限体积法开展房水流动的数值模拟。结果表明，前后房压差与虹膜–晶状体间隙的开度呈指数关系；此外，前房压力与小梁网的开度也呈指数关系，但小梁网闭合程度对前后房压差影响不大。研究结果有助于青光眼患者瞳孔阻滞以及小梁网病变程度的预测，对青光眼患者病情分析具有借鉴意义。

关键词 :青光眼，数值模拟，瞳孔阻滞，小梁网病变，房水流动特性

Copyright © 2020 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

青光眼是全球第二大致盲性眼病 [1]。病理性眼压增高是引发青光眼的主要危险因素 [2]，青光眼眼压增高多数是由于房水流出产生了障碍，如瞳孔阻滞、小梁硬化、前房角狭窄甚至关闭等，破坏了房水循环的动态平衡。瞳孔阻滞 [3] 是原发性闭角型青光眼(primary angle closure glaucoma, PACG)形成的常见原因，是指虹膜组织与晶状体粘连，房水经瞳孔流入前房时受到了额外的阻力，导致虹膜膨隆，当瞳孔阻滞力超过后房房水的压力时，将阻止房水由后房流入前房，虹膜的膨隆程度加剧，造成房角关闭，眼压进一步升高。周文炳 [4] 等临床研究显示，我国38.1%的PACG患者是由单纯的瞳孔阻滞引起的虹膜膨隆，临床表现为急性或亚急性发作。小梁网流出途径作为眼压调节的关键部位，是房水外流的主要途径，担负了人眼60%~95%的房水引流，也是原发性开角型青光眼(primary open-angle glaucoma, POAG)的病变部位。组织学检查中发现，随着年龄增加，小梁网正常细胞成分减少，而衰老细胞增加 [5]，同时伴有小梁网细胞外基质功能低下，这种结构变化会导致小梁间隙变小、房水外流阻力增高，最终导致眼压升高。其他POAG患者小梁网途径的特征包括小梁网硬度增加和Schlemm管内皮细胞层多孔性的降低 [6]。这些病理性改变都会增加POAG患者房水流出的阻力，从而造成眼压升高。综上所述，瞳孔阻滞以及小梁网病变是诱发青光眼的重要原因，因此研究两者对眼内房水流动状态的影响意义重大。

本研究依据超声生物显微镜(UBM)影像重建了人眼眼球的三维模型，采用环形通道模化虹膜–晶状体间隙以及小梁网结构并以环形通道的不同开度代表瞳孔阻滞和小梁网病变程度。应用有限体积法模拟人眼眼球房水流动，获取不同工况下房水的温度场、压力场以及速度场的分布。分析房水流动特性，揭示前后房压力与瞳孔阻滞和小梁网病变的关系。

2. 物理模型

三维人眼物理模型的构建过程如图1所示。首先采集复旦大学附属眼耳鼻喉科医院原发性闭角型青光眼患者的UBM影像，小梁网、虹膜、晶状体、睫状体及角膜的位置如图1(A)所示。其次，进行眼球轮廓提取，得到如图1(B)所示的二维模型图。图中标明模型的大部分参数尺寸，小梁网流出途径模化为四边形，边长分别为136 mm和100 mm。虹膜–晶状体的间隙太小未在图中注明，为4.67 mm。最终的三维模型由二维模型图旋转而得(图1(C))，将用于数值模拟计算。

图1. 三维人眼模型构建过程。(A) 患者UBM图像；(B) 人眼二维轮廓；(C) 人眼三维模型

图2. 虹膜–晶状体通道和小梁网通道的划分

图2为瞳孔阻滞和小梁网病变的模化。研究的人眼三维模型半剖图如图2(A)所示，房水流过的虹膜–晶状体间隙为虹膜–晶状体通道，房水从房角流出途经小梁网的部位为小梁网通道。为了模拟瞳孔阻滞的严重程度，将虹膜–晶状体环形通道沿周向划分为20等分，如图2(B)所示，每两个对称等分为一个对称组(保证前房房水流动的对称性)，用参数 $\alpha$ 表示虹膜–晶状体的开度，定义为开放的对称组数量与总对称组数量的比值。同理，小梁网的病变程度采用小梁网环形通道的周向开度ε和轴向开度θ表示，如图2(C)所示。将小梁网环形通道沿周向8等分，沿轴向(z轴方向) 4等分，ε定义为周向开放的数量与总的周向等分数量的比值，θ定义为轴向开放的数量与总的轴向等分数量的比值，其局部放大图如图2(D)所示。例，当轴向和周向都只有一个等分打开时，小梁网的开度表示为θ = 0.25，ε = 0.125。

3. CFD数值模拟

3.1. 控制方程

房水的流动遵循三个基本控制方程，即：质量方程、动量方程、能量方程。在稳态不可压缩流动情况下，质量守恒方程的形式如式(1)所示：

$\frac{\partial u}{\partial v}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$ (1)

式中，u，v，w是速度矢量在x，y，z方向的分量。该定律可表述为：单位时间内流体微元体中质量的增加，等于单位时间内流入该微元体的净质量。

纳维–斯托克斯方程(Navier-Stokes equation，N-S方程)是描述粘性不可压缩流体动量守恒的运动方程。它的矢量形式为(2)：

$\rho \frac{\text{d}V}{\text{d}t}=\rho g-\nabla p+\mu {\nabla }^{2}V$ (2)

式中， $\rho$ 是流体密度，p是压力，常数 $\mu$ 是动力粘性系数。N-S方程概括了粘性不可压缩流体流动的普遍规律，在流体力学中具有特殊意义。

由于虹膜和晶状体的温度与角膜的温度不同，需要求解能量方程来得到前房的温度分布。能量守恒方程形式如式(3)所示：

$\rho C_{p}v\cdot \nabla T=\lambda {\nabla }^{2}T$ (3)

式中，C_p为定压比热容，λ为导热系数，T为温度。

此外需考虑浮升力模拟温度差对房水流动的影响。本研究选取Boussinesq模型进行近似计算，这种近似只考虑了温度变化引起的密度变化，而忽略了压强变化引起的密度变化。关系式如式(4)所示：

$\rho ={\rho }_0\left[1-\beta \left(T-T_{ref}\right)\right]$ (4)

其中热膨胀率β为定压条件下密度变化率的函数：

$\beta =-\frac{1}{\rho }{\frac{\partial \rho }{\partial T}|}_p$ (5)

通过求解上述质量、动量、能量守恒的方程组，可获得给定条件下眼前节内部房水流场、温度场及压力场的详细信息。

3.2. 网格划分

网格的划分和生成是数值计算的重要环节。一般而言，网格划分越密得到的结果就越精确，但耗时也越多，网格划分粗糙虽然耗时少，但得到的结果不准确，因此选择合适的网格尺寸建立高质量的网格对计算精度和效率有重要影响。计算域的离散方法可分为非结构化和结构化网格，考虑到本研究中三维模型的不规则性，采用非结构化网格对其进行划分。由于晶状体与虹膜之间的间隙太小，为保证质量对其进行了局部加密，如图3所示。计算域网格总数大约为552万。

图3. 计算域网格划分

3.3. 边界条件及算法设置

设置角膜、虹膜和晶状体为固定的刚性无滑移边界 [7]。虹膜和晶状体温度与体温相同，设为37℃。由于角膜与空气接触，其温度通常低于体温，所以角膜的温度设定为34℃ [8]。在正常眼球中房水流量为2.5 ml/min，定义分泌房水的睫状体突为计算域的速度入口，根据模型尺寸的设定，简化入口速度为v = 2e−6 m/s [9]。房水经前房从小梁网流出，穿过小梁网引流至Schlemm管，研究表明，Schlemm管内压力为1200 Pa [10]，故将此出口设置为压力出口 [11]。房水假设为不可压缩牛顿流体，表1所示为房水的物性参数。本研究模拟的是人体平卧状态时的房水流动，向上方向为正，重力加速度的方向设置为向下。在模拟时采用了FLUENT 3D双精度压力求解器对房水流动域进行求解，欠松弛因子和离散格式均采用默认值，为保证计算精度，采用有限体积法对计算区域作离散化处理，动量和能量方程离散采用二阶迎风格式，压力和速度采用SIMPLE算法。

表1. 模型选用的物理参数值 [12] [13]

4. 结果

4.1. 模型验证

图4. 虹膜–晶状体环形通道示意图

为了验证仿真的准确性，将虹膜–晶状体环形通道(图4)的前后压差与Silver [14] 的压差计算公式中的理论计算值进行比较，其压差的计算公式是由连续性方程和动量方程推导而来：

$\Delta P=P\left({\theta }_B\right)-P\left({\theta }_A\right)=\frac{6\mu f\ln \left(\tan \left({\theta }_B/2\right)/\tan \left({\theta }_A/2\right)\right)}{\text{π}{g}^3}$ (6)

表2. 虹膜–晶状体环形通道参数

虹膜–晶状体环形通道的参数如表2所示，将这些参数带入式(6)可得：DP = 35.9 Pa。

压差的仿真结果为DP₁ = 36 Pa。从计算结果看，模拟结果与理论计算结果相差不大，证明模型选取的合理性。

4.2 眼球温度分布

图5中Case 0.1-0.25-0.125 表示虹膜–晶状体通道的开度为α = 0.1，小梁网通道的轴向开度为θ = 0.25，周向开度为ε = 0.125，瞳孔阻滞及小梁网病变程度用此方式表示。模拟结果显示，三种情况下虹膜到角膜的温度都呈梯度下降，靠近晶状体和虹膜的位置温度较高使得房水密度相对较小，靠近角膜位置的房水温度低密度较大，密度差异产生了垂直方向的浮力。密度大的房水由于重力因素向下流动，密度小温度高的房水由于浮力向上流动，形成了梯度状的温度分布。图5所示，改变虹膜–晶状体通道开度和小梁网通道开度对温度分布几乎没影响。因为在整个流域内房水的流速很低，对流换热对温度分布的影响可以忽略不计，导热在传热过程中起着主导作用。温度场的分布与Song [13] 文献中的一致，证明本研究中模型参数选取的合理性。

图5. 瞳孔阻滞和小梁网病变在不同工况下的温度场

图6. 速度云图。(A) 瞳孔阻滞下的房水流动分布；(B) 小梁网病变下的房水流动分布

4.3 房水流动特性

图6为瞳孔阻滞和小梁网病变下的房水流动状态。如图6(A)所示，房水通过睫状体突分泌，流过虹膜–晶状体通道，随着虹膜–晶状体通道开度的减小，房水的流动速度显著增加，图(a)中开度α = 0.1时虹膜–晶状体间隙附近的高速区明显扩大。当房水经瞳孔到达前房时，沿着虹膜表面流入前房角，在虹膜前部时流速较大，在流向房角的过程中速度逐渐减小。房水在瞳孔附近流动时空间体积较大，流速较小，当房水经过瞳孔流向房角时空间体积骤减，导致虹膜前部的房水流速较大。房水流动到房角附近时，一部分房水通过小梁网流出前房，另一部分沿角膜流动，其流速分布与虹膜表面相似，这部分房水与前房中央新的房水混合形成一个循环区域。此外，不同瞳孔阻滞程度下前房的房水流动基本相似，表明瞳孔阻滞对前房的房水流动影响不大。图6(B)所示，小梁网病变程度对后房房水流速几乎无影响，但对虹膜、角膜表面高速区有明显影响，随着开度的减小，虹膜表面的高速区面积增大，角膜表面的高速面积减小，出现此现象的原因与小梁网通道开放位置及开度大小有关。

4.4. 眼内压分析

表3为瞳孔阻滞和小梁网病变对前房压力(AP)和后房压力(QP)影响的计算数据。由表可知，瞳孔阻滞对前房压力大小的影响不大，但随着虹膜–晶状体通道开度逐渐减小时，后房压力不断升高，导致了前后房压差(DP)的增加。此外，数据表明，小梁网不同开度对前后房压差几乎没有影响，但随着开度的减小，前房的压力将增大，当前房房水流出受阻时导致后房房水流动也受到一定程度的影响，造成后房压力的升高，前后房压力同时升高使得压差变化不大。

图7所示为前后房压差和虹膜–晶状体通道开度之间的关系。当虹膜–晶状体通道的开度从1变化到0.6时，压差仅增加约20 Pa，可见，轻度的瞳孔阻滞不会引起前后房压差的显著增加。开度从0.6变化到0.2时，压差增加约100 Pa，仅从0.2到0.1，压差就增加了约130 Pa，瞳孔阻滞会导致前后房压差急剧增加，导致虹膜严重变形，造成房角闭合，从而产生高眼压，将增加PACG的风险。结果表明，前后房压差与虹膜–晶状体通道开度呈指数关系。

图8为前房压力与小梁网开度的关系。可以看出，在开度为α-θ-ε = 1-1-1时前房压力为1912 Pa (14.3 mmHg)，基本为正常眼压。在小梁网周向开度不变的情况下(ε恒定)，前房压力随轴向开度θ的减小而逐渐增大。当小梁网轴向开度不变时(θ恒定)，前房压力与周向开度ε呈指数相关。

表3. 不同工况下前后房压力的计算数据

图7. 前后房压差与虹膜–晶状体通道开度的关系

图8. 前房压力与小梁网开度的关系

5. 结论

本文采用三维建模软件对人眼进行了模型构建，应用有限体积数值方法模拟不同程度的瞳孔阻滞或小梁网病变的房水流动状态，模拟结果表明在不同工况下房水温度场的分布基本一致；瞳孔阻滞对前后房的压差影响较大，虹膜–晶状体通道的开度越小，压差越明显，且前后房压差与虹膜–晶状体通道的开度呈指数关系；小梁网病变对前后房压差的影响不大，但对前房压力和后房压力造成了明显的改变，研究表明，在小梁网通道轴向开度不变时，前房压力与周向开度呈指数相关。本研究获得的关系式可通过测量前后房压力来预测瞳孔阻滞或小梁网病变的严重程度，对青光眼患者病情分析具有借鉴意义。

文章引用

王 苗,陈 建,雷 苑,潘 雄. 青光眼房水流动的数值模拟分析
Numerical Simulation of Aqueous Humor Flow in Glaucoma[J]. 建模与仿真, 2020, 09(03): 264-273. https://doi.org/10.12677/MOS.2020.93027

参考文献