崔 巖， 王玉龍， 董榕波， 李 蕾， 李 海， 馬松云， 劉銀水

(1.華中科技大學(xué) 機(jī)械科學(xué)與工程學(xué)院， 湖北 武漢 430074；2.中國(guó)人民解放軍95829醫(yī)院， 湖北 武漢 430022；3.華中科技大學(xué) 同濟(jì)醫(yī)學(xué)院附屬協(xié)和醫(yī)院兒科， 湖北 武漢 430030；4.亞琛工業(yè)大學(xué) 機(jī)械工程學(xué)院通用力學(xué)研究所， 德國(guó) 亞琛 52062)

引言

吸入給藥是利用肺泡具有較大表面積、毛細(xì)血管豐富、物質(zhì)交換短和藥物吸收快等優(yōu)點(diǎn)，將藥物粉末/液滴分散，經(jīng)由口腔進(jìn)入呼吸道，最終進(jìn)入血液循環(huán)以達(dá)到預(yù)防和治療呼吸系統(tǒng)疾病的一種高效治療手段[1]。目前吸入給藥存在藥物輸送效率較低的突出問題，會(huì)急劇增加患者的醫(yī)療成本，若過量用藥則會(huì)損害患者的身體健康。藥物顆粒在呼吸道內(nèi)的沉積主要受慣性沉積、湍流擴(kuò)散沉積以及顆粒布朗運(yùn)動(dòng)的影響[2]，其中湍流擴(kuò)散沉積起主導(dǎo)作用。因此，精準(zhǔn)的求解人體呼吸道的湍動(dòng)能場(chǎng)，是研究提升吸入給藥效率的根本前提。

人體支氣管樹結(jié)構(gòu)復(fù)雜，呼吸道截面變化巨大，其流動(dòng)雷諾數(shù)存在較大的變化范圍(Re=500～8000)[3]。因此，需借助可靠的湍流模型預(yù)測(cè)包含層流到湍流過渡轉(zhuǎn)變、剪切流和回旋流等復(fù)雜流態(tài)。目前，何種湍流模型能夠高效精準(zhǔn)求解呼吸道流場(chǎng)依舊是一個(gè)開放性問題[4]。使用直接數(shù)值模擬(Direct Numerical Simulation, DNS)的方法求解瞬時(shí)Navier-Stokes(N-S)方程，以此對(duì)湍流進(jìn)行完全解析是最為精確的。但其存在仿真耗時(shí)長(zhǎng)，所需計(jì)算資源龐大等問題，難以高效的求解呼吸道流場(chǎng)[5]。若僅考慮對(duì)湍動(dòng)能場(chǎng)起主導(dǎo)作用的大尺度的渦，而對(duì)小尺度的渦采用數(shù)學(xué)模型加以刻畫，即大渦模擬仿真(Large Eddy Simulation, LES)，則可有效降低DNS仿真中的網(wǎng)格數(shù)量。JAYARAJU等人[6]和LAMBERT等人[7]通過使用LES模型對(duì)簡(jiǎn)化的上呼吸道模型仿真后發(fā)現(xiàn)，LES所計(jì)算的速度和湍動(dòng)能場(chǎng)與實(shí)驗(yàn)數(shù)據(jù)保持較高的一致性。LIN等人[8]分別使用DNS和LES模擬分析了同一呼吸道模型的內(nèi)部流場(chǎng)特性，結(jié)果表明LES能夠較好的模擬流場(chǎng)中的回旋流，且與DNS得到的仿真結(jié)果高度吻合。但是，LES本質(zhì)上還是一種DNS方法，其流場(chǎng)仍需大量網(wǎng)格來表征[9]，計(jì)算耗時(shí)較長(zhǎng)。然而，從臨床角度出發(fā)，醫(yī)生希望能在短時(shí)間內(nèi)得到呼吸道流場(chǎng)的仿真模擬結(jié)果，以此對(duì)呼吸系統(tǒng)疾病做出精準(zhǔn)的診斷。

