非接觸式眼壓計壓平角膜過程的數(shù)值模型研究

林正皓, 王 燦, 張忠立, 秦亭亭, 劉貝貝, 馮齊斌, 洪 扁

(上海市計量測試技術研究院, 上海市在線檢測與控制技術重點實驗室, 上海 201203)

1 引 言

非接觸式眼壓計(non-contact tonometer, NCT)是一種通過射出可控空氣脈沖壓平患者眼角膜中央特定面積,再由內(nèi)部計算機進行數(shù)據(jù)處理得到患者眼壓的儀器[1]。相較傳統(tǒng)接觸式眼壓計,NCT有著無需接觸人眼、操作便利、測量時間短等優(yōu)點,已在臨床上逐步成為主要的眼壓測量方式。為了判斷NCT測量眼壓時是否準確,對其工作情況及量值的溯源方法的研究極為必要。目前常用的溯源方法是通過制作標準模擬人眼,以內(nèi)部增壓或固定值眼壓2種形式,復現(xiàn)一個標準的眼壓值,再使用NCT去測量該標準模擬人眼,通過比較測量值和模擬人眼的標準值以實現(xiàn)溯源。針對模擬人眼復現(xiàn)的眼壓值,文獻[2]通過裝置物理接觸模擬人眼并壓平,獲得壓平時的力值后與壓平面積計算得到標準模擬人眼所復現(xiàn)的眼壓值。這個方法可行的關鍵是需要控制壓平模擬人眼面積與NCT空氣脈沖的壓平面積相同,以確保模擬人眼的變形量相同,即模擬人眼剛度的影響量一致。然而,對于NCT空氣脈沖外部流場及壓平角膜面積的研究仍為空缺。由于NCT空氣脈沖的作用時間短(ms級),采用傳統(tǒng)手段難以捕捉角膜被壓平的瞬間及當時的壓平面積,也無法對NCT的工作流場深入研究。隨著計算機技術的高速發(fā)展,計算機仿真已經(jīng)成為研究、設計、工程應用中不可或缺的重要手段[3～5]。通過雙向流固耦合的數(shù)值仿真技術,能夠?qū)崟r分析呈現(xiàn)流體域和固體域緊密結(jié)合的復雜動態(tài)變化問題。

本文通過采用雙向流固耦合的數(shù)值仿真技術,建立了NCT外部流場及角膜的三維瞬態(tài)數(shù)值仿真模型,研究NCT空氣脈沖的流場特征及角膜在NCT作用下的動態(tài)變化過程。

2 理論及數(shù)值建模

2.1 NCT理論研究方法

NCT的空氣脈沖符合流體力學中的射流理論[6],即高雷諾數(shù)流體自孔口、噴嘴等向外射入靜止、同密度的流體中所形成的流動。通過將染色的水經(jīng)過噴嘴射入靜止水槽中(圖1),可以看出射流自噴嘴進入外界靜止流體后,隨著射流行程的增加,射流邊界越寬,流量越大,且存在一個外邊界,在此邊界之外的流體基本不受射流影響。

圖1 染色水射流射入水槽中

對于射流的內(nèi)部(圖2),分為開始區(qū)域和基本區(qū)域2個部分,射流自噴嘴進入外界靜止流體,首先是開始區(qū)域,區(qū)域內(nèi)存在包含核心區(qū)的內(nèi)邊界,其內(nèi)部射流速度基本與噴嘴處一致,速度衰減極小(圖2噴嘴前三角區(qū)域);隨后射流進入基本區(qū)域,此時射流軸心速度隨距離衰減明顯,斷面速度分布逐漸扁平,各個射流斷面的無量綱速度和無量綱距離具有相似性。根據(jù)NCT使用要求,測量時操作距離需保持在噴嘴前方11 mm處,通過圓截面噴嘴的湍流系數(shù)a=0.08及核心區(qū)域長度Sn計算公式可大致計算得被測角膜的位置處于射流的基本區(qū)域Sn。

圖2 射流進入靜止流體示意圖

10.065 mm<11 mm

(1)

其中φ為噴嘴直徑,mm。

在射流的基本區(qū)域內(nèi),每個行程固定的斷面上,射流的速度分布僅僅為距中心線距離r的函數(shù),從而可以推得其接觸角膜后的壓力作用范圍也是一定的,進而可以研究角膜動態(tài)變化過程及壓平時狀態(tài)。

2.2 NCT的數(shù)值建模

對NCT空氣脈沖的流場及其中的角膜建立數(shù)值仿真模型。模型如圖3所示,由固體域(角膜與鞏膜)和流體域(噴嘴內(nèi)流體及外流場)2部分組合而成。

圖3 非接觸式眼壓計數(shù)值仿真模型

2.2.1 固體域設定

