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

(上海市計(jì)量測(cè)試技術(shù)研究院, 上海市在線檢測(cè)與控制技術(shù)重點(diǎn)實(shí)驗(yàn)室, 上海 201203)

1 引 言

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

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

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

2.1 NCT理論研究方法

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

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

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

圖2 射流進(jìn)入靜止流體示意圖

10.065 mm<11 mm

(1)

其中φ為噴嘴直徑,mm。

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

2.2 NCT的數(shù)值建模

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

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

2.2.1 固體域設(shè)定

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

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

2.2.2 流體域設(shè)定

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

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

圖6 非接觸式眼壓計(jì)噴嘴氣流速度隨時(shí)間變化曲線

2.2.3 耦合計(jì)算設(shè)定

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

圖7 雙向耦合計(jì)算設(shè)定框圖

3 結(jié)果及分析

3.1 流場(chǎng)分析

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

(2)

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

(3)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表1 不同眼壓下的角膜壓平時(shí)間

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

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

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

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

圖16 溯源實(shí)驗(yàn)裝置[2]

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

圖17 壓平頭

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

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

4 結(jié) 論

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

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

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