可吸入顆粒物在不同阻塞性呼吸道內運動與沉積特性

莊加瑋,刁永發,楚明浩,沈恒根 (東華大學環境科學與工程學院,上海 201620)

固體顆粒物是工業建筑中常見的污染物,可吸入顆粒物(PM10)進入呼吸道后會積聚在肺部,是工業塵肺病的主要誘因[1-2].慢性阻塞性肺病是塵肺患者常見且嚴重的合并癥之一,發病率和死亡率都很高[3].流行病學研究也表明,長期暴露于高濃度顆粒污染物中會對人體呼吸系統造成嚴重損害[4-5].因此,掌握顆粒在慢性阻塞性肺病工人呼吸道內運動和沉積特性將有助于定量評估其暴露風險,對于保障工人健康及后期的治療有實際應用價值.

顆粒污染物在人體呼吸道內傳輸和沉積涉及流體力學、環境學、生物醫學等多學科交叉,早期有學者通過實驗建立理想的局部呼吸道模型來研究顆粒的沉積率變化規律[6-8].隨著 3D打印技術發展,氣溶膠顆粒在真實人體呼吸道內沉積的研究也取得了進展[9-10].相關研究在Weibel-23級理想呼吸道模型基礎上,建立了G3～G5的數值計算幾何邊界,通過實驗證實了數值模型的準確性,為后期呼吸道兩相流數值計算奠定了基礎[11-12].目前CFD數值模擬已成為研究呼吸道內顆粒物運輸和沉積特性的主流方法,研究內容包括:顆粒物的沉積機制[13-14]、不同特性顆粒物的沉積規律[15-17]、人體呼吸機制對顆粒物沉積的影響[18-19]、呼吸道內霧化給藥[20-21]及各類參數(粒徑、呼吸強度、呼吸道溫度、濕度等)對顆粒物運輸的影響[22-24]等.為提高數值計算的準確性,Zhang等[25]發現呼吸的瞬時特性對顆粒運輸的影響不可忽視.近年來,有學者將 CFD-DEM 方法應用于呼吸道顆粒物運動研究,結果表明該方法可提高預測精度,有望成為解決濃稠型氣溶膠肺部動力學問題的新方法[26-27].此外,呼吸道病變會改變原有幾何結構,改變氣流分布,進而影響顆粒在呼吸道內運動和沉積形式[28-30].

現有研究主要聚焦顆粒在健康呼吸道內的沉積特性,關于顆粒在病變呼吸道內沉降特性的研究相對較少,而患者的實際病情不盡相同,關于呼吸道病變位置和阻塞率對顆粒沉積影響的研究更是鮮有報道.為此,本文擬構建兩類(G3～G6和 G11～G14)阻塞型呼吸道模型,通過 CFD-DPM 數值方法考察局部受限下呼吸道內兩相流運動和沉積特性,分析阻塞率(α)、阻塞位置、勞動強度等因素對呼吸道內流場、顆粒沉積形式和沉積率的影響,以期為職業病的預防和靶向給藥治療提供參考.

1 數學物理模型

1.1 氣固兩相流運動控制方程

人體吸氣過程氣體運動是典型的管內流動,考慮到呼吸道局部阻塞影響,流動將出現湍流和層流,故采用標準k-ε湍流模型進行數值計算,對應的連續性方程和動量方程為[29]:

一般情況,呼吸過程氣體中攜帶的顆粒濃度較低,顆粒相可視為稀疏相,即可忽略顆粒間的碰撞作用[13-17].假設顆粒為單分散球形粒子,僅考慮重力與曳力對顆粒運動的影響,依據牛頓第二定律,單一顆粒運動過程受的力方程為[30]:

式中:md為顆粒的質量,kg;dp為顆粒的直徑,μm;ρp為顆粒密度,kg/m3;p為顆粒的速度矢量,m/s;CD為曳力系數;為重力加速度,m/s2,方向沿Z軸正方向.

1.2 物理模型

