氣流作用下氣道組織變形的流固耦合研究1)

秦 源 陳 茜 魏 東 任曉勇 徐光魁,2)

* (西安交通大學航天航空學院多尺度力學-醫(yī)學交叉實驗室,西安 710049)

? (西安交通大學第二附屬醫(yī)院耳鼻咽喉頭頸外科,西安 710004)

** (空氣動力學國家重點實驗室,四川綿陽 621000)

引言

睡眠呼吸暫停是一種常見的呼吸系統(tǒng)疾病,影響著全球多達20%的人口,給人們帶來了健康和經(jīng)濟負擔[1-4].其中,阻塞性睡眠呼吸暫停(obstructive sleep apnea,OSA)是最常見的睡眠呼吸暫停類型[5].OSA 患者在睡眠過程中反復出現(xiàn)上氣道部分(或完全)阻塞,導致氣流陣發(fā)性減少或停止,持續(xù)時間超過10 s.這個現(xiàn)象將伴隨著打鼾以及嚴重的動脈低氧血癥,使得患者睡眠短暫,需要中途醒來以保證氣道的通暢.長期的OSA 將導致血管疾病、抑郁、糖尿病、日間嗜睡等臨床癥狀[4-6],危害著患者的健康與日常安全.目前OSA 最有效的治療方式為手術治療.懸雍垂腭咽成形術(uvulopalatopharynplasty,UPPP)是治療OSA 的常規(guī)外科手術之一[7],UPPP手術通過拓寬患者受限的氣道形態(tài),達到緩解OSA的目的[8].雖然先前的研究發(fā)現(xiàn)[9-10],上氣道最小橫截面積等形態(tài)學變化與手術效果高度相關,但UPPP手術的機制以及患者上氣道阻塞如何確切定位尚不清楚,導致長期手術有效率(大于6 個月)相對較低,僅為40%～ 50%.

近年來,計算流體動力學(computational fluid dynamics,CFD)和流固耦合(fluid solid interaction,FSI)方法被廣泛應用于呼吸系統(tǒng)相關的研究,如口腔矯治[11]、上頜下頜推進手術[12]、UPPP 手術[13]等.與臨床觀察相比,其模擬結果能夠提供詳細的流動特征,比傳統(tǒng)的醫(yī)學成像方法——如計算機斷層掃描(CT)和磁共振成像(MRI)——更具信息性: 由于OSA 患者常出現(xiàn)上氣道狹窄,為了在吸氣時獲得足夠的氣流,上呼吸道狹窄部分需要更多的負壓,反過來,負壓導致氣道的進一步收縮.從流體力學的角度計算上氣道氣流狀態(tài),得到氣道的壁壓分布,并基于此獲得氣道軟組織的彈性變形情況,可以更直觀地得到氣道阻塞的原因和有效的手術矯正方式.

為降低計算成本和復雜性,早期的一些研究[14-19]對上呼吸道進行了簡化建模: 例如,Chouly 等[14]開發(fā)了一個只考慮舌頭的幾何模型來研究呼氣時發(fā)生的氣道阻塞,Tetlow 等[15]將軟腭表示為放置在通道流中的懸臂柔性板,以描述OSA 中存在的動態(tài)不穩(wěn)定性.上述研究將幾何模型進行了簡化,但要想對氣道動力學的具體情況得出結論,就需要一個準確的上氣道成像.最近的研究多基于醫(yī)學成像方法構建精確的上氣道三維形態(tài)[13,20-23].Shang 等[13]以10 例患者特異性數(shù)據(jù)為基礎,構建了精確幾何的上氣道流動的CFD 模型;郭宇峰等[24]通過CFD 方法對比了兒童OSA 患者和正常兒童的上氣道吸氣流場差異,認為吸氣總阻力增大使患兒吸氣費勁,而腺樣體區(qū)巨大的阻力使氣體難以進入鼻氣道.

