后房壓力變化對虹膜變形的流固耦合有限元分析

曹月紅 蔡建程 陳雁翎 Volodymyr Brazhenko Ievgen Mochalin

青光眼是全球第一位不可逆性致盲眼病。據(jù)預(yù)計，2020 年全世界青光眼患者將達(dá)到7 960 萬，其中原發(fā)性閉角型青光眼（PACG）患者占26%[1]。房角關(guān)閉的發(fā)病機制中虹膜-晶狀體間隙變化引起房水由后房流向前房的阻力隨著相對性瞳孔阻滯的進(jìn)展而增加，致使前后房壓力差增大，周邊虹膜向前膨隆，堵塞小梁網(wǎng)的濾過部分，最終導(dǎo)致眼壓升高。急性閉角型青光眼早期有效的干預(yù)措施可以挽救患者的視功能[2]。醫(yī)學(xué)臨床觀察及定性分析PACG房角關(guān)閉機制已有一些文獻(xiàn)報道[3-6]，而定量分析虹膜變形及房角關(guān)閉的文獻(xiàn)相對較少。房水流動以及虹膜在房水作用下的變形定量分析，對于深入了解閉角型青光眼的發(fā)病機理、預(yù)防以及治療是非常重要的。

隨著計算流體動力學(xué)（Computational fluid dynamics，CFD）的發(fā)展，數(shù)值求解眼內(nèi)復(fù)雜流動越來越流行。Heys和Barocas[7]計算溫差驅(qū)動自然對流情況下人眼內(nèi)房水流動并解釋Krukenberg條紋形成的機理。郭競敏等[8]基于房水生理學(xué)相關(guān)理論及流體力學(xué)基本原理，對不同寬窄的前房進(jìn)行三維重建，并對不同瞳孔直徑及體位的房水流動進(jìn)行數(shù)值模擬。陳偉等[9]根據(jù)房水生理學(xué)基本理論，在超聲生物顯微鏡掃描結(jié)果的基礎(chǔ)上構(gòu)建了以人眼前房為主的二維模型，計算結(jié)果表明虹膜膨隆較為嚴(yán)重時，前后房的眼壓差值急劇上升。Wang等[10]數(shù)值研究了3D打印人眼模型在房水流場作用下的虹膜結(jié)構(gòu)模型，并用粒子圖像測速（Particle image velocimetry，PIV）技術(shù)進(jìn)行了實驗驗證，表明重力方向和溫度對房水流動影響很大，從而影響虹膜的變形。

在虹膜變形研究方面，陳琛等[11]對實驗動物的虹膜進(jìn)行整體加壓，研究瞳孔阻滯力作用下的虹膜面應(yīng)變同曲率半徑的變化與前后房眼壓強差（后文簡稱前后房壓差）的關(guān)系以及不同大小的瞳孔阻滯力與前后房眼壓強差的關(guān)系。薄雪峰等[12，13]設(shè)計模擬了瞳孔阻滯及虹膜膨隆的實驗系統(tǒng)和方法，并對兔眼虹膜在不同前后房壓差產(chǎn)生的膨隆變形及房角開放度進(jìn)行研究，研究結(jié)果與臨床觀察相一致。張向東等[14]依據(jù)力學(xué)原理建立了虹膜應(yīng)變與前后房壓差之間的函數(shù)關(guān)系，分析虹膜膨隆過程中虹膜變形與前后房壓差之間的定量關(guān)系。

筆者數(shù)值研究正常虹膜-晶狀體間隙（30 μm）及虹膜幾乎粘連（虹膜-晶狀體間隙分別為2 μm、5 μm）情況下的房水流動以及房水壓力作用下的虹膜變形，房水流動計算由CFD得到，并基于單向流固耦合技術(shù)，進(jìn)行虹膜在流場作用力下變形的有限元分析。嚴(yán)格來講房水與虹膜之間耦合為雙向耦合，虹膜在房水壓力作用下變形，虹膜的變形反過來影響房水流動。本研究先計算虹膜在房水壓力作用下變形，其目的是通過定量分析虹膜-晶狀體間隙對前后房壓差及虹膜變形的影響，為PACG瞳孔阻滯發(fā)病機制的探索提供參考。