固體域即角膜和鞏膜部分。模型參數(shù)以人體統(tǒng)計數(shù)據(jù)的平均值為基準設定[7,9],角膜視為各向同性的彈性體,曲率半徑R1=7.8 mm,密度ρ=1 000 kg/m3,彈性模量為0.86 MPa,泊松比為0.435[6],鞏膜的曲率半徑R2=11.5 mm,角膜和鞏膜的厚度設定為0.5 mm。角膜的外表面為流固耦合計算的數(shù)據(jù)交換面,內(nèi)表面施加了1個法向的壓力邊界條件以模擬人的眼壓,壓力大小取10、20、30、40、50 mmHg(目前眼科光學及儀器示值中仍在使用mmHg，1 mmHg=133.322 Pa,暫沿襲)，均勻覆蓋NCT的量程。。網(wǎng)格總單元數(shù)47 215,節(jié)點數(shù)72 213,如圖4所示。

圖4 固體域模型及其網(wǎng)格

2.2.2 流體域設定

流體域包括噴嘴內(nèi)一段流體及包裹固體域的外流場。NCT空氣脈沖的最大雷諾數(shù)約為Re=2.56×104,流動視為不可壓,計算采用SSTk-ω兩方程模型,介質(zhì)為空氣,其密度ρ=1.225 kg/m3,動力粘度μ=1.789 4×10-5Pa·s。流場模型如圖5所示,左側(cè)小圓柱為非接觸式眼壓計噴嘴內(nèi)的一段流體,直徑為2.4 mm,噴嘴出口至角膜中心的距離為11 mm?？諝饷}沖自流場最左側(cè)進入,經(jīng)過噴嘴噴出進入外流場,最終到達角膜(固體域)。空氣脈沖的速度曲線來源于型號為CoVis-ST的NCT實驗數(shù)據(jù)[10],如圖6所示。網(wǎng)格總單元數(shù)570 309,節(jié)點數(shù)103 862。

圖5 流體域模型、網(wǎng)格、邊界條件及關鍵參數(shù)的定義

圖6 非接觸式眼壓計噴嘴氣流速度隨時間變化曲線

2.2.3 耦合計算設定

流固耦合問題的求解方式主要有兩類：分離求解和直接求解[11,13],本文通過數(shù)據(jù)交換模塊聯(lián)合固體、流體2個求解器分離求解。流體求解器負責流場壓力、速度、溫度等物理量的計算;固體求解器負責位移、應力、應變等物理量的計算;數(shù)據(jù)交換模塊則負責以異步傳遞的形式更新固體、流體2個求解器中的共同變量,流體求解器計算解出力值,將其作為載荷傳遞給固體求解器;固體求解器計算解出位移值,再將其作為載荷傳遞給流體求解器,以達到雙向流固耦合求解。外流場與固體域的接合面為流固耦合計算的數(shù)據(jù)交換面,同時設定了動網(wǎng)格[14],使流體域數(shù)據(jù)交換面在固體域計算過程中變形時,仍然能夠時時“粘附”在固體域上,以保證流固耦合計算順利進行。時間步長由數(shù)據(jù)交換模塊統(tǒng)一控制,設定為0.05 ms,計算停止時間為16 ms。

圖7 雙向耦合計算設定框圖

3 結(jié)果及分析

3.1 流場分析

根據(jù)氣體射流理論,在整個射流范圍內(nèi),任意一點A的壓強等于周圍靜止流體的壓強,則沿軸向任意斷面間的動量守恒,即：

(2)

式中：ρ為流體密度;r0為噴嘴半徑;v0為噴嘴處流速;r為距離中心軸的距離;v為r處的速度;R為當前截面的半徑。因此可以由上式及無量綱速度、距離的相似性計算得到基本區(qū)域內(nèi)射流軸心速度vm的理論值：

(3)

式中：x為點A距噴嘴的距離;a為湍流系數(shù)。取眼壓為10 mmHg時的流固耦合計算結(jié)果。

16 ms時的射流軸心速度理論值和仿真值見圖8所示。從圖8中可以看出仿真結(jié)果很好體現(xiàn)出了在核心區(qū)內(nèi)射流的軸心速度基本沒有衰減的特性;從核心區(qū)過渡到基本區(qū)域時,仿真計算結(jié)果中曲線的變化趨勢也與理論計算的較為相符,但隨后仿真計算的速度曲線迅速衰減至0,這是由于流場中存在角膜這一障礙物導致的。

圖8 數(shù)值仿真與理論計算的軸心速度比較