目前,CFD 方法對上氣道的研究仍具有很大的局限性,具體體現(xiàn)在: 擴張肌的變形功能是決定OSA 嚴重程度和預測手術療效的關鍵因素[25-26],但CFD 方法忽略了這一問題,將上氣道簡單視為一個剛性不可變形的固體導管.Yu 等[21-22]根據(jù)患有OSA和鼻塞的患者建立了上氣道氣流和軟腭相互作用的解剖學準確模型.可觀察到鼻手術后氣道阻力顯著降低,特別是在咽腔擴大和上游阻力減少的腭咽區(qū).此研究僅將軟腭作為研究對象構建成可變形體,未考慮氣道壁的影響和變形.

因此,現(xiàn)今亟需建立可變形固體與流體雙向耦合作用模型,比較懸雍垂腭咽成形手術前、手術后患者的上呼吸道流動特征,闡明手術的生效機制,從而預測手術結果,指導外科的治療.此外,流動模式的選擇對上氣道內(nèi)的氣流場有顯著影響.為簡單起見,一些研究[15-16,19,27-30]假設上氣道內(nèi)為低雷諾數(shù)流動,故選擇層流模型.然而,上氣道存在狹窄的部分,在狹窄處氣流加速并從層流過渡到湍流.正確的湍流模型可以更好地模擬咽道內(nèi)的氣流特性.

目前,由于患者術后隨訪資料稀少、人體上氣道精確幾何結構復雜以及腭咽手術的操作困難,使用CFD 和FSI 方法對于上氣道研究的成果較少,并不足以指導外科醫(yī)生改進手術方案.本文旨在對比CFD 和雙向FSI 方法計算結果的差異,并使用雙向FSI 方法比較懸雍垂腭咽成形手術前后患者的上氣道流動及變形特征,以闡明手術的機制、預測手術效果;同時建立簡化的二維FSI 模型,研究軟腭運動對氣道氣流的影響.

1 方法

1.1 研究對象

經(jīng)西安交通大學第二附屬醫(yī)院倫理委員會批準,在知情同意的前提下,篩選于耳鼻咽喉頭頸外科病院確診OSA 并進行懸雍垂腭咽成形手術治療的患者.從符合上述條件的患者中選取一例手術成功案例A 及一例失敗案例B 作為研究對象,并對二者上氣道內(nèi)的氣流變化進行對比分析.

1.2 懸雍垂腭咽成形術前術后評估

確診OSA 的量化指標為呼吸暫停/低呼吸指數(shù)(apnea-hypopnea index,AHI),即睡眠中每小時呼吸暫停/低呼吸的次數(shù),由AHI 指數(shù)可判斷病情的嚴重程度.其中,當成人睡眠中AHI 指數(shù)在5～ 15 (次/h)范圍內(nèi)時,為輕度癥狀;當AHI 指數(shù)達到15～ 30(次/h)這一范圍時,為中度癥狀;重度OSA 患者的AHI 指數(shù)則大于30 (次/h)[31].對研究對象進行UPPP 手術治療,并在術后進行隨訪和口咽CT 平掃.比較OSA 患者術前、術后多導睡眠監(jiān)測結果和Epworth 嗜睡量表[32]評分,獲得OSA 患者術前、術后AHI 和缺氧程度等指標,驗證手術療效.

1.3 上氣道幾何模型構建

計算域以硬腭和甲狀軟骨板下緣為界,將CT掃描數(shù)據(jù)集導入商業(yè)軟件包(Mimics Medical 21.0)進行三維重建,如圖1(a) 所示.咽部軟組織厚度在2.0～ 4.0 mm 之間[33],本文抽取2 mm 彈性肌肉層進行建模,如圖1(b) 所示,同時,建立內(nèi)部流體域,如圖1(c)所示.在氣道模型中生成非結構化四面體體積網(wǎng)格,其中,流體域沿壁面邊界設置了5 層邊界層網(wǎng)格,并確保第一層y+＜1 以準確模擬邊界黏性層

其中,ρ=1.225 kg/m3為流體密度,y為單元質(zhì)心的壁距離,Uτ為摩擦速度,μ為流體的動力黏度.

1.4 雙向FSI 模型

雙向FSI 模型的流體區(qū)域使用ANSYS Fluent 2021R1 軟件求解如下流體控制方程

