王澍，姜成惠，姚佳烽，朱桂平

(1. 南京航空航天大學(xué) a. 航天學(xué)院； b. 機(jī)電學(xué)院，江蘇 南京 210016；2. 江蘇省口腔醫(yī)院，江蘇 南京 210029)

0 引言

腭咽閉合不全和阻塞性睡眠呼吸暫停低通氣綜合征(OSAHS)都是由于腭咽閉合障礙導(dǎo)致的呼吸道疾病，腭咽閉合障礙會嚴(yán)重影響患者的生活質(zhì)量和身心健康[1]。腭咽閉合功能受到很多的因素影響，一般認(rèn)為上氣道組織結(jié)構(gòu)異常是導(dǎo)致患者的腭咽閉合障礙的重要原因[2]。上氣道結(jié)構(gòu)異常導(dǎo)致的氣道狹窄是引發(fā)OSAHS的病理學(xué)基礎(chǔ)，而軟腭結(jié)構(gòu)缺陷則會導(dǎo)致腭咽閉合不全。除此之外，神經(jīng)障礙和學(xué)習(xí)型障礙也會導(dǎo)致腭咽閉合障礙。

在早期，人們通過從尸體鑄造而成的模型或者動物實(shí)驗(yàn)的方式來研究人體上呼吸道的流動特性[3]。但是，尸源性研究受組織離體的特性變化影響較大，而動物的氣道結(jié)構(gòu)和人類的呼吸道結(jié)構(gòu)也有較大差異，故而以上研究方法仍然有較大的局限性。此后，影像技術(shù)得到了巨大的發(fā)展，研究人員開始使用CT和MRI技術(shù)掃描得到人體在正常生理狀態(tài)下的上氣道結(jié)構(gòu)數(shù)據(jù)，并構(gòu)建對應(yīng)的實(shí)體模型用于相應(yīng)的實(shí)驗(yàn)研究[4]。這種實(shí)驗(yàn)的針對性較強(qiáng)，一般只能用于測試一種參數(shù)。

近些年來，隨著高性能計(jì)算機(jī)的出現(xiàn)以及數(shù)值分析軟件的不斷成熟，計(jì)算流體力學(xué)被廣泛應(yīng)用于上呼吸道的流場模擬，一般的研究流程為先通過CT或MRI等圖像掃描技術(shù)獲取研究對象的氣道模型，再通過仿真軟件計(jì)算所需要的數(shù)據(jù)。于馳[5]分別建立了健康人與OSAHS患者的咽腔模型，對比了兩者呼吸道壓力以及軟腭的位移狀況。CHOULY F等[6]建立了簡化的上氣道模型，并使用流固耦合的方法模擬了在呼氣過程中舌頭和咽腔壁與氣流的相互關(guān)系，研究OSAHS患者打鼾現(xiàn)象的產(chǎn)生原因。使用數(shù)值模擬分析的方法可以有效地克服呼吸道內(nèi)表面復(fù)雜，空間狹小的限制。

1 模型的構(gòu)建與仿真

1.1 上呼吸道三維模型的構(gòu)建

對健康無呼吸道疾病的人體在靜息狀態(tài)下的上呼吸道進(jìn)行CT掃描，得到了人體從聲門到嘴唇和前鼻孔的三維模型，如圖1(a)所示。為了更好地進(jìn)行相應(yīng)的流固耦合計(jì)算，通過測量上呼吸道關(guān)鍵部位的尺寸，重建了簡化的上呼吸道模型，如圖1(b)所示。模型的流體域由口腔、鼻腔和咽腔組成，模型的固體域?yàn)檐涬竦哪Ｐ汀?/p>

圖1 上呼吸道三維模型重建

1.2 仿真模型

1)流體仿真模型

計(jì)算所使用的流體介質(zhì)為常溫下的空氣，密度ρ=1.225 kg/m3，動力黏性系數(shù)μ=1.8×10-5kg/ms。在流體仿真過程中，希望通過計(jì)算得到在某一特定狀態(tài)下的上氣道壓力分布，所以忽略上氣道和軟腭的變形以及相應(yīng)肌肉的調(diào)節(jié)變化，將氣道壁和軟腭視作不可變形的剛性體，將空氣視作不可壓縮的流體且忽略呼吸過程中的溫度變化。所以通過上氣道的空氣要滿足連續(xù)性方程和Navier-Stokes方程，可表示為

(1)

(2)

式中：ρ為密度；t為時(shí)間；ux、uy、uz分別為速度在3個(gè)方向上的分量；V為速度矢量；p為動水壓強(qiáng)；μ為動力黏性系數(shù)；f為單位質(zhì)量的質(zhì)量力。