1 數(shù)學(xué)模型

1.1 房水流動

與血液不同，房水可以建模為Newton流體，所以其流動受控于Navier-Stokes方程[15]。同時由于流速很?。ㄗ畲罅魉僭趍m/s量級），所以可視為不可壓縮流動，其定常流動的連續(xù)方程及動量方程為：

式中v為速度，ρ、p、μ為流體密度、壓力及動力黏度，g為重力加速度矢量。

由于溫度梯度的存在，流體密度發(fā)生變化，產(chǎn)生浮力，形成自然對流。密度變化與Boussinesq假設(shè)近似[16]：

式中ρ0為參考密度、β為體積膨脹系數(shù)，T和Tref分別是溫度和參考溫度。

對于不可壓縮流動，忽略黏性引起的耗散項后，描述穩(wěn)態(tài)溫度場的方程為：

其中cP為定壓比熱容，κ為熱傳導(dǎo)系數(shù)。房水的特性見表1。

表1.房水物性[7]Table 1.Properties of aqueous humor used in simulations

判斷房水流動的流態(tài)（層流或湍流），需使用特征Reynolds數(shù)和Rayleigh數(shù)：

D代表特征尺寸，如前房徑向尺寸（10-2米量級），ν為特征速度，最大在10-3m/s量級。本研究的流場的溫度△T=3 ℃，估算出Re～20<<2 300，Ra～50<<109，所以可以判斷房水流動為層流流動。

前房房角出口為小梁網(wǎng)，本研究建模中使用了0.1 mm厚度的各向同性多孔介質(zhì)區(qū)域代表一部分小梁網(wǎng)，并把其內(nèi)的壓力梯度表示為黏性損失項（Darcy定理）和慣性損失項：

式中α為滲透率，C2為慣性阻力系數(shù)。它們分別由下面公式[17]估計：

其中Dp是多孔介質(zhì)填料床的平均粒徑，它與孔徑d之間的關(guān)系可以由Dp=d表示；ε為空隙率，定義為空隙體積除以填充床體積。

1.2 虹膜變形

虹膜結(jié)構(gòu)材料為超彈性材料，學(xué)者們對其進(jìn)行了測量研究，得出的結(jié)果有一定差異特別是在大變形情況下[11，14，18-20]。當(dāng)虹膜變形較小，虹膜結(jié)構(gòu)建模成線彈性體，參見文獻(xiàn)[10，21]。對于線彈性體，用位移u表示的變形微分方程為：

式中E和ν分別為楊氏模量和泊松比，對于虹膜E=0.027 MPa，ν=0.49[10，21]；f為體積力，在本研究f=ρirisg，虹膜密度ρiris=1 050 kg/m3[22]。

其邊界條件為虹膜根部固定（Dirichlet條件），即

而其他表面為應(yīng)力邊界條件（Neumann條件），即

由線彈線體的本構(gòu)關(guān)系得到應(yīng)力邊界條件：

本研究使用有限體積法對房水流動進(jìn)行計算，使用有限元法進(jìn)行虹膜變形的計算。

2 數(shù)值計算模型

2.1 房水CFD模型

本研究的房水流道幾何模型參考了文獻(xiàn)[7，9]中的模型。圖1顯示了房水流動的模型，其中圖1A圖為正常虹膜形狀，虹膜-晶狀體間隙為30 μm。文獻(xiàn)[7] 認(rèn)為10 μm基本上為正常虹膜-晶狀體間隙的最 小值，本研究在處理虹膜幾乎后粘連的情況，分別設(shè)置虹膜-晶狀體間隙為5 μm及2 μm，見圖1B和C。圖1中眼角膜的曲率半徑為8.3 mm，晶狀體曲率半徑為10 mm，在眼中心軸上晶狀體到角膜的距離為3 mm。后房直徑為12.8 mm，瞳孔直徑為3 mm。圖1C中虹膜的瞳孔端部形狀與圖1A和B稍不一樣，端部曲率半徑從0.2 mm增加至0.29 mm，這主要為下文網(wǎng)格生成考慮。當(dāng)虹膜-晶狀體間隙為2 μm時，對較小曲率虹膜瞳孔端部進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分具有困難。