雷諾平均N-S方程(Reynolds-averaged Navier-Stokes， RANS)則不求解瞬時(shí)的N-S方程，轉(zhuǎn)而計(jì)算時(shí)間平均湍流場(chǎng)。如此，流場(chǎng)表征所需的網(wǎng)格數(shù)量可顯著降低，進(jìn)而有效提升了流場(chǎng)的計(jì)算效率[10]。常用的RANS模型包括渦黏模型(兩方程的k-ε模型和k-ωSST模型，以及四方程的Transition SST模型等)和雷諾應(yīng)力模型(七方程的RSM-SSG模型)。ISLAM等[11]采用不同渦粘模型模擬呼吸道內(nèi)部流場(chǎng)后發(fā)現(xiàn)，k-ωSST模型計(jì)算上呼吸道內(nèi)的湍流時(shí)要優(yōu)于標(biāo)準(zhǔn)的k-ε模型。XU等[12]分別考慮了吸氣和呼氣兩個(gè)過程，采用不同的RANS模型得到的結(jié)果與PIV實(shí)驗(yàn)結(jié)果進(jìn)行比較后發(fā)現(xiàn)Transition SST模型更適用于吸入過程的湍流場(chǎng)仿真，而k-ωSST模型則更適用于呼出過程的仿真。ELCNER[13]通過對(duì)簡(jiǎn)化的呼吸道模型進(jìn)行仿真模擬，分析比較了k-ε、k-ω和LES三種湍流模型。結(jié)果表明LES的結(jié)果與實(shí)驗(yàn)結(jié)果的吻合度最高，其余湍流模型計(jì)算得到的喉部下游流動(dòng)區(qū)域的流場(chǎng)與實(shí)驗(yàn)結(jié)果存在較大區(qū)別。雖然k-ε和k-ω湍流模型均能得到流場(chǎng)的全局特征，但難以較好的求解流速較低時(shí)的低雷諾數(shù)區(qū)域和過渡流區(qū)域。

結(jié)合Spezial-Sarka-Catski(SSG)模型[14]的雷諾應(yīng)力模型(RSM)摒棄了各項(xiàng)同性渦粘假設(shè)，通過求解雷諾應(yīng)力傳輸方程和耗散率方程，以封閉雷諾平均N-S方程。RSM模型比兩方程模型更嚴(yán)格地計(jì)算了流線曲率、渦流、回旋流的變化對(duì)流場(chǎng)的影響，因此它在計(jì)算復(fù)雜湍流流動(dòng)時(shí)具有更高優(yōu)勢(shì)。SOMMERFELD等人[15]通過分析對(duì)比RSM模型和RANS模型后發(fā)現(xiàn)，使用k-ε模型計(jì)算得到的呼吸道流場(chǎng)難以精準(zhǔn)的表征湍流邊界層，使用k-ωSST計(jì)算的流場(chǎng)能夠很好的捕捉到粘性底層與湍流層之間的混合流動(dòng)區(qū)域的細(xì)節(jié)，各向異性的雷諾應(yīng)力(RSM-SSG)湍流模型則能更好的計(jì)算近壁面的湍流場(chǎng)。

此外，由于真實(shí)人體支氣管樹結(jié)構(gòu)復(fù)雜，過往相關(guān)研究均采用混合網(wǎng)格的方法劃分呼吸道內(nèi)部流場(chǎng)，需要消耗巨大的計(jì)算時(shí)間和資源。然而，從臨床角度出發(fā)，醫(yī)生希望在48 h內(nèi)給出基于CFD仿真的吸入給藥效率分析結(jié)果，先前的研究并不能滿足臨床需求。鑒于上述矛盾，若能采用全結(jié)構(gòu)化網(wǎng)格劃分其流場(chǎng)，則能極大的降低網(wǎng)格的數(shù)量(相同尺寸的結(jié)構(gòu)化網(wǎng)格數(shù)量?jī)H為非結(jié)構(gòu)化網(wǎng)格數(shù)量的1/3)，并結(jié)合更適合真實(shí)人體呼吸道內(nèi)湍流仿真模型，進(jìn)而可高效的提升仿真模擬計(jì)算的精度和速度。