雷諾數(shù)的計(jì)算公式為

(3)

式中：ρ、v、μ分別為流體的密度、流速和黏性系數(shù)；d為特征長度。根據(jù)公式可知在低速的空氣流動下，腭咽處為層流流動?；谝陨吓袛啵狙芯窟x擇的計(jì)算模型為不可壓縮的黏性湍流Realizablek-ε瞬態(tài)模型。

2)流固耦合仿真模型

通過流固耦合仿真可以直觀地得到在不同條件下的軟腭狀態(tài)，由于軟腭的運(yùn)動會對流場產(chǎn)生不可忽視的影響，所以使用雙向流固耦合對軟腭的變形做分析。在流固耦合計(jì)算中，流體域依然按照之前的流體仿真模型設(shè)置，將軟腭設(shè)置為線彈性可變形的材料，軟腭與流體域的接觸面設(shè)置為流固耦合面。流固耦合區(qū)域需要滿足以下方程：

(4)

以上方程表示流體域與固體域的應(yīng)力、位移的大小相等或者守恒。

3)邊界條件設(shè)置

使用ANSYS workbench中的雙向流固耦合模塊計(jì)算。在仿真計(jì)算中，忽略在實(shí)際的腭咽閉合過程中氣體成分與空氣之間的差異，并且完全不考慮肌肉作用對軟腭變形的影響，只計(jì)算軟腭在氣流流場作用下的變形情況。

在流體仿真模型中，入口邊界為速度入口，出口設(shè)置為壓力出口，其余壁面設(shè)置為固定不可變形壁面。通過調(diào)整入口的氣流速度大小計(jì)算在不同氣流作用下的上呼吸道壓力分布，分別記錄下從0～2 m/s均勻分布的幾組數(shù)據(jù)。為了模擬軟腭在呼吸和發(fā)音下不同時(shí)間的狀態(tài)，還設(shè)置了一組氣流的速度變化為按照正弦波形變化。正方向的速度用于模擬呼氣過程，負(fù)方向的速度則用于模擬吸氣過程。

在本研究中，軟腭的彈性模量大小被考慮為影響軟腭變形的因素，在計(jì)算中使用了5 000 Pa、10 000 Pa、15 000 Pa、20 000 Pa和25 000 Pa的5組材料來分析軟腭的變形與材料彈性模量之間的關(guān)系，泊松比設(shè)置為0.45，材料密度設(shè)置為1 050 kg/m3。

2 仿真結(jié)果與討論

2.1 上呼吸道內(nèi)壓力分布仿真結(jié)果

通過計(jì)算可以得到在不同速度入口條件下的上呼吸道內(nèi)壓力分布，在呼氣狀態(tài)下基于該模型的壓力分布狀態(tài)如圖2(a)所示，可知在軟腭最靠近咽壁處有最大負(fù)壓，在口腔的入口處有最大的正壓。在吸氣狀態(tài)下的壓力分布狀態(tài)如圖2(b)所示?？梢詮臍獾辣诘膲毫Ψ植贾锌闯觯诤魵鉅顟B(tài)下，口腔處的壓力高于咽腔處，口腔的最大壓力分布在口腔入口處。咽腔處有較大范圍的負(fù)壓分布，而且最大的負(fù)壓位置處于軟腭與咽壁的最狹窄處。在吸氣狀態(tài)下，咽腔的壓力高于口腔。

圖2 流體仿真結(jié)果

2.2 軟腭變形隨入口速度和彈性模量的變化

分別計(jì)算入口速度為0.5、1.0、1.5和2.0 m/s時(shí)的上呼吸道的壓力分布，記錄其最大正壓與最大負(fù)壓隨入口速度變化如圖3所示?？梢钥闯錾蠚獾乐械恼龎汉拓?fù)壓以及壓差與入口處的速度呈現(xiàn)某種指數(shù)函數(shù)的相關(guān)關(guān)系。根據(jù)之前得到的呼氣時(shí)的壓力分布狀況，可以得知在軟腭的兩側(cè)的壓差變化隨著入口速度的變化在逐漸增大。

圖3 上呼吸道壓差隨速度變化

使用雙向流固耦合分別計(jì)算從0～2 m/s的幾組穩(wěn)定入口速度下的變形情況，根據(jù)之前得到的軟腭各個(gè)位置處的變形特點(diǎn)，使用軟腭處的最大變形來代表軟腭的變形程度，得到的結(jié)果如圖4所示。從圖中可以看出在其余條件一致的情況下，軟腭處的最大變形與速度也呈指數(shù)相關(guān)關(guān)系。同樣可以在圖中看出軟腭的彈性模量對軟腭的變形有較大的影響，在相同情況下，軟腭的彈性模量越小，軟腭的運(yùn)動幅度越大。