圖9給出了不同時間下的射流軸心壓力、速度曲線,從軸心速度曲線中可以看出,不同時間下射流都存在核心區(qū),且核心區(qū)內(nèi)流體速度基本沒有衰減,到達基本區(qū)域后,流體遇角膜而速度急劇衰減至0,與理論一致。不同的是,射流初始核心區(qū)域較理論值小,但隨著時間的推移,軸心速度曲線終點不斷向右推移,射流的核心區(qū)域不斷增加,這是由于在角膜變形較小甚至沒有變形時,角膜前端中心位置位于射流的理論核心區(qū)附近,致使射流遭遇障礙物速度迅速衰減至0,核心區(qū)長度受此影響變小;而在16 ms時,眼壓為10 mmHg的角膜受壓變形嚴重,給予射流充分發(fā)展的空間,此時的核心區(qū)長度與理論相符。軸心壓力曲線圖與速度圖類似,射流核心區(qū)內(nèi)壓力分布都為0,與理論一致;到達基本區(qū)域后,壓力迅速升高,是因為射流遇角膜后速度滯止為0,使壓力升高。同樣,由于角膜變形逐漸加劇,壓力曲線的終點隨時間向右推移。

圖9 不同時間的軸心速度、壓力曲線比較

由于噴嘴發(fā)出的空氣脈沖速度隨時間而變大,因此圖9中壓力、速度曲線隨時間推移而增高,速度圖中各個曲線起點的軸心速度即為對應時間噴嘴激發(fā)的空氣脈沖速度。圖10為T=9.35 ms時流場中的壓力、速度分布云圖,可以直觀地體現(xiàn)到圖9的曲線變化趨勢。

圖10 流場壓力、速度云圖 (T=9.35 ms)

圖11為眼壓10 mmHg,不同時間的流場速度灰度云圖,可以清晰看出射流在不同時間下的內(nèi)外邊界。射流的外邊界基本不隨著時間變化,角度約為8.4°;內(nèi)部邊界角度在9.35 ms時約為4.9°,由于16 ms時角膜形變嚴重,核心區(qū)長度增加,內(nèi)邊界角度收窄,約為4.6°。

圖11 不同時間的射流區(qū)域比較

3.2 角膜的動態(tài)變化過程

通過流固耦合的仿真計算完整獲得了不同眼壓下,0～16 ms時間段內(nèi)角膜受NCT空氣脈沖影響的動態(tài)變化過程,以眼壓為10 mmHg的情況為例(圖12),可以觀察到空氣脈沖自噴嘴發(fā)出擊中角膜,并于T=9.35 ms時吹平角膜,之后由于空氣脈沖繼續(xù)變強,角膜逐漸受壓內(nèi)凹這一完整的變化過程。

圖12 角膜受空氣脈沖影響的動態(tài)變化

圖13為角膜在不同時間下的變形及表面壓力分布曲線,變形曲線圖縱軸的變形量指角膜上的點x軸坐標的負值,變形量為正,角膜外凸,變形量為負,角膜內(nèi)凹。首先觀察考慮16 ms以外的結(jié)果,可以看到角膜的中心區(qū)域變形最大,其變形量隨時間的推移而增加,在距離角膜中心約2 mm以外的部分基本不受空氣脈沖的影響而變形。從角膜表面壓力分布圖中同樣可以看到,距離中心越近,角膜受到的壓力越大,隨著距離的增加,壓力迅速衰減,在距離角膜中心約2 mm處壓力衰減至0,并開始產(chǎn)生一小段負壓區(qū)域。隨時間推移,壓力曲線中心峰值不斷增加,但壓力在中心區(qū)域內(nèi)的衰減速度也相應增加,皆在距離角膜中心約2 mm處壓力衰減至0,隨后產(chǎn)生負壓區(qū)域。而對于16 ms的結(jié)果,由于角膜內(nèi)凹變形量過大,帶動邊緣部分變形,使變形范圍有所擴大,同樣壓力曲線的峰范圍也相應擴大,零壓負壓區(qū)域向外移動,整體趨勢仍然與其余結(jié)果一致。

圖13 不同時間下的角膜變形及表面壓力分布

3.3 不同眼壓下的結(jié)果比較

為比較選擇的5個眼壓點的角膜壓平時的狀態(tài),將0～16 ms時間段內(nèi),角膜中心區(qū)域即將開始產(chǎn)生凹陷,且凹陷隨后將不再恢復原位的最終時刻作為角膜的壓平時間進行分析,各個眼壓點的壓平時間見表1所示。由表1可以看到隨著眼壓的增大,空氣脈沖壓平角膜所需要的時間越長,這與眼壓計噴嘴出口空氣脈沖速度隨時間增大的特性相符。

表1 不同眼壓下的角膜壓平時間