綜上所述，本研究提出使用全結(jié)構(gòu)化網(wǎng)格劃分真實(shí)人體支氣管樹模型，并分別采用k-ωSST模型、Transition SST模型、雷諾應(yīng)力模型(RSM-SSG)和大渦模擬(LES)計(jì)算呼吸道內(nèi)部流場(chǎng)，并與實(shí)驗(yàn)測(cè)量結(jié)果相對(duì)比，在確?？煽坑?jì)算精度的前提下，探尋最適用于真實(shí)人體呼吸道流場(chǎng)的快速且精確的計(jì)算方法，為呼吸道計(jì)算流體力學(xué)仿真模擬的臨床應(yīng)用提供參考和依據(jù)。

1 湍流模型

1.1 Shear-Stress Transport (SST) k-ω模型

SST模型是The Baseline (BSL)k-ω模型的改進(jìn)，并在湍流黏度的定義中考慮了湍流剪切應(yīng)力。這使得該算法在近壁區(qū)繞流和旋流方面比標(biāo)準(zhǔn)k-ω和BSLk-ω模型更精確。其湍動(dòng)能k及耗散率ω的方程如下[16]：

(1)

(2)

式中，ui——i方向的平均速度

uj——j方向的平均速度

Γk——k的有效擴(kuò)散系數(shù)

Γω——ω的有效擴(kuò)散系數(shù)

Yk——k因湍流引起的耗散

Yω——ω因湍流引起的耗散

Gk——k引起的湍動(dòng)能

Gω——ω引起的湍動(dòng)能

Dω—— 交叉擴(kuò)散率

有效擴(kuò)散系數(shù)公式為：

(3)

(4)

式中，σk和σω——k和ω的紊流普朗特?cái)?shù)

μt—— 修正后的黏度系數(shù)

(5)

式中，α*—— 抑制湍流黏度的低雷諾修正系數(shù)

S—— 剪切力張量的常數(shù)項(xiàng)

F2—— 湍流普朗特?cái)?shù)的融合項(xiàng)

α1—— 常數(shù)

1.2 Transition SST模型

(6)

轉(zhuǎn)捩源項(xiàng)定義如下：

Pγ1=Flengthca1ρS[γFonset]cγ3

(7)

Eγ1=ce1Pγ1γ

(8)

Pγ2=ca2ρΩγFturb

(9)

Eγ2=ce2Pγ2γ

(10)

式中，S—— 應(yīng)變率張量的模

Flength—— 轉(zhuǎn)捩區(qū)長(zhǎng)度

Ω —— 漩渦強(qiáng)度

Fonset—— 渦量雷諾數(shù)ReV的函數(shù)

Fturb—— 黏性系數(shù)系數(shù)比的函數(shù)

(11)

(12)

T —— 時(shí)間尺度

Fθt—— 由源項(xiàng)控制參數(shù)的混合參數(shù)

1.3 RSM模型

本研究采用的RSM模型為Spezial-Sarka-Catski的SSG模型，該模型已被證明在平面應(yīng)變、旋轉(zhuǎn)平面剪切和軸對(duì)稱膨脹/收縮等一系列基本剪切流中具有優(yōu)越的性能。典型的線性雷諾應(yīng)力模型如下[20]：

(13)

Dij—— 擴(kuò)散項(xiàng)

Pij—— 剪切應(yīng)力產(chǎn)生項(xiàng)

Φij—— 壓力應(yīng)變項(xiàng)

εij—— 和黏性耗散項(xiàng)

1.4 LES模型

對(duì)于LES模型，針對(duì)不可壓縮流動(dòng)，經(jīng)過濾波函數(shù)處理后的N-S方程和連續(xù)方程為[21]：

(14)

(15)

表1 四種湍流模型的優(yōu)劣

2 仿真及實(shí)驗(yàn)設(shè)置

2.1 真實(shí)人體支氣管樹模型

本研究所研究的呼吸道模型是根據(jù)一名成年男子的高分辨率計(jì)算機(jī)斷層掃描(CT圖片)獲得的呼吸道醫(yī)學(xué)影像，通過醫(yī)學(xué)影像三維重構(gòu)軟件，生成的立體光刻(STL)三維幾何模型[26]，如圖1a所示。三維重構(gòu)后的支氣管模型包括進(jìn)氣管、口腔、氣管和支氣管(7級(jí))。氣管處的水力直徑約16～18 mm，第七級(jí)支氣管處的水力直徑約為2～3 mm，氣管橫截面面積變化之大和支氣管樹結(jié)構(gòu)之復(fù)雜是呼吸道內(nèi)部產(chǎn)生不同流動(dòng)狀態(tài)的根本原因。由于實(shí)驗(yàn)的需要，根據(jù)人體肺部五葉十八段的生理結(jié)構(gòu)，末端支氣管被劃分為10個(gè)區(qū)域并與相應(yīng)的出口連接。