慢性阻塞性肺病患者由于長期慢性炎癥刺激和人體自我修復,使得呼吸道邊壁增厚,形成了不可逆的氣道阻塞[3].基于 Weibel呼吸道模型,人體呼吸道由氣管(G0)到肺泡(G23)共 24代組.其中,G0～G16和G17～G23分別為氣體的傳輸區和交換區[31].本文考察顆粒在傳輸區的沉積特性,因此選取了 2個不同呼吸道區域,即G3～G6和G11～G14,均為四級三分叉共面對稱支氣管,分叉角為 70°.假設氣道阻塞都發生在第二級支氣管(即G4及G12),阻塞率(α)由第二級支氣管中部直徑與初始直徑比值表示.圖 1給出了呼吸道結構平面示意[31,35],模型通過Solidworks三維繪圖軟件構建,呼吸道具體尺寸如表1[32].

表1 呼吸道幾何參數(mm)Table 1 Geometric parameters of the respiratory tract mm(mm)

圖1 受COPD影響的呼吸道模型平面示意Fig.1 Schematic view of the airway models affected by COPD

1.3 數值方法與邊界條件

采用有限體積法求解控制方程,通過標準格式對壓力項進行離散,動量方程利用二階迎風格式進行離散,用SIMPLE算法進行壓力速度耦合求解,連續性方程和動量方程的收斂殘差均設置為 10-6.流場到達穩態后,利用離散相拉格朗日法求解顆粒運動,氣體和顆粒間為單向耦合,即僅考慮流場對顆粒運動的作用.

為模擬呼吸道內氣體流動,入口處設置為速度進口邊界條件(velocity-inlet),速度大小沿入口圓形截面成完全拋物線分布,并通過自定義函數(UDF)進行定義.出口均設置為壓力出口(pressure-outlet),四周管壁都采用無滑移邊界(no-slip wall).本研究中,人體的呼吸流量(Qin)分別為 15,30,60L/min,分別代表休息、輕度和中度勞動.入口處平均速度大小()可通過下式計算得到[30]:

式中: n為呼吸道級數;Dn為呼吸道直徑,mm.

顆粒與氣體同時進入呼吸道,呼吸道壁面的黏液層可有效黏附與其碰撞的粒子,故認為顆粒與壁面發生接觸即沉積.因此,所有的管壁都設為捕集(trap)邊界,入口與出口均設為逃逸(escape)邊界.為使顆粒注入的條件與實際相似,入口處顆粒采用隨機拋物線分布,即某一位置處顆粒的速度與氣流速度一致,顆粒的頻率分布與其速度分布一致[29].為排除注入顆粒數量對研究結果的影響,分別對顆粒總數NT= 1000,5000,10000,20000;Qin= 30L/min,dp=5μm的4個工況進行了數值計算,發現當顆粒數NT=10000時,顆粒總沉積率已基本到達穩定,且同 NT=20000時差異甚小,誤差不超過3%,故本文在數值計算時注入呼吸道的顆粒數量為10000.粒子的隨機拋物線分布可通過下式表達[33]:

式中: C(r)和 C分別為入口截面局部區域及整體平均顆粒數濃度,個/mm;R為入口處呼吸道半徑,mm;Int[ ]為取整函數;Np為ra≤r≤rb圓環內顆粒數量.粒子注入位置的隨機生成通過一個編譯的內置 Matlab代碼實現.圖2給出了G3～G6呼吸道入口處顆粒數量沿徑向的頻率分布.

圖2 顆粒數量沿徑向的頻率分布Fig.2 Frequency distribution of particle number along radial direction

1.4 主要物理量計算方法

呼吸道入口處氣流的平均雷諾數(Rein),定義為:

斯托克斯數(St)是一個無量綱數,表征顆粒在流體中的跟隨性,斯托克斯數越大,顆粒的慣性越大,跟隨性越差.其定義為:

顆粒沉降過程中空氣流過呼吸道總的壓力損失為:

式中:pin為計算域進口處壓力,Pa;out為計算域出口處平均壓力,Pa.

顆粒物的沉積狀況通過沉積效率(η)定量描述,定義為:

式中: Npd為沉積在所選區域的顆粒數量;Npt為進入到整個呼吸道區域的顆粒數總量.

為研究不同機制對顆粒沉積的作用,將呼吸道內總的沉積效率(ηt)分為兩部分,具體可表達為[30]:

式中: ηi和ηg分別為因慣性碰撞和重力作用發生的沉積效率.

2 結果與討論

2.1 網格無關性與數值方法驗證

在數值計算中,采用非結構性四面體網格劃分計算域,并對邊壁區域進行局部加密處理.為驗證網格獨立性,以健康的G3～G6呼吸道為例,構建了總數分別為669087,970077,1501412和1825582 4類網格,比較了不同網格數量下呼吸道內氣流總的壓力損失(Δpt)和顆粒沉積效率(ηt).

如圖3所示示,Δpt隨Rein的增大而增加,不同雷諾數下,網格數 669087,1501412和 1825582間 Δpt的最大差異分別為3.54%和0.68%.此外,不同網格數量下,ηt都隨 St的增大而提高,且沉積效率的最大差異都小于 3%.為平衡計算的精度與效率,最終采用網格數為 1501412的網格進行后續模擬計算.經網格獨立性驗證后,本文所有計算區域網格總數在110～180萬之間.

圖3 不同網格數量下呼吸道內總壓力損失和沉積效率Fig.3 Total pressure drop and deposition efficiency in respiratory tractfor different grid numbers

Kim 等[34]通過實驗方法研究了不同粒徑和氣體流量下顆粒在G3～G5呼吸道內的沉積特性.為驗證數值模型的可靠性,構建了與該實驗相同的呼吸道模型,圖4給出了本文模擬結果與Kim等[34]實驗測量結果及Feng等[27]、Chen等[29]數值計算結果的比較.

由圖 4可以看出,不論第一級分叉或第二級分叉,本文模擬結果與前人的實驗結果及數值結果均保持相似的變化趨勢,且具備較高的吻合度.其中,第一級分叉處顆粒沉積效率的模擬值要稍高于實驗值,而第二級分叉處顆粒沉積效率的模擬值要稍低于實驗值.盡管如此,本文模擬結果與實驗結果的吻合度要稍優于 Feng等[27]、Chen等[29].綜上表明,本文采用的數值模型可以準確預測顆粒在呼吸道內的運動與沉積規律.

圖4 呼吸道內顆粒物沉積率實驗結果和數值結果的比較Fig.4 Comparisons between simulated and measured results of particle deposition efficiency in respiratory tract

2.2 呼吸道內流場分布

圖5和圖6分別給出了Qin為15和60L/min時,α分別為0、0.2、0.5和0.8情況下G3～G6和G11～G14呼吸道內無量綱速度分布、受 COPD影響的 G4.2和G12.2的壓力、前后徑向截面(A-A’和B-B’)的速度和二次流矢量分布.

圖5 休息狀態下不同阻塞率呼吸道內氣體流動特性Fig.5 Airflow characteristics in respiratory tract with different obstruction rate under rest state

圖5表明,健康呼吸道內的氣體會在B3和B11處均勻地分為兩支,且速度分布在平面上都關于 Z軸對稱.G3 ～ G6內氣體軸向上的速度中心會明顯偏向于呼吸道內側,形成了較高的壓力,在徑向平面(A-A’和 B-B’)產生了二次流;而 G11 ～ G14 軸向速度中心基本未發生偏移,且徑向平面的二次流比較微弱.此外,徑向平面的速度分布與入口處相似,即沿軸向成拋物線分布.造成上述現象的主要原因是,B3和B11分叉處曲率的存在導致氣流向呼吸道內側集中,且氣流流速越大,氣體偏移愈明顯.