圖4 不同入口速度下軟腭的最大位移

為了更好地模擬人真實(shí)的呼吸和發(fā)音狀況，使用正弦速度函數(shù)代替之前計(jì)算使用的穩(wěn)定速度入口，入口處的邊界速度滿足函數(shù)v=sin(2πt)，每個(gè)周期為1 s，最大峰值速度為1 m/s。對5組不同的彈性模量的軟腭材料進(jìn)行計(jì)算，得到在正弦速度作用下軟腭最大變形隨時(shí)間的變化如圖5所示。呼氣與吸氣狀況軟腭的變形方向相反，在呼氣狀況下，軟腭向咽壁變形；在吸氣狀況下，軟腭向舌根處變形，此處將向咽壁方向的變形定義為正方向，向舌根處的變形定義為負(fù)方向。軟腭的變形與入口速度相關(guān)，軟腭在一個(gè)周期內(nèi)的正向最大變形出現(xiàn)在速度的波峰時(shí)，軟腭在一個(gè)周期內(nèi)反向最大變形出現(xiàn)在速度的波谷時(shí)，且在同樣的速度下，呼氣狀態(tài)的變形量明顯大于吸氣狀況下的變形。

圖5 正弦速度條件下的軟腭最大變形

2.3 不同呼吸方式對軟腭變形的影響

正常情況下，人主要使用鼻呼吸，口腔氣流占總呼吸氣流比例大約為4%～8%。但是在某些情況下，人體會適應(yīng)性地增大口腔氣流的占比，口腔氣流甚至可達(dá)到總氣流的70%[7]。不同的呼吸方式對上氣道的壓降、氣流速度等都有很大的影響，也會影響到軟腭的變形。在以下研究中，將呼吸方式簡化為鼻腔呼氣、鼻腔吸氣、口鼻腔共同呼氣、口鼻腔共同吸氣4個(gè)過程。

圖6(a)為上呼吸道在使用鼻呼氣情況下的壓力分布，與圖2(a)對比可發(fā)現(xiàn)壓力分布大致相同，最大正壓位于口腔處，且最大負(fù)壓位于咽腔最狹窄處。對比使用口鼻同時(shí)呼氣和只使用鼻呼氣的最大正負(fù)壓和壓差如圖7所示。很明顯，在只使用鼻呼氣的情況下，軟腭兩側(cè)的壓差略大于同時(shí)使用口鼻呼氣。圖6(b)為吸氣階段的氣道壓力分布，對比圖2(b)可以發(fā)現(xiàn)其壓力分布也大致相同。

圖6 鼻呼吸情況下的壓力分布

圖7 鼻呼氣與口鼻同時(shí)呼氣壓力對比

對以上兩種呼吸狀況分別做流固耦合分析，使用正弦邊界條件模擬呼吸過程，得到結(jié)果如圖8所示。由圖可知在使用鼻呼吸的情況下軟腭的變形略大于同時(shí)使用口鼻。

圖8 鼻呼吸和口鼻呼吸模式最大變形對比

3 結(jié)語

本文嘗試通過流體分析和流固耦合計(jì)算找到軟腭的受迫運(yùn)動的影響因素。首先基于CT掃描得到的人體上呼吸道結(jié)構(gòu)圖，重建了人體上呼吸道和軟腭的模型，并通過流體仿真和流固耦合方法計(jì)算了不同的入口速度、軟腭彈性強(qiáng)度以及呼吸方式下上呼吸道的壓力分布和軟腭的變形程度，得到以下結(jié)論。

1)不同的呼吸方式對流場的壓力分布有較大的影響。在同等速度下，只使用鼻呼氣軟腭兩側(cè)最大壓差要比同時(shí)使用口鼻呼氣大20%～30%。

2)軟腭兩側(cè)的壓力差大小隨著氣流速度的增大不斷增大；且在呼氣狀態(tài)下，口腔側(cè)壓力高于咽腔側(cè)，在吸氣時(shí)，咽腔側(cè)壓力高于口腔側(cè)。

3)氣流流經(jīng)口腔和咽腔處產(chǎn)生的壓力差是導(dǎo)致軟腭變形的一個(gè)重要原因。因此，在呼氣狀態(tài)下，軟腭向咽后壁移動，而在吸氣時(shí)向舌后根移動。

4)軟腭的變形隨其彈性模量的增加而不斷減小。

通過仿真分析的方式，分析了多種不同的因素對軟腭變形的影響，通過對軟腭變形影響因素的不斷研究，有望能夠更深入地了解腭咽閉合機(jī)制并輔助治療與腭咽閉合相關(guān)的疾病。