圖1 采用的真實(shí)人體呼吸道模型[26]

2.2 網(wǎng)格劃分

網(wǎng)格劃分是呼吸道仿真的重要環(huán)節(jié)之一，若網(wǎng)格計(jì)算單元質(zhì)量低(網(wǎng)格非正交性和歪斜度偏高)， 則會(huì)引發(fā)數(shù)值計(jì)算發(fā)散，進(jìn)而影響流場(chǎng)的求解精度[27]。由于口腔的模型較為復(fù)雜，此處采用非結(jié)構(gòu)化網(wǎng)格劃分。其余部位則充分利用ICEM CFD在結(jié)構(gòu)化網(wǎng)格劃分方面的優(yōu)勢(shì)，將支氣管樹及10個(gè)收集器全部采用結(jié)構(gòu)化網(wǎng)格劃分。KOULLAPIS等[9]通過對(duì)比不同網(wǎng)格數(shù)量(700萬(wàn)～5000萬(wàn))對(duì)結(jié)果的影響，發(fā)現(xiàn)在保持較好精度的條件下，采用700萬(wàn)非結(jié)構(gòu)化網(wǎng)格就能夠較好的解析呼吸道內(nèi)的流動(dòng)。對(duì)于同一個(gè)模型，相比于非結(jié)構(gòu)化網(wǎng)格，使用結(jié)構(gòu)化網(wǎng)格能夠大幅降低網(wǎng)格數(shù)量，且能更好的控制網(wǎng)格生成質(zhì)量，使計(jì)算結(jié)果更易收斂。近壁面流動(dòng)是影響顆粒沉積的主要原因[28]。為避免近壁面邊界層網(wǎng)格厚度過低從而影響壁面壓力梯度的求解精度，在網(wǎng)格劃分時(shí)近壁面貼合的網(wǎng)格層數(shù)保持在15層以上(如圖2所示)。該支氣管樹的網(wǎng)格總數(shù)約為1200萬(wàn)。

圖2 采用結(jié)構(gòu)化網(wǎng)格劃分支氣管樹模型：結(jié)構(gòu)化網(wǎng)格總數(shù)為900萬(wàn)，呼吸道壁面邊界層網(wǎng)格層數(shù)大于15層

2.3 實(shí)驗(yàn)設(shè)置與邊界條件