隨著α增大,通過病變呼吸道側的空氣流量逐步減少,呼吸道內速度分布不再對稱,徑向平面(A-A’和 B-B’)上的二次流也逐漸減弱.此外,在病變側(G4.2和 G12.2)產生了滯止和回流區(如圖 5中C1-C3所示).這是因為,呼吸道受阻后,阻塞位置上游的壓力升高,流體通過 A-A’向下游運動時需克服壓力,速度減低并在壁面附近減速至零形成滯止區 C1;另一方面,因噴管效應,氣流通過狹窄區域后形成了射流,壓力迅速降低,流體需克服反向壓力向前流動,形成回流區 C2和 C3.值得注意的是,C1-C3區域范圍會隨著α的增大而發展,從速度分布平面上看,當α=0.8,G4.2區C1已經超過管徑的 3/4,而 G12.2的 C1更是達到了 1,且 C2和 C3有相似的情況.此外,G3～G6的滯止和回流區起初均為月牙形,而G11 ～G14的滯止和回流區均為環形,隨著 α的增大,都會不斷向管道中心凸起,使得滯止和回流區不斷擴大.

對比圖5和圖6發現,勞動強度提高后呼吸道內流體的絕對速度雖有增大,但并未改變無量綱速度分布.由于呼吸量的增加,G3～G6徑向平面(A-A’和B-B’)的二次流顯著增強,并形成了關于中平面對稱的兩個渦結構;且射流現象更加明顯,盡管阻塞位置上游的壓力有所增大.上述現象導致,當α相同時,呼吸道內滯止和回流區范圍略有縮小.

從圖5和圖6可以得出,阻塞型呼吸道內氣體流量分配呈明顯的不均勻性,故而病人會出現局部缺氧癥狀,且同α關系密切.為定量分析病變位置、α及勞動強度對流量分配的影響,表2給出了受COPD影響下病變側 G14分支處氣體的平均流量,同時列出了流量的相對誤差(σ),定義為:

圖6 中等勞動強度下不同阻塞率呼吸道內氣體流動特性Fig.6 Airflow characteristics in respiratory tract with different obstruction rate under moderate labor intensity

式中: Qc和 Qh分別為病變和健康呼吸道下游 G14分支處氣體的平均流量,L/min.

表2表明,受COPD影響,呼吸道下游氣體流量會有不同程度的減少,且α越大,氣體流量相對誤差越大,病人局部缺氧愈嚴重.在病變位置與病變程度相同的前提下,呼吸道G14分支處氣體絕對流量會隨勞動強度的提升而增加,但氣流流量的相對誤差基本保持不變.這意味著,COPD 病人從事相同強度體力勞動時呼吸可能會更急促,高強度勞動下也有發生哮喘的可能性.此外,對比不同病變位置下G14分支處氣體流量可以得到,呼吸道病變位置越深,氣體流量相對誤差稍大,表明病人的局部缺氧越重,而缺氧范圍卻越小.

表2 受COPD影響下呼吸道G14分支處平均體積流量對比Table 2 Average flow rate of lower respiratory tract G14 affected by COPD

2.3 顆粒的沉積形式與機制

依據以上分析結果,由于阻塞管的存在,呼吸道內流場分布發生了明顯變化,這在一定程度上可能會改變管道內顆粒物的沉積形式.

由圖7可以看出,對于健康或阻塞程度比較低的G3 ～ G6呼吸道,受慣性碰撞作用,顆粒會在分叉處形成沉積熱點,且大多沉積在熱點附近壁面.比較圖7(a)和圖 7(b)發現,α增大使得更多氣體流入 G4.1支管,顆粒受到的慣性力得到增強,故而有更多的顆粒沉積在 G4.1支管下游.部分顆粒在噴管射流作用下通過阻塞區域時發生沉積,通過病變位置的顆粒在重力和慣性力作用下繼續沉積,但沉積顆粒數量明顯減少.比較圖 7(b)和圖 7(c)可以得到,呼吸道內顆粒的沉積量隨著dp的變大明顯增加,尤其在G4.1支管下游,這是因為大顆粒運動過程中受到的重力和慣性力增強,沉積率變大,分布范圍也會更大.此外,通過圖 7(b)和圖7(d)的對比還可觀察到,Qin提升同樣會增大顆粒的沉積率和沉積范圍.

圖7 不同工況下G3 ～ G6呼吸道內顆粒物沉積形式的可視化Fig.7 Visualization of particle deposition patterns in G3 ～ G6 respiratory tract under different working conditions