使用結(jié)構(gòu)化六面體網(wǎng)格進(jìn)行區(qū)域離散。六面體網(wǎng)格的優(yōu)點是正交性好，離散誤差比四面體網(wǎng)格離散誤差要小，且同樣大小使用網(wǎng)格量少。虹膜-晶狀體間隙處使用11 層網(wǎng)格，以便更好地捕捉該處房水詳細(xì)流場特別是物理量的梯度。圖2顯示了房水流動CFD模型，圖中虹膜-晶狀體間隙為圖1B中的5 μm。計算模型共計使用約194萬個單元，文獻(xiàn)[17]的研究表明30 萬六面體網(wǎng)格已經(jīng)可以達(dá)到房水流動仿真的網(wǎng)格無關(guān)性要求。事實上眼內(nèi)房水流動為層流流動，所以網(wǎng)格精細(xì)要求程度比湍流流動低。

圖1.房水流動計算區(qū)域虹膜-晶狀體間隙為30 μm (A)、5 μm (B)和2 μm (C)Figure 1.The computational domain of aqueous humor flow.The iris-lens gap is 30 μm (A),5 μm (B) and 2 μm (C).

圖2.房水流動動力的計算模型Figure 2.The computational model of aqueous humor flow.

前房房角出口位置設(shè)置厚度為0.1 mm的多孔介質(zhì)區(qū)域，代表小梁網(wǎng)的一小部分。完整小梁網(wǎng)分為3 個部分：與前房相接的為葡萄膜網(wǎng)，然后靠外側(cè)是角鞏膜網(wǎng)，占小梁網(wǎng)大部分，最后是連接Schlemm管的鄰管區(qū)薄層組織[23]。前2個區(qū)域小梁網(wǎng)組織孔徑較大，約為100 μm，而鄰管區(qū)組織孔徑 為0.6 μm，小梁網(wǎng)空隙率ε≈0.5。本研究取孔徑d=1 μm，空隙率ε≈0.5 作為多孔介質(zhì)區(qū)域的參數(shù)，代表小梁網(wǎng)一小部分區(qū)域。該區(qū)域網(wǎng)格層數(shù)為11層。如果不設(shè)置該多孔計算區(qū)域，在房角出口處使用壓力出口邊界，計算時會遇到倒流現(xiàn)象，影響計算精度。

房水進(jìn)口設(shè)為速度進(jìn)口條件，進(jìn)口速度由房水分泌量2.4 μm/min除以進(jìn)口面積得到，進(jìn)口處房水溫度設(shè)為37℃。出口設(shè)置為壓力出口，設(shè)置壓力為15 mmHg，即1 995 Pa。壁面設(shè)置為無滑移壁面，角膜壁面溫度設(shè)置為34℃，其余壁面溫度設(shè)為37℃。溫差引起自然對流為前房流動的主要原因。重力加速度為y方向。

在通用CFD軟件ANSYS FLUENT中，采用有限體積法對房水流動的偏微分方程進(jìn)行離散：動量方程和能量方程離散采用二階迎風(fēng)格式，壓力速度耦合迭代采用SIMPLE算法[24]。迭代計算時，監(jiān)測房水的體積平均壓力及速度。迭代3 000 步以后，平均壓力和速度均趨于常數(shù)，表明迭代趨于收斂，再繼續(xù)迭代5 000步結(jié)束。

2.2 虹膜變形FEA模型