人體在不同的運(yùn)動(dòng)狀態(tài)或健康狀態(tài)下，呼吸速率和呼吸量會(huì)有較大的差別。根據(jù)吸入氣體的速度是否恒定(吸入速度是否隨時(shí)間變化，吸入速度不隨時(shí)間變化為穩(wěn)態(tài)呼吸，反之)，可以將呼吸狀態(tài)分為穩(wěn)態(tài)呼吸和非穩(wěn)態(tài)呼吸。本研究主要研究在穩(wěn)定呼吸狀態(tài)下持續(xù)吸入0.5 s的過程，即Q(t) = 60 L/min，(0

表2 吸入流量Q=31.1 L/min時(shí)，10個(gè)出口的流速

表3 仿真邊界條件和數(shù)值方法設(shè)置

圖3 實(shí)驗(yàn)裝置：模型的進(jìn)氣管和口腔完全浸沒于裝有水-甘油混合液體的玻璃水箱中，模型的10個(gè)出口置于水箱下部并分別與10個(gè)活塞隔膜泵連接[29]。

3 計(jì)算結(jié)果與分析

圖4顯示了當(dāng)t=0.5 s時(shí)，采用LES算法得到的進(jìn)氣口到其中一個(gè)支氣管出口在呼吸道冠狀面和矢狀面上的時(shí)間平均流速分布(仿真中每10個(gè)時(shí)間步長(zhǎng)求一次平均作為該時(shí)間區(qū)間內(nèi)的時(shí)間平均值)。其中Ⅰ，Ⅱ，Ⅲ和Ⅳ分別表示口腔和氣管的矢狀面、氣管分叉處的冠狀面、左主支氣管的冠狀面和右主支氣管的冠狀面的具體位置。線段A，B，C，D，G，H和J分別表示支氣管樹內(nèi)矢狀面或冠狀面上的截線。由圖4可知，隨著支氣管橫截面面積逐漸減小，氣流速度逐漸增大。每個(gè)支氣管分叉處的速度變化最為劇烈，會(huì)出現(xiàn)較為明顯的剪切層和回旋流區(qū)。為了更詳細(xì)的對(duì)比不同湍流模型計(jì)算結(jié)果的差異，取t=0.5 s這一時(shí)刻點(diǎn)的相關(guān)云圖和數(shù)據(jù)(不再贅述)，對(duì)支氣管樹的不同區(qū)域進(jìn)行詳細(xì)分析。

圖4 通過LES模型計(jì)算得到的呼吸道內(nèi)從進(jìn)口到其中一個(gè)出口的冠狀面和矢狀面時(shí)間平均流速分布；A,B,C,D,G,H和J分別表示支氣管樹內(nèi)冠狀面或矢狀面與相應(yīng)位置處橫斷面的交線；Ⅰ,Ⅱ,Ⅲ和Ⅳ分別展示了口腔和氣管的矢狀面、氣管分叉處的冠狀面、左主支氣管的冠狀面和右主支氣管的冠狀面的具體位置。

圖5a顯示了不同RANS湍流模型和LES模型在口腔和氣管矢狀面上的計(jì)算得到的結(jié)果。該處支氣管橫截面面積較大，氣體流動(dòng)較慢，雷諾數(shù)常在3000～5000。從圖5a中可看出，氣體經(jīng)過口腔(交線B處)和咽部(交線C處)時(shí)速度較低，在口腔上壁出現(xiàn)了流動(dòng)分離現(xiàn)象。由于咽喉收縮，氣體流經(jīng)此處時(shí)速度突然增大，并伴隨著剪切流和回旋流。隨著氣管結(jié)構(gòu)產(chǎn)生彎曲，呼吸道內(nèi)部高速流從氣管后壁轉(zhuǎn)移到氣管前壁，并且在前壁產(chǎn)生了一個(gè)較小的速度分離區(qū)域，該現(xiàn)象在LES模型的仿真結(jié)果中較為明顯。氣管中段區(qū)域，LES模型仿真計(jì)算的速度從氣管前壁到氣管后壁過渡的較為均勻。k-ωSST模型的計(jì)算結(jié)果和RSM-SSG模型的計(jì)算結(jié)果較為一致。在Transition SST模型的仿真中，并沒有在靠近氣管前壁面區(qū)域產(chǎn)生較大的流速，流速最大處主要集中在氣管中部。

對(duì)于上呼吸道內(nèi)部的流場(chǎng)，LES模型能夠捕捉到氣管內(nèi)更詳細(xì)的剪切流和回旋流。不同RANS湍流模型得到的總體速度分布具有較好的一致性，但是在預(yù)測(cè)喉部的剪切層和回旋流時(shí)有較大的差異。k-ωSST模型與RSM-SSG模型在近壁面區(qū)域計(jì)算的結(jié)果與LES模型計(jì)算的結(jié)果吻合性較好。Transition SST的結(jié)果與LES的結(jié)果差異較大。這是由于口腔和氣管中常存在回旋流和自由剪切流，Transition SST轉(zhuǎn)捩流模型計(jì)算這些流動(dòng)時(shí)效果不佳[32]?？谇缓蜌夤艿臅r(shí)間平均湍動(dòng)能如圖5b所示，湍動(dòng)能在入口處較小，隨著流動(dòng)的發(fā)展，在口腔的上部開始升高。最大的湍動(dòng)能水平位于喉部氣管前壁和喉部下游區(qū)域的氣管后壁。k-ωSST模型和RSM-SSG模型明顯低估了該區(qū)域的湍動(dòng)能水平。

圖5 使用不同RANS湍流模型和LES模型在口腔和氣管矢狀面上的計(jì)算結(jié)果(t=0.5 s)

為了進(jìn)一步對(duì)比結(jié)果，取圖4中上呼吸道中的A,B,C和D四條交線處的實(shí)驗(yàn)測(cè)量值、不同RANS湍流模型計(jì)算結(jié)果和LES模型計(jì)算結(jié)果進(jìn)行比較。為了處理數(shù)據(jù)方便和更直觀地對(duì)比各湍流模型計(jì)算結(jié)果與實(shí)驗(yàn)的差異，在每條交線上平均取200個(gè)離散點(diǎn)并編號(hào)，進(jìn)而提取各個(gè)點(diǎn)上的時(shí)間平均速度和時(shí)間平均湍動(dòng)能并進(jìn)行無量綱化，最后將編號(hào)進(jìn)行歸一化處理，均映射到0～1范圍之內(nèi)的區(qū)間上，下同，如圖6所示。圖中D*表示交線上點(diǎn)的歸一化坐標(biāo)。umag和k分別表示仿真得到的交線上的時(shí)間平均速度和時(shí)間平均湍動(dòng)能。uT=2.28 m/s表示水-甘油混合液體在Q=31.1 L/min 時(shí)，進(jìn)氣管橫截面的時(shí)間平均速度。4種湍流模型在交線A處得到的速度分布曲線較為一致，但與實(shí)驗(yàn)有明顯差異。原因是實(shí)驗(yàn)在裝滿水-甘油混合液體的水箱中進(jìn)行的，水箱體積有限、壁面離進(jìn)氣口較近等因素對(duì)實(shí)驗(yàn)結(jié)果影響較大[15]，從而使得入口處速度不對(duì)稱。氣體進(jìn)入口腔后，所有湍流模型在B處均計(jì)算出了剪切流。隨著流動(dòng)進(jìn)一步發(fā)展，剪切層逐漸轉(zhuǎn)移到上咽(交線C處)內(nèi)側(cè)，Transition SST模型的結(jié)果和RSM-SSG模型的結(jié)果與實(shí)驗(yàn)結(jié)果更為吻合。交線D處的呼吸道橫截面突然增大，速度變化較為劇烈，四種模型得到的速度曲線與實(shí)驗(yàn)結(jié)果具有較高的一致性，但是均不能精確的捕捉到剪切層的位置，其中Transition SST的表現(xiàn)最為不佳。

圖6 使用不同RANS湍流模型和LES模型在A,B,C和D四條交線處計(jì)算得到的無量綱化時(shí)間平均速度曲線(左圖)和無量綱化時(shí)間平均湍動(dòng)能曲線(右圖)(t=0.5 s)；a),b),c)和d)分別表示A,B,C和D四條交線；D*表示交線上點(diǎn)的歸一化坐標(biāo)；umag和k分別表示仿真得到的交線上的時(shí)間平均速度和時(shí)間平均湍動(dòng)能；uT=2.28 m/s表示水-甘油混合液體在Q=31.1 L/min時(shí)，進(jìn)氣管橫截面的時(shí)間平均速度