其中,下標i和j表示笛卡爾坐標,v為速度,t為時間,p為壓力,f為體作用力.

考慮到吸入氣流的馬赫數(shù)較低(＜ 0.05),將其建模為不可壓縮流體,根據(jù)流體雷諾數(shù)(Re)估值,上氣道內(nèi)的流動為層流或過渡流.選擇標準剪力傳運(standard shear stress transport,SST)的k-ω湍流模型來模擬湍流[34],速度場和壓力場之間的耦合采用couple 算法實現(xiàn).在入口邊界處設置靜壓力條件,在出口邊界處設置半周期正弦函數(shù)質(zhì)量流率出口,其周期為4 s,幅值為300 mL/s.

固體區(qū)域使用ANSYS Mechanical 2021R1 對軟組織進行瞬態(tài)大變形結構分析.模型的固體域的運動方程為

其中σ表示應力,u表示位移,ρ=920 kg/m3為材料密度.軟組織材料采用neo-hooken 超彈性模型模擬本構關系,取喉癌病人手術取出的癌旁組織進行拉伸試驗得到材料參數(shù).雙向的FSI 模型將流體模型與軟組織結構模型相結合,耦合邊界條件為

其中,上標s表示固體區(qū)域的變量,f表示流體區(qū)域的變量,n表示邊界法向.

本文采用順序耦合方法對氣流域與結構域(軟腭)的相互作用進行了數(shù)值模擬(如圖2 所示).分別求解兩組控制方程,并通過ANSYS Workbench 中的system coupling 接口交換位移與壓力信息,得到了流動和氣道變形的耦合效應.

圖2 FSI 計算流程Fig.2 Flow of FSI calculation

1.5 構建軟腭運動二維FSI 模型

除了氣道的阻塞塌陷,軟腭運動也是改變氣道流體流動狀態(tài)的重要因素.本文基于上氣道真實結構建立了軟腭的二維簡化模型,如圖3 所示,以探討軟腭彈性模量對氣道流體流動狀態(tài)的影響.FSI 模型的計算流程如1.4 節(jié)所述,流體區(qū)域進出口之間壓差定義為9 c m H2O,軟腭分別采用彈性模量為1.5 MPa,1 MPa 和0.5 MPa 的線彈性材料.同時在軟腭和其他軟組織之間定義基于無摩擦的接觸作用,在接觸模型公式中加入1 mm 的幾何偏移量,以防止模型在大變形情況下網(wǎng)格互相穿透.

圖3 基于醫(yī)學圖像的二維軟腭建模Fig.3 Two-dimensional soft palate model

2 結果

2.1 三維精確氣道模擬

案例A 和案例B 的患者在術前、術后經(jīng)多導睡眠監(jiān)測和Epworth 嗜睡量表評估,得到如表1 和表2 所示的OSA 嚴重程度指標.

表1 術前OSA 指標Table 1 Pre-operative OSA index

表2 術后OSA 指標Table 2 Post-operative OSA index

其中案例A 的OSA 程度有所改善,而案例B 的OSA 程度加重,可看作失敗案例.經(jīng)FSI 分析,案例A 術后的氣流流動及氣道變形情況與術前有顯著區(qū)別,而案例B 的手術治療對氣道狀態(tài)沒有明顯影響,其壓力分布情況較術前更加不利.

2.1.1 氣流流速及流線分布

