重建人眼前房三維模型, 探究不同角膜溫度對前房熱傳遞過程的影響。基于正常志愿者的光學相干斷層掃描影像, 采用UG參數化造型的方法對人眼前房進行三維建模, 并用有限體積分析方法對不同角膜溫度下整個前節的熱傳遞和房水流動過程進行數值模擬。結果顯示角膜溫度對前房熱傳遞和房水流動有顯著影響:前房溫度呈線性分布, 沿瞳孔軸對稱, 溫差越大則對稱性越差; 房水在低溫側上浮, 在高溫側下沉, 呈渦流樣流動, 溫差越大則流速越快; 角膜、晶體和虹膜的熱通量主要受房水流速影響, 流速越快則上述壁面熱通量絕對值越大。基于光學相干斷層掃描影像的人眼前房數值模擬切實可行, 研究結果有助于深入認識角膜溫度差異在前房熱傳遞和房水動力學中的重要作用。
引用本文： 郭競敏, 張虹, 王軍明. 不同角膜溫度下人眼前房熱傳遞的數值模擬. 生物醫學工程學雜志, 2015, 32(6): 1191-1195. doi: 10.7507/1001-5515.20150211 復制
0 引言
人體熱傳遞是生物醫學工程非常重要的研究方向。熱量的傳遞伴隨溫度的變化，影響著組織的新陳代謝和結構功能。目前紅外成像技術可獲得角膜表面溫度，但很難掌握眼內溫度分布和熱傳遞細節。事實上，角膜溫度不僅影響其自身結構功能，亦可改變眼內尤其是前房的熱傳遞過程，甚至影響前房內的房水循環和物質代謝，參與疾病進程。近些年，生物醫學迅速發展，研究手段不斷更新，尤其是數值模擬的方法可對流場進行定量分析，獲得臨床難以觀測的數據。這種方法具有仿真度高、運算快等優點，可為深入認識前房的溫度分布以及相應的房水循環機制提供數理依據。
本研究根據正常志愿者的光學相干斷層掃描影像資料重建人眼前房模型，通過數值模擬的方法對不同角膜溫度下前房流線、溫度分布和熱通量進行分析，現報告如下。
1 方法
1.1 前房影像檢查與參數測量
選取33名正常志愿者，其中男8名，女25名，年齡(56.15±10.84)歲。暗室(＜1 lux)下適應30 s后利用VisanteTM OCT1000型進行前房斷層掃描。采用高分辨模式、眼前段單線掃描程序，間隔45°掃描4個橫截面，待掃描線中點與角膜中央的反光點重合時采集儲存圖像，利用美國邁斯沃克軟件有限公司提供的MatLab(7.10.0.499 R2010a)設計圖像處理軟件進行分析，獲得如圖 1(a)所示的表征前房特征的主要參數：瞳孔直徑(pupil diameter, PD)指瞳孔緣間距離；前房寬度(anterior chamber width, ACW)指雙側鞏膜突間距；中央前房深度(anterior chamber depth, ACD)指角膜后表面至晶狀體前表面的距離；虹膜曲率(iris curvature, I-Curv)指虹膜后表面根部與瞳孔緣頂點連線，虹膜后表面與該直線的最大垂直距離; 房角開放距離500(angle opening distance 500, AOD500)，距鞏膜突500μm處小梁網作垂直于角膜的直線并與虹膜相交，兩點之間的距離為AOD500。取得上述平均值為最終建模參數，具體數值如下：PD為4.593 mm，ACW為11.49 mm，ACD為2.44 mm，I-Curv為0.229 mm，AOD500為0.365 mm。
圖1
前房三維模型建立示意圖
(a)光學相干斷層掃描影像；(b)UG參數化造模；(c)三維網格圖
Figure1. Schematic diagram of 3D anterior chamber model(a) optical coherence tomography imaging; (b) UG model; (c) 3D mesh
1.2 數值模擬原理及邊界條件設置
前房內房水循環遵循基本的物理規律——質量守恒、動量守恒和能量守恒定律。聯合求解上述方程組可以獲得房水流速、流向等信息。由于房水的物理特性與理想狀態下水的特性類似，且流速較慢，故層流條件下可忽略湍流的非定常效應，不需要額外增加湍流模型[1]。通過求解三維不可壓層流條件下的Navier-Stoke方程(1)可獲得前房內部的流場細節。
| $\frac{{\partial {U_i}}}{{\partial t}} + {U_i}\frac{{\partial {U_i}}}{{\partial {x_j}}}=-\frac{1}{\rho }\frac{{\partial P}}{{\partial {x_i}}} + v\frac{{{\partial ^2}{U_i}}}{{\partial {x_j}\partial {x_j}}}$ |
其中Ui為張量形式的速度，表征x、y、z三個方向的速度；P為流體受到的壓力；ρ為流體的密度；v為流體的運動粘性系數。在求解動量方程的基礎上需同時求解能量方程以獲得前房內部的溫度分布。對于不可壓縮流體，在不考慮內熱源的情況下，流體的內能變化等于流體吸收的熱量和流體受力做功產生的熱量之和。熱量主要通過對流、傳導和輻射三種方式進行傳遞，本文中涉及的熱交換過程主要以熱對流的方式進行，如式(2)所示：
| $Q=h \times \left({{T_{\text{f}}}-{T_{\text{w}}}} \right)$ |
其中Q為流經物體表面的總熱量，Tf為房水的溫度，Tw為壁面的溫度，h為對流傳熱系數。這樣通過求解表征質量，動量和能量守恒的方程組便可以獲得給定條件下前房內部流場、溫度場的詳細信息。
采用UG參數化造型的方法等比例還原光學相干斷層掃描影像測量的前房尺寸, 如圖 1(b)所示。角膜后表面、虹膜及晶狀體前表面設定為無滑移壁面。虹膜和晶狀體溫度為37℃[2]。角膜溫度分別為33、35、37、39、41℃。角膜與晶體的曲率半徑設定為6.2、10 mm。AOD500作為出口壓力與眼內壓等同, 設定為15 mm Hg(1 kPa＝7.5 mm Hg)。入口定義為虹膜晶狀體間隙，高度設定為50μm，流向為晶體表面切線方向，流速為房水生成速度3×10-6 L·min-1(質量流量為5.0×10-5 g·s-1)[2-3]。采用ANSYS內置的ICEM軟件進行網格劃分，并對局部網格進行加密，構建前房三維模型，網格總數量達2.53萬，如圖 1(c)所示。
而后，利用ANSYS公司的CFX求解器對網絡節點上的控制方程進行離散求解。該求解器采用隱式有限體積算法，精度高且收斂速度快[4]。收斂殘差標準設定為1×10-6。
2 結果
2.1 前房內部房水的流動特征
圖 2示立位時不同角膜溫度下前房水流線和流速分布。當角膜溫度(33、35℃)低于虹膜晶體溫度(37℃)時，房水沿晶體虹膜表面向上方流動，少部分從上方的房角流出，多數沿角膜后表面向下流，部分從下方的房角流出，部分構成新的房水循環，形成逆時針渦流。當角膜溫度(39、41℃)高于虹膜晶體溫度時房水循環反向，呈順時針渦流。而當角膜溫度(37℃)等于虹膜晶體溫度時房水渦流被打破，少部分從上方的房角流出，多數受重力作用經前房流入下方的房角。房水流速主要受前房溫度差影響，溫差越大，流速越大，當無溫差時房水流速明顯減小。具體數值見表 1。
圖2
不同角膜溫度下房水流線和流速分布
Figure2.
Streamline and velocity distributions of aqueous humor flow at different corneal temperatures
表1
不同角膜溫度下房水流速和流向
Table1.
Velocity and direction of aqueous humor flow at different cornea temperatures
2.2 前房內部房水的傳熱特征
圖 3示立位時不同角膜溫度下前房溫度分布。為便于不同溫度間比較，圖 3定義了一個歸一化的溫度差異函數η(T)=(Tf-T1)/(T2-T1)。T1和T2分別為角膜和虹膜溫度，Tf為房水溫度，η(T)為房水與角膜和虹膜的溫度差異，該值接近1則表示此區域的房水溫度被加熱，反之則被冷卻。當前房存在溫度差時，溫度剖面從高溫區域往低溫區域呈線性分布，且沿瞳孔軸基本對稱，但溫差越大，對稱性越差。當角膜溫度小于虹膜溫度時，房水被虹膜加熱并將熱量傳遞給角膜；當角膜溫度大于虹膜溫度時，房水被角膜加熱并將熱量傳導給虹膜。晶體、虹膜和角膜表面的熱通量比較如表 2所示。溫差越大、房水流速越大，則熱通量絕對值越大。當不存在溫度差時，熱量的傳遞不復存在，各壁面熱通量值均為零。
圖3
不同角膜溫度下前房溫度分布
Figure3.
Temperature distributions of anterior chamber at different corneal temperatures
表2
不同角膜溫度下各壁面熱通量
Table2.
Wall heat fluxes at different corneal temperatures
3 討論
目前研究人眼溫度分布的非侵入式儀器主要為輻射熱測量計和紅外成像技術，通過角膜表面溫度的測量，有關不同疾病、環境、手術條件下角膜溫度差異的研究逐漸深入[5]。但是上述技術難以獲得眼內溫度分布和熱傳遞細節，而角膜差異對前房的熱傳遞尤其是房水循環可能產生重要影響。眾所周知，房水循環在眼內物質代謝、神經遞質傳遞、藥物濃度分布以及眼壓調控中均有重要作用，故研究角膜溫度差異對前房的熱傳遞和房水動力學影響均有積極指導意義。
早在20世紀70年代，研究者嘗試通過計算機數值模擬來探索眼內熱傳遞方式和影響因素。Canning等[6]基于簡單的前房二維模型，在熱力學相關方程的基礎上明確了前房積膿、積血、角膜后沉著物的形態分布規律，并且證實角膜與虹膜溫度差對房水循環有重要影響，房水流動最大速率與溫度差呈正相關。Heys等[7]以人眼超聲影像構建二維前節模型，用27個節點單元初步探索了人眼前節熱傳遞與房水循環的關系。結果顯示溫度差使近角膜面冷而密度大的房水下移，近晶體虹膜面熱而密度小的房水上移，從而驅動房水在前房循環運動。考慮到激光、輻射等情況下角膜溫度變異大，甚至可超過虹膜晶體溫度，而人眼虹膜晶體溫度相對穩定，故本研究角膜溫度的設定范圍更為廣泛，有利于全面掌握房水循環規律。結果顯示房水在前房內總是沿高溫側向上流，低溫側下沉，構成渦流樣對流，當這種溫度差異因為閉瞼等因素不存在時，房水渦流樣對流消失，房水流速明顯減小。該現象也進一步證實前房房水運動主要受溫度差產生的浮升力驅動。
Ng等[8]通過有限元模型探索不同因素對人眼熱傳遞的影響，結果顯示人眼溫度分布從眼后極部向角膜呈現遞減分布且沿瞳孔軸對稱。該模型顯示預測的角膜溫度與實際紅外成像測得的角膜溫度有67%相符度。但該研究未考慮房水循環的影響。而后該研究組將房水循環納入到熱傳遞二維模型進行分析，發現房水流動對眼后節的熱傳遞影響微小，但可改變前房沿瞳孔軸對稱的溫度分布特征，對前房尤其是中央角膜區存在“加熱”作用[9]。隨著計算機技術飛速發展，人眼模型逐漸從二維走向三維，Karampatzakis等[10]利用66 000個四面體元素構建人眼有限元三維模型，發現房水循環對上方前房和角膜具有更明顯的“加熱”作用。以上研究與本研究顯示的前房溫度分布特征類似，前房內房水流動對溫度分布的對稱性有一定影響。本研究進一步對比發現溫度差越大，房水循環越快，對瞳孔軸上方的前房加熱作用越充分，高溫區域上移，低溫區域下移。而當前房溫度差為0時前房溫度梯度的分布特征也不復存在。
熱通量反映房水從接觸壁面吸收熱量(以正值表示)或釋放熱量(以負值表示)的大小。既往研究證明前房溫度差對角膜熱通量存在影響。Kumar等[11]以兔眼為模型研究發現前房溫度差為2℃時角膜平均熱通量為550 W·m-2, 溫度差為4℃時角膜平均熱通量為1 230 W·m-2。而Heys等[7]在溫度差為3℃時算得角膜平均熱通量為800 W·m-2。本研究基于不同溫度差得到的角膜熱通量結果與以上研究有可比性，且進一步獲得了虹膜和晶體的熱通量大小。研究證實房水流速越快，前房熱交換過程越強，熱量越難以富集，壁面熱通量絕對值越大。
4 結論
本研究首次借助光學相干斷層掃描在人眼前房成像方面的優勢，捕捉反映前房特征的關鍵參數，對正常人前房進行三維重建和數值模擬分析，不僅證實了該方法的可行性，也獲得臨床相關研究很難觀測的熱力學現象和規律。隨著計算方法和參數設計細化完善，相信在未來數值模擬作為臨床研究的補充，將被更廣泛地應用在眼部疾病機制探索、輻射安全管理、激光治療和內眼手術等多個領域。
0 引言
人體熱傳遞是生物醫學工程非常重要的研究方向。熱量的傳遞伴隨溫度的變化，影響著組織的新陳代謝和結構功能。目前紅外成像技術可獲得角膜表面溫度，但很難掌握眼內溫度分布和熱傳遞細節。事實上，角膜溫度不僅影響其自身結構功能，亦可改變眼內尤其是前房的熱傳遞過程，甚至影響前房內的房水循環和物質代謝，參與疾病進程。近些年，生物醫學迅速發展，研究手段不斷更新，尤其是數值模擬的方法可對流場進行定量分析，獲得臨床難以觀測的數據。這種方法具有仿真度高、運算快等優點，可為深入認識前房的溫度分布以及相應的房水循環機制提供數理依據。
本研究根據正常志愿者的光學相干斷層掃描影像資料重建人眼前房模型，通過數值模擬的方法對不同角膜溫度下前房流線、溫度分布和熱通量進行分析，現報告如下。
1 方法
1.1 前房影像檢查與參數測量
選取33名正常志愿者，其中男8名，女25名，年齡(56.15±10.84)歲。暗室(＜1 lux)下適應30 s后利用VisanteTM OCT1000型進行前房斷層掃描。采用高分辨模式、眼前段單線掃描程序，間隔45°掃描4個橫截面，待掃描線中點與角膜中央的反光點重合時采集儲存圖像，利用美國邁斯沃克軟件有限公司提供的MatLab(7.10.0.499 R2010a)設計圖像處理軟件進行分析，獲得如圖 1(a)所示的表征前房特征的主要參數：瞳孔直徑(pupil diameter, PD)指瞳孔緣間距離；前房寬度(anterior chamber width, ACW)指雙側鞏膜突間距；中央前房深度(anterior chamber depth, ACD)指角膜后表面至晶狀體前表面的距離；虹膜曲率(iris curvature, I-Curv)指虹膜后表面根部與瞳孔緣頂點連線，虹膜后表面與該直線的最大垂直距離; 房角開放距離500(angle opening distance 500, AOD500)，距鞏膜突500μm處小梁網作垂直于角膜的直線并與虹膜相交，兩點之間的距離為AOD500。取得上述平均值為最終建模參數，具體數值如下：PD為4.593 mm，ACW為11.49 mm，ACD為2.44 mm，I-Curv為0.229 mm，AOD500為0.365 mm。
圖1
前房三維模型建立示意圖
(a)光學相干斷層掃描影像；(b)UG參數化造模；(c)三維網格圖
Figure1. Schematic diagram of 3D anterior chamber model(a) optical coherence tomography imaging; (b) UG model; (c) 3D mesh
1.2 數值模擬原理及邊界條件設置
前房內房水循環遵循基本的物理規律——質量守恒、動量守恒和能量守恒定律。聯合求解上述方程組可以獲得房水流速、流向等信息。由于房水的物理特性與理想狀態下水的特性類似，且流速較慢，故層流條件下可忽略湍流的非定常效應，不需要額外增加湍流模型[1]。通過求解三維不可壓層流條件下的Navier-Stoke方程(1)可獲得前房內部的流場細節。
| $\frac{{\partial {U_i}}}{{\partial t}} + {U_i}\frac{{\partial {U_i}}}{{\partial {x_j}}}=-\frac{1}{\rho }\frac{{\partial P}}{{\partial {x_i}}} + v\frac{{{\partial ^2}{U_i}}}{{\partial {x_j}\partial {x_j}}}$ |
其中Ui為張量形式的速度，表征x、y、z三個方向的速度；P為流體受到的壓力；ρ為流體的密度；v為流體的運動粘性系數。在求解動量方程的基礎上需同時求解能量方程以獲得前房內部的溫度分布。對于不可壓縮流體，在不考慮內熱源的情況下，流體的內能變化等于流體吸收的熱量和流體受力做功產生的熱量之和。熱量主要通過對流、傳導和輻射三種方式進行傳遞，本文中涉及的熱交換過程主要以熱對流的方式進行，如式(2)所示：
| $Q=h \times \left({{T_{\text{f}}}-{T_{\text{w}}}} \right)$ |
其中Q為流經物體表面的總熱量，Tf為房水的溫度，Tw為壁面的溫度，h為對流傳熱系數。這樣通過求解表征質量，動量和能量守恒的方程組便可以獲得給定條件下前房內部流場、溫度場的詳細信息。
采用UG參數化造型的方法等比例還原光學相干斷層掃描影像測量的前房尺寸, 如圖 1(b)所示。角膜后表面、虹膜及晶狀體前表面設定為無滑移壁面。虹膜和晶狀體溫度為37℃[2]。角膜溫度分別為33、35、37、39、41℃。角膜與晶體的曲率半徑設定為6.2、10 mm。AOD500作為出口壓力與眼內壓等同, 設定為15 mm Hg(1 kPa＝7.5 mm Hg)。入口定義為虹膜晶狀體間隙，高度設定為50μm，流向為晶體表面切線方向，流速為房水生成速度3×10-6 L·min-1(質量流量為5.0×10-5 g·s-1)[2-3]。采用ANSYS內置的ICEM軟件進行網格劃分，并對局部網格進行加密，構建前房三維模型，網格總數量達2.53萬，如圖 1(c)所示。
而后，利用ANSYS公司的CFX求解器對網絡節點上的控制方程進行離散求解。該求解器采用隱式有限體積算法，精度高且收斂速度快[4]。收斂殘差標準設定為1×10-6。
2 結果
2.1 前房內部房水的流動特征
圖 2示立位時不同角膜溫度下前房水流線和流速分布。當角膜溫度(33、35℃)低于虹膜晶體溫度(37℃)時，房水沿晶體虹膜表面向上方流動，少部分從上方的房角流出，多數沿角膜后表面向下流，部分從下方的房角流出，部分構成新的房水循環，形成逆時針渦流。當角膜溫度(39、41℃)高于虹膜晶體溫度時房水循環反向，呈順時針渦流。而當角膜溫度(37℃)等于虹膜晶體溫度時房水渦流被打破，少部分從上方的房角流出，多數受重力作用經前房流入下方的房角。房水流速主要受前房溫度差影響，溫差越大，流速越大，當無溫差時房水流速明顯減小。具體數值見表 1。
圖2
不同角膜溫度下房水流線和流速分布
Figure2.
Streamline and velocity distributions of aqueous humor flow at different corneal temperatures
表1
不同角膜溫度下房水流速和流向
Table1.
Velocity and direction of aqueous humor flow at different cornea temperatures
2.2 前房內部房水的傳熱特征
圖 3示立位時不同角膜溫度下前房溫度分布。為便于不同溫度間比較，圖 3定義了一個歸一化的溫度差異函數η(T)=(Tf-T1)/(T2-T1)。T1和T2分別為角膜和虹膜溫度，Tf為房水溫度，η(T)為房水與角膜和虹膜的溫度差異，該值接近1則表示此區域的房水溫度被加熱，反之則被冷卻。當前房存在溫度差時，溫度剖面從高溫區域往低溫區域呈線性分布，且沿瞳孔軸基本對稱，但溫差越大，對稱性越差。當角膜溫度小于虹膜溫度時，房水被虹膜加熱并將熱量傳遞給角膜；當角膜溫度大于虹膜溫度時，房水被角膜加熱并將熱量傳導給虹膜。晶體、虹膜和角膜表面的熱通量比較如表 2所示。溫差越大、房水流速越大，則熱通量絕對值越大。當不存在溫度差時，熱量的傳遞不復存在，各壁面熱通量值均為零。
圖3
不同角膜溫度下前房溫度分布
Figure3.
Temperature distributions of anterior chamber at different corneal temperatures
表2
不同角膜溫度下各壁面熱通量
Table2.
Wall heat fluxes at different corneal temperatures
3 討論
目前研究人眼溫度分布的非侵入式儀器主要為輻射熱測量計和紅外成像技術，通過角膜表面溫度的測量，有關不同疾病、環境、手術條件下角膜溫度差異的研究逐漸深入[5]。但是上述技術難以獲得眼內溫度分布和熱傳遞細節，而角膜差異對前房的熱傳遞尤其是房水循環可能產生重要影響。眾所周知，房水循環在眼內物質代謝、神經遞質傳遞、藥物濃度分布以及眼壓調控中均有重要作用，故研究角膜溫度差異對前房的熱傳遞和房水動力學影響均有積極指導意義。
早在20世紀70年代，研究者嘗試通過計算機數值模擬來探索眼內熱傳遞方式和影響因素。Canning等[6]基于簡單的前房二維模型，在熱力學相關方程的基礎上明確了前房積膿、積血、角膜后沉著物的形態分布規律，并且證實角膜與虹膜溫度差對房水循環有重要影響，房水流動最大速率與溫度差呈正相關。Heys等[7]以人眼超聲影像構建二維前節模型，用27個節點單元初步探索了人眼前節熱傳遞與房水循環的關系。結果顯示溫度差使近角膜面冷而密度大的房水下移，近晶體虹膜面熱而密度小的房水上移，從而驅動房水在前房循環運動。考慮到激光、輻射等情況下角膜溫度變異大，甚至可超過虹膜晶體溫度，而人眼虹膜晶體溫度相對穩定，故本研究角膜溫度的設定范圍更為廣泛，有利于全面掌握房水循環規律。結果顯示房水在前房內總是沿高溫側向上流，低溫側下沉，構成渦流樣對流，當這種溫度差異因為閉瞼等因素不存在時，房水渦流樣對流消失，房水流速明顯減小。該現象也進一步證實前房房水運動主要受溫度差產生的浮升力驅動。
Ng等[8]通過有限元模型探索不同因素對人眼熱傳遞的影響，結果顯示人眼溫度分布從眼后極部向角膜呈現遞減分布且沿瞳孔軸對稱。該模型顯示預測的角膜溫度與實際紅外成像測得的角膜溫度有67%相符度。但該研究未考慮房水循環的影響。而后該研究組將房水循環納入到熱傳遞二維模型進行分析，發現房水流動對眼后節的熱傳遞影響微小，但可改變前房沿瞳孔軸對稱的溫度分布特征，對前房尤其是中央角膜區存在“加熱”作用[9]。隨著計算機技術飛速發展，人眼模型逐漸從二維走向三維，Karampatzakis等[10]利用66 000個四面體元素構建人眼有限元三維模型，發現房水循環對上方前房和角膜具有更明顯的“加熱”作用。以上研究與本研究顯示的前房溫度分布特征類似，前房內房水流動對溫度分布的對稱性有一定影響。本研究進一步對比發現溫度差越大，房水循環越快，對瞳孔軸上方的前房加熱作用越充分，高溫區域上移，低溫區域下移。而當前房溫度差為0時前房溫度梯度的分布特征也不復存在。
熱通量反映房水從接觸壁面吸收熱量(以正值表示)或釋放熱量(以負值表示)的大小。既往研究證明前房溫度差對角膜熱通量存在影響。Kumar等[11]以兔眼為模型研究發現前房溫度差為2℃時角膜平均熱通量為550 W·m-2, 溫度差為4℃時角膜平均熱通量為1 230 W·m-2。而Heys等[7]在溫度差為3℃時算得角膜平均熱通量為800 W·m-2。本研究基于不同溫度差得到的角膜熱通量結果與以上研究有可比性，且進一步獲得了虹膜和晶體的熱通量大小。研究證實房水流速越快，前房熱交換過程越強，熱量越難以富集，壁面熱通量絕對值越大。
4 結論
本研究首次借助光學相干斷層掃描在人眼前房成像方面的優勢，捕捉反映前房特征的關鍵參數，對正常人前房進行三維重建和數值模擬分析，不僅證實了該方法的可行性，也獲得臨床相關研究很難觀測的熱力學現象和規律。隨著計算方法和參數設計細化完善，相信在未來數值模擬作為臨床研究的補充，將被更廣泛地應用在眼部疾病機制探索、輻射安全管理、激光治療和內眼手術等多個領域。