與時(shí)均平均速度場(chǎng)相比，各個(gè)模型計(jì)算的湍動(dòng)能明顯不同。在進(jìn)氣管處，LES模型低估了湍動(dòng)能水平，而k-ωSST、 Transition SST和RSM-SSG模型在進(jìn)氣管中心區(qū)域得到地湍動(dòng)能與實(shí)驗(yàn)結(jié)果較吻合。其中，RSM-SSG模型在近壁面處捕捉到了較大的湍動(dòng)能，這是由于仿真設(shè)置的進(jìn)口邊界條件為大氣壓力，與進(jìn)口周圍靜態(tài)環(huán)境相比進(jìn)器管口產(chǎn)生較大的速度變化，從而表現(xiàn)出更高的湍動(dòng)能。在B,C和D三條交線處，LES和Transition SST模型得到的湍流強(qiáng)度較為接近，能夠捕捉到湍動(dòng)能的整體分布和大小。k-ωSST模型和RSM-SSG模型得到的湍流強(qiáng)度均低于實(shí)驗(yàn)。RSM-SSG模型計(jì)算的湍動(dòng)能沿著氣管衰減最快，但其仍能夠在近壁面捕捉到與實(shí)驗(yàn)更接近的湍流強(qiáng)度。

隨著氣流進(jìn)入左右支氣管，氣管分叉處和左右主支氣管(結(jié)構(gòu)特征由人體生物學(xué)定義)的時(shí)間平均速度分布云圖如圖7所示。在Lizal[26]的體外測(cè)量實(shí)驗(yàn)中，左肺和右肺的通氣量之比為3∶7，這種不相等的通氣量導(dǎo)致流速在氣管分叉處的左右區(qū)域差異巨大。在氣管的右壁處可以看到一個(gè)高速流動(dòng)區(qū)域，此處的雷諾數(shù)能達(dá)到6000以上，而左側(cè)分叉處出現(xiàn)了一個(gè)較大的回旋流區(qū)域。隨著左主支氣管下游變窄和變短，流速逐漸增大，該區(qū)域產(chǎn)生了較大的剪切流區(qū)域。右主支氣管與下一級(jí)支氣管分支角度較大，這種幾何形狀的不對(duì)稱導(dǎo)致右主支氣管內(nèi)出現(xiàn)射流現(xiàn)象。