由圖 8可以發現,G11～G14呼吸道內顆粒沉積相對分散,且顆粒大多沉積在各支管的內側,這是因為,G11～G14呼吸道內氣體流速比較低,顆粒運動過程主要受重力作用發生沉降.當呼吸道病變加重,即α=0.8時,除少部分顆粒在通過阻塞區域發生沉積,病變呼吸道G12.2下游顆粒沉積數量也會明顯減少;此外,由于呼吸道健康側流量增加,顆粒受到的慣性力增強,顆粒沉積位置會向下游移動.比較圖 8(b)和圖8(c)可以得到,dp增大使得顆粒在呼吸道內局部沉積密度明顯變大,但并未改變沉積形態.相反,通過對比圖8(b)和圖8(d)發現,增大Qin會在一定程度降低顆粒的局部沉積密度,這是因為,呼吸道內氣流速度增大,顆粒受到的慣性作用增強,而重力作用削弱,導致更多的顆粒流向下游.

為定量描述顆粒沉積分布特征和機制,將呼吸道劃分為7個不同的沉積區域(圖1),根據顆粒沉積的坐標信息,統計出各區域顆粒沉積數量Npd和呼吸道內顆粒總的沉積數量 Npz,由式(10)和(11)分別計算η,ηi和ηg.采用與圖7和圖8相同的計算工況,呼吸道內顆粒的沉積分布律如圖9所示.

圖8 不同工況下G11 ～ G14呼吸道內顆粒物沉積形式的可視化Fig.8 Visualization of particle deposition patterns in G11 ～ G14respiratory tract under different working conditions

圖9 不同工況下G3～G6和G11～G14呼吸道表面顆粒沉積分布Fig.9 Particle distribution on the G3～G6and G11～ G14respiratory tractunder different working conditions

圖 9的結果定量反映出,阻塞型呼吸道內粒子沉積分布的不對稱性,且呼吸道健康側粒子沉積數量要明顯高于病變側.對于G3 ～ G6呼吸道,當α = 0.2,Qin=30L/min,dp= 2.5μm時,顆粒沉積分布的不對稱性較小,主要表現為區域4和5的差異;當α增大到0.8時,粒子沉積分布的不對稱性有所提高,表現為區域2和3及區域4和5的不同,沉積率差異分別為1.1% 和1.34%;當 dp增大到 5.0μm,α = 0.8,Qin= 30L/min 時,區域2 和4顆粒沉積數量迅速增加,導致不同區域間沉積率差異顯著提高,區域2和3沉積率差異高達9.72%;同樣,呼吸量提升可獲得與增大粒徑相似的結果.

由圖9還可發現,顆粒在G11 ～ G14呼吸道健康側不同位置處的沉積會更均勻,對應病變側相同位置處的沉積率相差也不大,因而G11 ～ G14呼吸道內粒子沉積分布的整體不對稱性要稍小.此外,還發現G11 ～ G14呼吸道內粒子沉積的不對稱性對粒徑大小十分敏感,而對呼吸量的敏感性稍差,這一結果同G3 ～ G6呼吸道也有不同.

造成上述呼吸道內顆粒物沉積分布特征的差異,主要是因為顆粒在呼吸道內不同位置處的沉積機制不同.G3 ～ G6呼吸道內ηi均明顯高于ηg,這意味著,慣性碰撞在顆粒沉積過程中起主導作用.α的增大使得 ηi和 ηg都有下降,但對總的沉積率影響不大;相反,ηi和 ηg會隨著粒徑的增大而增加,且 ηi會有顯著的提高;此外,呼吸量增大會增強顆粒的慣性碰撞作用,但在一定程度也會削弱顆粒因重力作用發生的沉積.對比圖 9(a)和(b)不難發現,G11 ～ G14呼吸道內顆粒的沉積機制表現相反,重力在顆粒沉積過程中起主導作用.

2.4 顆粒的總沉積率