圖4 為案例A 和案例B 氣道內(nèi)流場的速度及壓力分布情況.圖4(a)和圖4(b)分別展示了兩案例術前、術后的速度及流線分布.可以看出,當氣流進入收縮區(qū)域后出現(xiàn)高流速.在通過最窄的部分后,流體會減速并與壁面分離,從而形成比較強烈的射流.根據(jù)文丘里效應,由于流率恒定,在流道橫斷面收縮的區(qū)域,流體會被加速,其流速與過流斷面大小成反比,而高流速流體附近會產(chǎn)生低壓,吸附周圍組織使其閉合.而由于軟組織的閉合,氣道更加狹窄,進一步加速了氣流,形成了惡性循環(huán),產(chǎn)生了局部的高流速、低壓力區(qū)域.因此,患者氣道內(nèi)氣流流速激增,對周圍的軟組織產(chǎn)生吸附作用,這將引起軟組織的塌陷,形成呼吸阻塞.其中,圖4(a)所示的案例A 術前最大流速為32.9 m/s,術后最大流速下降為22.4 m/s;圖4(b)所示的案例B 術前、術后氣道內(nèi)的最大流速值分別為21.3 m/s 和22.4 m/s,氣道內(nèi)最大流速值無明顯差異.由此可見,對氣道流速分布的改變關系著手術的治療效果.

圖4 流速和壓力分布Fig.4 Velocity and pressure distribution

2.1.2 壁壓分布

圖4(c)和圖4(d) 分別展示了兩案例術前、術后的壁壓分布,收縮區(qū)域的高速氣流附近出現(xiàn)極低的負壁壓.

氣道狹窄的區(qū)域由于流速加快,出現(xiàn)較嚴重的負壓力區(qū),因此被吸附收縮,導致狹窄區(qū)域的橫截面積進一步減小,從而在吸氣過程中逐步加重局部負壓力程度,最終達到一個極低的負壓力值,造成氣道塌陷與呼吸阻塞.

在手術成功的案例A 中,術前術后的壁壓分布呈現(xiàn)出明顯變化,最小負壓由術前的-940.3 Pa 改變?yōu)樾g后的-393.4 Pa;而手術失敗的案例B 中,最小負壓由術前的-423.7 Pa 下降為術后的-785.5 Pa,氣道整體的負壓力程度加重.進出口壓差影響著吸氣時的阻力,更大的壓差表明吸氣的過程受到更大的阻力.由結果來看,OSA 患者吸氣時受到更大的阻力,同時存在一些高負壓的不利部位,而有效的手術可以減輕局部的高負壓和吸氣時所受阻力.

2.2 二維軟腭運動模擬

軟腭區(qū)域在吸氣過程中流場的速度和壓力分布如圖5 所示.在軟腭與扁桃體之間形成的狹小通道處,氣流被加速,氣道壁出現(xiàn)較低的負壓力,同時,由于軟腭逐漸遠離舌,附近氣流逐漸形成輕微的擾動,E=1 MPa 時,軟腭在氣道中出現(xiàn)明顯撲動,這將導致氣流的不穩(wěn)定以及打鼾的現(xiàn)象.軟腭的彈性模量E進一步影響著流場內(nèi)流速與壓力的分布和極值,如表3 所示.在這一范圍內(nèi),隨著軟腭彈性模量的降低,氣道中最大流速降低,最低負壓值升高,即剛度更小的軟腭將使得氣道內(nèi)的流動狀態(tài)更有利于呼吸.

表3 不同軟腭彈性模量E 對應的流速和壓力極值Table 3 Maximum velocity and minimum pressure corresponding to different values of E

圖5 軟腭區(qū)域流速和壓力分布Fig.5 Velocity and pressure distribution in the soft palate region

3 討論

3.1 CFD 與FSI 計算結果的對比

對同一氣道模型分別采用CFD 方法和FSI 方法計算.由于非線性接觸公式和咽壁的大變形,FSI 模型通常更容易出現(xiàn)穩(wěn)定性和收斂問題,模型的接觸方式以及流體域動網(wǎng)格的更新方式都影響著模型的收斂性,同時,相比于CFD 方法,FSI 方法需要在流固耦合截面的邊界條件上達到收斂,計算效率較低.但CFD 方法忽略了氣道的彈性變形,僅保留流體計算域,FSI 方法則考慮了氣流與氣道軟組織的相互作用.二者計算結果如圖6 所示.定義氣道入口與出口之間的壓降為 ΔP,表4 比較了兩種方法在計算壁壓力時產(chǎn)生的差異.

表4 CFD 和FSI 計算的壁壓力結果Table 4 Calculation results of wall pressure by CFD and FSI methods