圖14與圖15分別為不同眼壓的角膜初始狀態(tài)曲線與壓平狀態(tài)曲線。由于此次角膜設定為線彈性模型,角膜初始狀態(tài)會受到眼壓的影響,圖14中50 mmHg的初始狀態(tài)曲線相比其他幾個小眼壓點有較明顯的凸出。因此,圖15中50 mmHg角膜壓平時的軸向位移與其余曲線相比值更小,是因為其在初始狀態(tài)下比其余眼壓的角膜更為凸出。從圖15中還可以看出,角膜的壓平直徑在不同眼壓點下變化不大,在所選的5個眼壓點下角膜壓平直徑基本處于3.2～3.6 mm的區(qū)間內(nèi),與美國青光眼研究基金會(GRF)的NCT壓平直徑[15]一致。

圖14 不同眼壓的角膜初始狀態(tài)的比較

圖15 不同眼壓的角膜壓平狀態(tài)的比較

為了驗證NCT外部流場及角膜的三維瞬態(tài)數(shù)值仿真模型，上海市計量測試技術研究院研制了一套模擬人眼角膜及配套裝置,并將其經(jīng)過特制的3.6 mm直徑的圓形壓平頭壓平溯源,且將溯源后的裝置在各醫(yī)療機構(gòu)中實地測試,獲得了與臨床結(jié)果相符的結(jié)論[2]。具體溯源過程采用了轉(zhuǎn)矩平衡的方法,通過研制的L型平衡支架，將壓平頭壓平模擬眼時的壓力傳遞至支架另一端的高精度電子天平上,實時、穩(wěn)定地獲得不同眼壓狀態(tài)模擬眼與外部壓平壓力的對應關系。

圖16 溯源實驗裝置[2]

為進一步驗證溯源時壓平直徑對測量眼壓的影響,采用2個不同直徑的壓平頭,以相同方法對裝置進行溯源。壓平頭直徑為3.06 mm和3.2 mm, 3.06 mm壓平頭采用Goldmann眼壓計的標注壓平直徑, 3.2 mm壓平頭采用仿真結(jié)果區(qū)間的最小值,以驗證與直徑為3.6 mm的壓平頭的結(jié)果差異。圖17中從左至右依次為3.06、3.2、3.6 mm的壓平頭。

圖17 壓平頭

圖18為3個不同壓平頭溯源得到的測量眼壓數(shù)據(jù)及變化曲線。由圖18可以看出壓平頭直徑的縮小對測量眼壓有增大的作用,3.06 mm壓平頭溯源得到的測量眼壓與3.6 mm的數(shù)據(jù)相比,其在每個測點上已超出約5.7～9.6 mmHg眼壓值,同時其偏差隨被測眼壓的增大而增大,平均變化率達到31.4%;而3.2 mm壓平頭溯源得到的各個測點的數(shù)據(jù)與3.6 mm的數(shù)據(jù)相比較大,但仍然比較接近,其測量眼壓在最初的測點上偏差最大,達到3 mmHg,隨后其差值隨被測眼壓的增大逐漸減小,最終測點上偏差為1.2 mmHg,平均變化率11.8%。溯源用壓平頭直徑的減小總體上會使測量眼壓增大,然而當直徑從3.6 mm減小至3.2 mm時,測量眼壓隨被測眼壓的增加會逐漸"收斂"到3.6 mm壓平頭的測量眼壓上,而當直徑進一步減小至3.06 mm時,測量眼壓將不斷增大,被測眼壓越大,其誤差越大。驗證了溯源時控制壓平頭直徑的重要性,當壓平頭直徑在3.2～3.6 mm時,對測量眼壓的影響較小,同時驗證了仿真結(jié)果中的壓平直徑區(qū)間。

圖18 不同壓平頭的眼壓數(shù)據(jù)比較

4 結(jié) 論

1) 建立的數(shù)值模型能夠有效研究NCT空氣脈沖壓平不同眼壓的角膜的動態(tài)變形過程,角膜的壓平時間隨眼壓的增加而增加,但不同眼壓的角膜壓平直徑變化不大,角膜壓平直徑3.2～3.6 mm時,結(jié)果與GRF的相關文獻中的數(shù)據(jù)一致,驗證了模型的有效;

2) NCT的流場符合射流理論,存在核心區(qū)、內(nèi)外邊界等特征。外邊界約8.4°,內(nèi)邊界由于角膜受壓變形,其角度隨時間變小;由于流場射流的特性,角膜表面的壓力集中在中間區(qū)域,在距角膜中心2 mm處時,壓力基本衰減至0,隨后會產(chǎn)生一小段負壓區(qū)域;

3) 根據(jù)數(shù)值模型的結(jié)果進行了相關的實驗,當壓平直徑為3.2 mm時,測量眼壓與壓平直徑為3.6 mm時差異不大,偏差隨著被測眼壓的增大而減小,平均變化率為11.8%;當壓平直徑為3.06 mm時,其測量眼壓與壓平直徑為3.6 mm時有較大差異,同時其偏差隨被測眼壓增大而增大,平均變化率達到31.4%,側(cè)面驗證了上述角膜壓平范圍的合理性。