在ANSYS結(jié)構(gòu)分析模塊使用Solid186 單元進(jìn)行虹膜結(jié)構(gòu)的網(wǎng)格劃分，共使用了134 080 個單元607 040個節(jié)點。Solid186為三維固體結(jié)構(gòu)單元，具有二次位移模式，單元通過20個節(jié)點來定義，每個節(jié)點有3個沿著xyz方向平動自由度。圖3顯示了虹膜有限元網(wǎng)格。

虹膜根部設(shè)為固支約束，即設(shè)置該表面上節(jié)點xyz方向的位移為0；對于虹膜瞳孔靠近晶狀體側(cè)的邊線分別設(shè)為自由及約束2 種情況，見圖3。在其余表面上，使用表面效應(yīng)單元Surf154從房水壓力場插值得到對應(yīng)單元表面的壓力數(shù)據(jù)，使用ANSYS中的sfe命令進(jìn)行表面載荷施加。然后進(jìn)行虹膜結(jié)構(gòu)變形及應(yīng)力分析。有限元代數(shù)方程組的求解使用預(yù)處理共軛梯度（Preconditioned conjugate gradient，PCG）求解器。

圖3.虹膜有限元網(wǎng)格分析Figure 3.Finite element analysis mesh of the iris.

3 結(jié)果

3.1 房水流場

圖4.房水速度場(mm/s)虹膜-晶狀體間隙為30 μm（A）、5 μm（B）和2 μm（C）Figure 4.The intraocular velocity fields.The iris-lens gap is 30 μm (A),5 μm (B) and 2 μm (C).

在模型z方向作中截面（xy平面），其速度大小云圖分布見圖4?？梢钥吹胶蠓康牧魉俜浅Ｐ。讦蘭/s的量級；在虹膜-晶狀體間隙處速度增加，此處為壓差驅(qū)動的Poiseuille流動；前房房水呈整體的自然對流流動：瞳孔附近的房水溫度較角膜壁面的高，由Boussinesq近似可知房水會因密度減小而形成浮力，到達(dá)上部房角時由于小梁網(wǎng)處的流動阻力較大，絕大部分房水拐彎流到角膜側(cè)，并受冷后下沉，然后在下部虹膜側(cè)受熱后又上升，形成溫差驅(qū)動的自然對流。整體上，房水流動的最大速度在mm/s量級，并且房角下部流出去的房水比上部的多。該房水流動模式與文獻(xiàn)[7，8，25]的結(jié)果相似，表明本研究房水流動模型的正確性。

眼內(nèi)壓力分布如圖5所示，可以看出對于正常虹膜形狀，前后房壓力幾乎相等（后房平均壓力比前房高0.476 Pa），而主要壓力降發(fā)生在房角出口的小梁網(wǎng)（多孔介質(zhì)區(qū)域）處。虹膜-晶狀體間隙為 5 μm 時，后房壓力比前房壓力高31 Pa。對于2 μm的虹膜-晶狀體間隙，前后房的壓力差達(dá)到815 Pa。利用非線性乘冪函數(shù)對前后房壓差與虹膜-晶狀體間隙之間關(guān)系進(jìn)行擬合，得到的函數(shù)見圖6。從圖中可以發(fā)現(xiàn)當(dāng)虹膜-晶狀體間隙小于 5 μm后，前后房壓差開始迅速增加。

3.2 虹膜變形

把前后房的壓力加載到虹膜表面上，先進(jìn)行瞳孔無后粘連時的有限元分析，得到的虹膜變形結(jié)果如圖7 所示，圖中為從后房側(cè)觀察到的變形分布，黑色線框為變形前虹膜位置。從圖7A中可以看出對于正常虹膜形狀，虹膜在房水壓力作用下的變形是非常小的，最大位移發(fā)生在瞳孔附近，其值為51.2 μm，變形方向主要為從后房向前房。

圖5.房水壓力場（Pa）虹膜-晶狀體間隙為30 μm（A）、5 μm（B）和2 μm（C）Figure 5.The intraocular pressure fields.The iris-lens gap is 30 μm (A),5 μm (B) and 2 μm (C).