圖6 CFD 和FSI 計算結果對比Fig.6 Comparison of CFD and FSI calculation results

FSI 方法計算得到的氣道壁最小壓力值略小于CFD 方法的計算結果,氣道入口與出口之間的壓降明顯大于CFD 的計算結果.由圖6 可知,CFD 方法與FSI 方法得出的氣流流速與壓力分布存在較大差異.產(chǎn)生這種差異是由于雙向FSI 計算考慮了氣流對氣道產(chǎn)生的負壓力所造成的氣道組織彈性變形,在流速較高的區(qū)域,由于壓力較低,氣道壁收縮形成了更狹窄的區(qū)域,使得流速進一步加大,由圖6(a)可以看出,FSI 計算結果中的高流速區(qū)域更大;CFD 則不考慮氣道收縮在吸氣過程中帶來的阻力,因此計算結果中,高流速區(qū)域更小,吸入等質(zhì)量氣體所需要的壓差也更小.

依據(jù)2.1 和2.2 節(jié)所述的計算結果,可知手術療效與氣流流動狀態(tài)的改變緊密相關,因此,僅根據(jù)CFD 方法得到的計算結果探討影響手術成功與否是不準確的,氣道軟組織的彈性變形在氣道阻塞的計算中是不能忽略的.

3.2 流速和壓力分布與手術療效的關系

目前UPPP 手術主要目的是拓寬氣道受限的橫截面積,但僅僅靠氣道橫截面積的變化來判斷手術療效是不夠的.圖7 展示了案例A 和案例B 術前術后的氣道橫截面積,在手術成功的案例A 中,術前氣道的最小橫截面積為2.53 mm2,如圖7(a)所示,手術拓寬了氣道的C段,然而D段橫截面積變小,最小橫截面積為3.40 mm2,與術前相比未發(fā)生明顯變化;最終案例A 表現(xiàn)為術后的AHI 指數(shù)下降了31,手術有效.這說明橫截面積并不是判定UPPP 手術成功與否的關鍵因素,只簡單地考慮橫截面積的影響或許是UPPP 手術失敗率高的原因.

圖7 橫截面積Fig.7 Cross-sectional area

定義AHI 的下降值為 ΔAHI,表5 列出了成功案例A 和失敗案例B 的流動參數(shù).圖7(a) 所示的C段橫截面積擴大后,降低了吸氣過程的阻力,也使得術前的狹窄區(qū)域局部負壓力絕對值下降,也由此減小了氣道收縮進一步帶來的阻力;而圖7(c)所示的E段是術后形成的狹窄區(qū)域,氣流通過該區(qū)域所產(chǎn)生的高速流使得其下方區(qū)域收縮,與術前沒有明顯差別.由此可見,成功手術明顯改變了最大氣流速度和氣道壁上的壓力值分布,手術后進出口之間的壓差更小,計算區(qū)域的負壓程度減輕,最小負壓值增大,有效地減輕了氣道的過度收縮.圖7(c)說明,失敗手術案例B 手術也同樣擴充了氣道F段的截面積,但在E段形成了一個更狹窄的區(qū)域,流動出現(xiàn)最大流速增大、最小負壓值降低的現(xiàn)象,使得OSA 程度進一步加重,表現(xiàn)為AHI 指數(shù)由術前的17.8 上升為術后的52.5,同時缺氧程度加重.

表5 案例A、案例B 的流動參數(shù)Table 5 Flow parameters of case A and case B

圖8 和表6 展示了不同進出口壓差下的流動狀態(tài).進出口的壓差較大時,高流速區(qū)域附近出現(xiàn)極低的負壁壓,這部分氣道將在變形中嚴重塌陷,造成呼吸阻塞.對比A,B 兩案例可以發(fā)現(xiàn),有效的手術減小了進出口之間的壓降程度,當手術降低氣道內(nèi)的最大氣流速度和進出口壓降值,同時增加最小壁壓時,手術表現(xiàn)為成功;當以上指標未發(fā)生明顯改變或反向變化時,手術表現(xiàn)為無效或失敗.