為考察阻塞性呼吸道內顆粒總沉積率的變化規律,分別統計了不同病變程度下 G3～G6和G11～G14呼吸道內顆粒物總沉積率(ηt)隨顆粒粒徑的變化,結果如圖 10所示.容易看出,不同病變程度下呼吸道表面 ηt都會隨著顆粒粒徑的增大而提高.

如圖 10(a)～(c)所示,勞動強度增強,ηt隨粒徑的上升趨勢會變得更為迅急,這是由于G3～G6呼吸道表面顆粒沉積機理由慣性碰撞主導,增大呼吸量將提升慣性作用和二次流碰撞機率.此外,α不同也會一定程度上影響 ηt的變化趨勢和大小,總體上,α越大,ηt越小.值得注意的是,不同勞動強度下,受 COPD影響造成的 ηt差異表現不同,低呼吸強度下(Qin=15L/min),ηt差異會隨著粒徑的增大逐漸凸顯,當dp=10μm 時,最大差異可達 20%.隨著呼吸量的提升,ηt差異主要體現在dp>5μm的大粒子上,當呼吸量增大至60L/min時,α對ηt的影響已十分有限.

圖10 不同阻塞率下G3 ～ G6和G11 ～ G14呼吸道內顆粒物的沉積效率Fig.10 Particle deposition efficiency in G3 ～ G6and G11 ～ G14respiratory tract under different levels of constriction

這一現象可以解釋為,對于慣性作用比較強的大顆粒,呼吸道阻塞很大程度降低了病變側顆粒的沉降數量,而健康側由于氣體流量增大,二次流對顆粒沉積的貢獻變大,顆粒沉積數量上升,當 Qin增大,健康側二次流效應會進一步增強,顆粒沉積密度進一步變大,使得呼吸道內ηt基本保持不變.

相比于G3～G6呼吸道,α對G11～G14呼吸道ηt的影響更為顯著,表現為 α增加可明顯降低顆粒總沉積率ηt,增大ηt差異,最大超過30%.這種差異的變化有兩個特點,一是其會隨著 dp的增大先逐漸變大而后減小;二是其在低勞動強度下不會出現明顯改變,隨著 Qin進一步提升,α 對 ηt的影響得到削弱.上述現象主要因為,對于 dp<7μm 的小顆粒,G11～G14呼吸道內重力沉積占主導,α越大,病變側顆粒的沉降數量越少;隨著 Qin和 dp的增大,顆粒在健康側慣性沉積得到增強,使得呼吸道內ηt差異減小.

3 結論

3.1 受COPD影響,呼吸道內流場分布表現出非對稱性,并在受阻塞位置的上下游形成滯止和回流區.α增大,滯止和回流區范圍擴大,病人局部缺氧越嚴重,當α = 0.8時,相對缺氧率可達90%以上.同等勞動強度下,病人的呼吸會更急促;勞動強度越強,病人絕對缺氧量變大,發生哮喘的可能性愈高.此外,受阻塞位置越深,局部缺氧癥狀越重,范圍也越集中.

3.2 呼吸道受阻不會改變顆粒的沉積機制,但對其沉積形式有顯著影響.具體為顆粒在健康側沉積數量增加,而在病變側則相反,顆粒呈不對稱分布.對于G3～G6呼吸道,α增大,勞動強度越強,dp越大,均會導致粒子沉積分布不對稱性的提高,當 α=0.8,Qin=30L/min,dp=5.0μm時,局部沉積率差異已高達9.72%.相比之下,G11 ～ G14呼吸道內粒子沉積分布會更均勻,整體不對稱性要更小,且對呼吸量的敏感性稍差.

3.3 阻塞型呼吸道內總沉降率 ηt會減小,且 α越大,ηt越小.呼吸道阻塞造成的 ηt差異在低勞動強度下更為明顯.對于G3～G6呼吸道,ηt差異主要體現在dp>5μm 的大顆粒上,最大可達 20%左右.此外,α對G11～G14呼吸內 ηt影響要更為突出,ηt差異也越顯著,并會隨dp的增大先增大而后減小,在dp=7μm附近達到最值,且最大差異超過30%.