圖6.前后房壓差與虹膜-晶狀體間隙的關(guān)系Figure 6.The relationship between the anterior and posterior chamber pressure difference and the iris-lens gap.CFD,computational fluid dynamics.

對于瞳孔幾乎后粘連（虹膜-晶狀體間隙為5 μm和2 μm）的情況，虹膜變形分布見圖7B和C，可以看出在后房較大壓力作用下，虹膜整體脫開晶狀體側(cè)向前房變形，對于5 μm虹膜-晶狀體間隙，虹膜的最大變形量達(dá)到2.597 mm，而對于2 μm虹膜-晶狀體間隙，其最大變形量達(dá)到7.479 cm，這在實際中顯然是不可能的。這是因為本研究的計算使用流固單向耦合，即只計算房水壓力作用下的虹膜變形，而事實上房水與虹膜的互相作用是動態(tài)的，當(dāng)虹膜離開晶狀體一定距離時前后房壓差會減小，從而虹膜的變形量會減小。這一動態(tài)過程的模型需要雙向耦合，一般使用順序耦合方法計算，如沈雙等[26]進(jìn)行人內(nèi)耳前庭系統(tǒng)膜迷路的流固耦合數(shù)值模擬時就使用了該方法。

圖7.虹膜瞳孔無后粘連時的變形 虹膜-晶狀體間隙為30 μm（A）、5 μm（B）和2 μm（C）Figure 7.The iris deformation when it is free at the pupil boundary.The iris-lens gap is 30 μm (A),5 μm (B) and 2 μm (C).

圖8.虹膜瞳孔后粘連時的變形虹膜-晶狀體間隙為30 μm（A）、5 μm（B）和2 μm（C）Figure 8.The iris deformation when it is fixed at the pupil boundary.The iris-lens gap is 30 μm (A),5 μm (B) and 2 μm (C).

下面計算假設(shè)瞳孔后粘連，從而在虹膜有限元模型中把后房側(cè)虹膜瞳孔一周的節(jié)點自由度約束為0，如圖3所示。分析得到的虹膜變形結(jié)果如圖8所示，可以看出在30 μm的虹膜-晶狀體間隙的前后房壓差作用下，虹膜的變形量非常小，在微米量級。而在5 μm虹膜-晶狀體間隙的前后房壓差作用下，虹膜已經(jīng)有明顯的膨隆，虹膜中間部位的最大變形量已經(jīng)達(dá)到了0.342 mm。而在極端的2 μm虹膜-晶狀體間隙時，虹膜膨隆程度已經(jīng)超過常規(guī)前房深度（2.5～3.0 mm）。

對以上情況的虹膜最大變形量與虹膜-晶狀體間隙之間關(guān)系進(jìn)行乘冪函數(shù)非線性擬合，得到結(jié)果見圖9。同樣可以發(fā)現(xiàn)當(dāng)虹膜-晶狀體間隙小于5 μm 后，虹膜膨隆程度將加速。如果認(rèn)為虹膜平面距離角膜約為2 mm，如圖1C中所示，那么從圖9 可以看出當(dāng)虹膜-晶狀體間隙在3 μm時，虹膜最大變形量已經(jīng)超過了2 mm，此時虹膜將與角膜接觸，阻礙房水從前房角排出，從而引發(fā)閉角型青光眼。所以可以定量認(rèn)為3 μm虹膜-晶狀體間隙是后房眼壓升高及前房角關(guān)閉的指導(dǎo)參數(shù)。

4 討論

圖9.虹膜最大變形與虹膜-晶狀體間隙的關(guān)系Figure 9.The relationship of the maximum iris deformation with the iris-lens gap.