圖7 使用不同RANS湍流模型和LES模型計(jì)算得到的時(shí)間平均速度場(chǎng)(t=0.5 s)

氣管分叉處和左右主支氣管的湍動(dòng)能云圖如8所示，盡管湍流強(qiáng)度在氣管中出現(xiàn)了衰減，但是在剪切流和回旋流區(qū)域之間出現(xiàn)了較高的湍流水平，最大的湍流強(qiáng)度出現(xiàn)在右主支氣管。LES和Transition SST模型均能夠較好計(jì)算出左右主支氣管處的湍動(dòng)能大小，k-ωSST模型和RSM-SSG模型依舊低估了此處的湍動(dòng)能水平。如上所述，RSM-SSG雖然不能準(zhǔn)確預(yù)測(cè)呼吸道中心區(qū)域的湍流水平，但是在近壁面處能夠捕捉到與實(shí)驗(yàn)較為吻合的湍流強(qiáng)度。

圖8 使用不同RANS湍流模型和LES模型計(jì)算的時(shí)間平均湍動(dòng)能場(chǎng)(t=0.5 s)

圖9顯示了G,H和J三條交線處的無量綱化的時(shí)間平均速度和時(shí)間平均湍動(dòng)能曲線，在G處可以發(fā)現(xiàn)，LES模型計(jì)算的結(jié)果和實(shí)驗(yàn)測(cè)量結(jié)果具有更好的一致性，Transition SST模型得到的結(jié)果與LES計(jì)算的結(jié)果最為接近。k-ωSST和RSM-SSG模型在此處的仿真結(jié)果均略高于實(shí)驗(yàn)結(jié)果。隨著氣流進(jìn)入左右主支氣管，氣管內(nèi)開始出現(xiàn)了較為強(qiáng)烈的剪切流動(dòng)，LES模型的計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果最為接近。湍動(dòng)能方面，LES模型與Transition SST模型在G處計(jì)算得到的湍動(dòng)能分布與實(shí)驗(yàn)較為一致，二者均能捕捉到此處湍動(dòng)能的最大值和最小值。k-ωSST模型和RSM-SSG模型依舊低估了呼吸道中心處的湍動(dòng)能水平，但是RSM-SSG仍然能很好的捕捉到近壁面的湍流強(qiáng)度。綜合比較三處交線的湍動(dòng)能曲線可以發(fā)現(xiàn)，LES計(jì)算低級(jí)支氣管處的湍動(dòng)能具有很大的優(yōu)勢(shì)，其計(jì)算的湍動(dòng)能與實(shí)驗(yàn)測(cè)量值保持較高的一致性。k-ωSST模型在H處計(jì)算的湍動(dòng)能與LES模型的結(jié)果最接近，Transition SST模型和RSM-SSG模型計(jì)算的湍動(dòng)能與實(shí)驗(yàn)結(jié)果的吻合度較低。

圖9 使用不同RANS湍流模型和LES模型在交線G,H和J三處計(jì)算的無量綱化時(shí)間平均速度曲線(左圖)和無量綱化時(shí)間平均湍動(dòng)能曲線(右圖)(t=0.5 s)；a),b)和c)分別表示G,H和J三條截線；D*表示橫坐標(biāo)歸一化區(qū)間；umag和k分別表示仿真得到的截線上的時(shí)間平均速度和時(shí)間平均湍動(dòng)能；uT=2.28 m/s表示水-甘油混合液體在Q=31.1 L/min時(shí)，進(jìn)氣管橫截面的時(shí)間平均速度。