表6 不同 ΔP 對應的流動參數(shù)Table 6 Flow parameters corresponding to different values of ΔP

3.2 軟腭彈性模量對氣流流量的影響

軟腭的運動也是上氣道氣流不穩(wěn)定的驅動因素之一.在入口和出口之間9 cm H2O 壓差的邊界條件下,不同彈性模量的軟腭在變形過程中對氣流狀態(tài)的影響不同.軟腭在吸氣過程中逐漸遠離凸起的扁桃體部位,以保證氣流的進入.當E=0.5 MPa 時,軟腭在吸氣過程中的最大偏移量為1.96 mm,當E=1.5 MPa 時,軟腭的最大偏移量為1.34 mm.當彈性模量E在0.5～ 1.5 MPa 的范圍內(nèi)變化時,從氣流流速與氣道壓力分布來看,更小的剛度對于氣道的影響是有利的,但從整個吸氣過程中氣道出口流出的氣體總質(zhì)量來看,卻并非如此: 當E=0.5 MPa時,吸入氣體的總質(zhì)量為4.85×10-4kg,當E=1.5 MPa時,吸入氣體的總質(zhì)量為5.44×10-4kg.由此可見,盡管更大剛度的軟腭會造成氣流流速和壓力分布情況的不利,但這種情況下氣道仍然能夠吸入更充足的氣體.這是因為,盡管更小的負壁壓會使氣道的受力狀況更加不利,但較大的剛度卻使氣道不易變形塌陷,保證了氣流的進入.

軟腭部分的脂肪含量將影響軟腭的彈性模量大小,當軟腭部分彈性模量較小而引起吸入氣體質(zhì)量不足時,可以通過手術來進行矯正.同時,通過CT 掃描得到的二維實際結構對軟腭區(qū)域的氣流狀態(tài)進行模擬,以指導手術對該區(qū)域結構的改進.

4 總結

本文對人體三維精確氣道進行了準確建模,比較了UPPP 手術前后患者上氣道流動及變形特征,同時建立簡化的軟腭區(qū)域二維模型,研究不同彈性模量下軟腭運動對氣道氣流的影響.過去對于上氣道流動狀態(tài)的研究主要采用CFD 方法,然而該方法忽略了氣道軟組織的變形,FSI 方法則很好地彌補了這一缺陷.通過CFD 和FSI 計算結果的對比,可以看到兩者對流場流速和壁壓分布的計算有一定區(qū)別,因而使用FSI 方法考慮軟組織的影響是非常必要的.

長期以來,UPPP 手術僅簡單地通過擴大上氣道橫截面積來治療OSA,但最小橫截面積大小并不是UPPP 手術成功與否的決定性因素.通過比較手術成功案例和失敗案例的雙向FSI 模型計算結果,可以發(fā)現(xiàn),兩者的主要區(qū)別在于手術對于氣道最大氣流速度、進出口壓降值以及最小壁壓值的改變: 成功的手術使氣道流場中的最大速度值和氣道進出口壓差降低,同時增大了最小壁壓,避免氣道壁出現(xiàn)嚴重的負值;而失敗的手術往往沒有改變或反向改變上述參數(shù).

軟腭的運動對上氣道氣流也存在一定影響.當軟腭的彈性模量E在0.5～ 1.5 MPa 的范圍內(nèi)變化時,盡管剛度更小的軟腭會改善流場內(nèi)壁壓和流速的分布,但這種情況下軟腭相較于大剛度的情況更易變形塌陷,因此從吸入氣體的總質(zhì)量來看,彈性模量為1.5 MPa 的軟腭更有利于氣體吸入.但是,本文工作仍有一定的局限性,如未考慮鼻腔和口腔結構對于氣道流動狀態(tài)的影響等.

綜上所述,本文通過對UPPP 手術前后OSA 患者的參數(shù)分析,發(fā)現(xiàn)UPPP 手術主要改變了流場中的最大速度值和氣道進出口壓差值,這可以幫助臨床醫(yī)生進一步理解UPPP 手術的作用機制;同樣,所建立的FSI 模型也為預測手術效果提供了更好的工具.