臨床上認(rèn)為瞳孔阻滯是PACG發(fā)作的一種重要機理，由于虹膜-晶狀體位置不當(dāng)引起房水由后房流向前房的阻力隨著瞳孔阻滯的發(fā)展而增加，致使前后房壓力差增大，虹膜向前膨隆并與角膜接觸，致使房角關(guān)閉，最終導(dǎo)致眼壓升高。對于這一機理的定量分析研究目前尚不多，本研究正是在這一背景下，利用CFD數(shù)值計算了虹膜-晶狀體間隙分別為30 μm、5 μm、2 μm情況下的房水流動，基于單向流固耦合把房水壓力施加到虹膜表面上，利用有限元分析計算了虹膜變形，以研究瞳孔阻滯引發(fā)的閉角型青光眼。

流場的分析表明，對于30 μm正常虹膜-晶狀體間隙，前后房的壓力幾乎相等；對于5 μm和2 μm虹膜-晶狀體間隙（瞳孔幾乎后粘連），后房比前房壓力高出30 Pa和810 Pa。后房的流速很?。考墳棣蘭/s），虹膜-晶狀體間隙的流動為壓差驅(qū)動的Poiseuille流動，而前房流動主要為溫差驅(qū)動的自然對流，最大流速在mm/s量級。利用乘冪函數(shù)對前后房壓差與虹膜-晶狀體間隙之間關(guān)系進(jìn)行非線性擬合，發(fā)現(xiàn)當(dāng)該間隙小于5 μm后，前后房壓差將異速增加。

虹膜變形的有限元分析結(jié)果表明，在瞳孔無后粘連時，在較大的后房壓力作用下虹膜脫開晶狀體向前房變形，最大位移發(fā)生在瞳孔位置，有利于緩解前后房壓差。如果瞳孔發(fā)生后粘連，虹膜中部隆起，在5 μm虹膜-晶狀體間隙的前后房壓差作用下，虹膜已經(jīng)有明顯的膨隆。數(shù)據(jù)擬合結(jié)果表明當(dāng)虹膜-晶狀體間隙小于3 μm后，能形成顯著虹膜膨隆，導(dǎo)致虹膜前表面與角膜內(nèi)表面接觸（房角關(guān)閉），從而引起閉角型青光眼。

本研究的定量研究結(jié)果有助于從生物力學(xué)的角度認(rèn)識人眼虹膜與房水流動的相互作用力學(xué)特性，為閉角型青光眼發(fā)病機理及診斷治療提供指導(dǎo)。已有的研究通常進(jìn)行定性分析瞳孔阻滯，臨床中也做了大量研究，如人工晶狀體置換改變虹膜-晶狀體間隙，使后房壓力得以下降，從而有效緩解高眼壓對眼部的損害。本研究通過改變虹膜-晶狀體間隙進(jìn)行定量研究，通過計算數(shù)據(jù)擬合發(fā)現(xiàn)3 μm虹膜-晶狀體間隙是后房眼壓升高及前房角關(guān)閉的指導(dǎo)參數(shù)，用于指導(dǎo)臨床，適時進(jìn)行手術(shù)干預(yù)。

本研究為房水流場壓力作用下的虹膜變形分析的流固單向耦合分析，研究重點為虹膜-晶狀體間隙對前后房眼壓差以及虹膜變形的影響，虹膜結(jié)構(gòu)建模為線彈性體，其生物力學(xué)參數(shù)選取為目前普遍認(rèn)可的數(shù)值。本研究并未涉及疾病狀態(tài)下個體化數(shù)據(jù)研究，今后將在該方面以及虹膜非線性模型領(lǐng)域進(jìn)一步開展研究。

利益沖突申明本研究無任何利益沖突

作者貢獻(xiàn)聲明曹月紅：課題設(shè)計，收集數(shù)據(jù)，資料分析及解釋，撰寫論文。蔡建程：參與收集數(shù)據(jù)，資料分析及解釋，論文撰寫及修改。陳雁翎：參與收集數(shù)據(jù)，參與修改論文中關(guān)鍵性結(jié)果、結(jié)論。Volodymyr Brazhenko：參與研究，論文修改。Ievgen Mochalin：數(shù)值計算指導(dǎo)，論文修改