4 結(jié)論

本研究使用真實(shí)人體呼吸道模型，采用全結(jié)構(gòu)化網(wǎng)格劃分支氣管樹，并分別使用k-ωSST模型、Transition SST、雷諾應(yīng)力模型(RSM-SSG)和大渦模擬(LES)計(jì)算呼吸道流場(chǎng)，通過與實(shí)驗(yàn)測(cè)量結(jié)果比較發(fā)現(xiàn)：

(1) 采用結(jié)構(gòu)化網(wǎng)格劃分支氣管樹能夠極大的降低網(wǎng)格數(shù)量，進(jìn)而可有效提升計(jì)算效率。

(2) 在上呼吸道區(qū)域，由于上呼吸道的橫截面面積較大，氣流在此處的流動(dòng)處于發(fā)展階段，雷諾數(shù)相對(duì)較低。使用不同湍流模型在該區(qū)域計(jì)算得到的速度場(chǎng)中，Transition SST模型在流動(dòng)混合區(qū)域的計(jì)算與實(shí)驗(yàn)結(jié)果存在差異，k-ωSST模型、雷諾應(yīng)力模型(RSM-SSG)和大渦模擬(LES)的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的吻合度較高。在該區(qū)域湍動(dòng)能場(chǎng)計(jì)算結(jié)果的比較中發(fā)現(xiàn)，LES不僅能較為準(zhǔn)確地預(yù)測(cè)湍動(dòng)能的分布狀態(tài)，而且能捕捉到喉部下游的湍流細(xì)節(jié)；Transition SST模型雖在預(yù)測(cè)近壁面湍動(dòng)能時(shí)表現(xiàn)欠佳，但在預(yù)測(cè)呼吸道中心區(qū)域的湍動(dòng)能水平時(shí)與LES最接近；k-ωSST和RSM-SSG湍流模型的計(jì)算結(jié)果對(duì)湍動(dòng)能的預(yù)測(cè)均低于實(shí)驗(yàn)測(cè)量結(jié)果，但RSM-SSG模型在近壁面處的湍動(dòng)能計(jì)算結(jié)果與實(shí)驗(yàn)保持一致。由于藥物顆粒的湍流擴(kuò)散沉積主要與近壁面湍流場(chǎng)計(jì)算精度相關(guān)，因此RSM-SSG的計(jì)算結(jié)果是可接受的。

(3) 在支氣管區(qū)域，由于該區(qū)域的速度變化劇烈，常伴隨較強(qiáng)的剪切流和回旋流，雷諾數(shù)相對(duì)較高。LES和RSM-SSG在該區(qū)域的流動(dòng)狀態(tài)和湍動(dòng)能計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的差異較小。k-ωSST和Transition SST在此區(qū)域的速度場(chǎng)計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果較為接近，但近壁湍動(dòng)能場(chǎng)計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果差異較大。

綜上所述，在計(jì)算真實(shí)人體呼吸道內(nèi)部流場(chǎng)時(shí)，LES的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果最為接近。但當(dāng)受限于計(jì)算效率時(shí)，結(jié)合結(jié)構(gòu)化網(wǎng)格的巨大優(yōu)勢(shì)，使用RSM-SSG模型代替LES模型，可以有效節(jié)約計(jì)算時(shí)間，能夠最大程度的減小計(jì)算精度與LES模型的差異對(duì)顆粒湍流擴(kuò)散沉積的影響，從而為計(jì)算流體力學(xué)仿真模擬大規(guī)模應(yīng)用于呼吸道疾病的臨床應(yīng)用提供參考和奠定基礎(chǔ)。

致謝

作者十分感謝科技部國(guó)際合作司中國(guó)與斯洛文尼亞科技人員交流項(xiàng)目對(duì)該研究的支持,以及COST (European Cooperation in Science and Technology; Simulation and pharmaceutical technologies for advanced patient-tailored inhaled medicines; www.cost.eu)對(duì)本研究提供的幫助。