人體呼吸道內(nèi)可吸入顆粒徑向平面運(yùn)動(dòng)的CFD-DEM模擬

陳曉樂(lè) 鐘文琪 孫寶賓 金保昇 周獻(xiàn)光

(1東南大學(xué)能源熱轉(zhuǎn)換及其過(guò)程測(cè)控教育部重點(diǎn)實(shí)驗(yàn)室,南京 210096)(2東南大學(xué)醫(yī)學(xué)院,南京 210096)

可吸入顆粒物(PM10)作為大氣主要污染物,對(duì)人體健康尤其是呼吸系統(tǒng)和心血管具有嚴(yán)重危害.近年來(lái),我國(guó)霧霾天數(shù)猛增,據(jù)中國(guó)氣象局?jǐn)?shù)據(jù),2013年3月以來(lái)全國(guó)霧霾天數(shù)創(chuàng)52年新高[1].流行病學(xué)研究結(jié)果表明，可吸入顆粒物尤其是細(xì)顆粒物(PM2.5)與心血管和呼吸系統(tǒng)疾病的死亡率有密切聯(lián)系[2-3].我國(guó)首次發(fā)布的腫瘤發(fā)病率登記年報(bào)中亦指出,我國(guó)惡性腫瘤發(fā)病第一位的是肺癌,并且明顯呈現(xiàn)年輕化趨勢(shì)[4].

對(duì)可吸入顆粒物在呼吸道內(nèi)運(yùn)動(dòng)機(jī)理的研究有利于掌握由此引發(fā)的呼吸道疾病的病理,能夠?yàn)榭晌胨巹┭邪l(fā)提供重要理論依據(jù)和基礎(chǔ)數(shù)據(jù).然而我國(guó)在該領(lǐng)域起步較晚,開(kāi)展研究的單位也較少[5-10],因此我國(guó)人體可吸入顆粒物的深入研究迫在眉睫.

隨著呼吸道內(nèi)顆粒物運(yùn)動(dòng)與沉積研究的不斷深入,目前的研究重點(diǎn)已經(jīng)由局部呼吸道內(nèi)球形顆粒物轉(zhuǎn)向?qū)μ禺愋院粑馈⒎闻輩^(qū)、非球形顆粒物以及氣固交互作用等方面的分析[11].計(jì)算流體力學(xué)-離散元方法(computational fluid dynamics-discrete element method,CFD-DEM)與傳統(tǒng)的計(jì)算流體力學(xué)-離散相模型(computational fluid dynamics-discrete phase model,CFD-DPM)相比在顆粒-顆粒碰撞、顆粒旋轉(zhuǎn)、非球形顆粒物構(gòu)建等方面有突出優(yōu)勢(shì)[12-18],因此CFD-DEM方法在肺泡區(qū)和非球形顆粒物模擬方面可能具有較大優(yōu)勢(shì),文獻(xiàn)[19]對(duì)CFD-DEM方法進(jìn)行了驗(yàn)證,并針對(duì)顆粒物起始位置與顆粒物軌跡關(guān)系進(jìn)行了分析.

在過(guò)去的研究中,學(xué)者較多地關(guān)注了各種因素對(duì)顆粒物沉積的影響,如非對(duì)稱(chēng)速度進(jìn)口[20]、顆粒物入口分布[21]、彎曲進(jìn)口[22]、非穩(wěn)態(tài)呼吸[23-25]等,但對(duì)顆粒物在呼吸道內(nèi)輸運(yùn)過(guò)程中的分布和運(yùn)動(dòng)研究較少.多數(shù)模擬研究都采取了一個(gè)或幾個(gè)時(shí)間步長(zhǎng)間隔在入口處投入大量顆粒物,取部分截面觀察顆粒物在截面內(nèi)的分布情況.該方法可以獲得顆粒物在呼吸道內(nèi)分布的濃度和運(yùn)動(dòng)趨勢(shì),但因?yàn)楹粑纼?nèi)流場(chǎng)分布的不均勻使顆粒物速度并不一致,在某一時(shí)刻截面上出現(xiàn)的顆粒物是不同時(shí)間投入控制體內(nèi)的,所以不利于分析顆粒投入位置與輸運(yùn)過(guò)程中位置和速度之間的關(guān)系.

本文為進(jìn)一步積累可吸入藥劑研發(fā)的基礎(chǔ)理論依據(jù)，分析呼吸道內(nèi)可吸入顆粒物初始位置與輸運(yùn)特性之間的關(guān)系,構(gòu)建了基于Weibel-23級(jí)肺結(jié)構(gòu)的G3～G5級(jí)支氣管模型,利用CFD-DEM方法模擬呼吸道模型內(nèi)流場(chǎng)及顆粒物的運(yùn)動(dòng)與沉積,分析呼吸道入口中部和壁面附近顆粒物在G4與G5級(jí)分區(qū)內(nèi)的輸運(yùn)特性.

1 模擬方法

1.1 控制方程

假設(shè)流體為不可壓縮層流流動(dòng),控制方程為

(1)

·(εμ(u+(u)tr))+ρεg-F

(2)

式中,ρ為空氣密度，ρ=1.225 kg/m3;ε為顆粒所在網(wǎng)格的空隙率;u為空氣速度矢量;p為空氣壓力;g為重力;(u)tr為u的轉(zhuǎn)置矩陣;F為流體與顆粒物的相互作用力,定義為

(3)

式中,fD,i為顆粒i所受到的曳力;n為該網(wǎng)格內(nèi)顆粒物的數(shù)量;ΔV為該網(wǎng)格的體積.

單個(gè)顆粒的運(yùn)動(dòng)方程為

(4)

(5)

式中,mp,i為顆粒i的質(zhì)量;up,i和ωp,i為顆粒i的平移和轉(zhuǎn)動(dòng)速度矢量;fc,i為顆粒i所受到的碰撞力;Ii為顆粒的轉(zhuǎn)動(dòng)慣量;∑Ti為顆粒所受轉(zhuǎn)矩的矢量和.

曳力fD作為可吸入顆粒物在呼吸道內(nèi)運(yùn)動(dòng)和沉積的主導(dǎo)作用力，可由下式確定:

(6)

式中,dp為顆粒物直徑;CD為曳力系數(shù),計(jì)算式為

(7)

根據(jù)Hertz-Mindlin接觸理論[27-28]計(jì)算顆粒物間的碰撞,法向彈性力fcn為

(8)

式中,δn為法向的變形量;等效楊氏模量Eeq和等效半徑Req定義為

(9)

(10)

式中,Ei,υi,Ri和Ej,υj,Rj分別為顆粒i和j的楊氏模量、泊松比和半徑.

(11)

(12)

(13)

式中,mi和mj分別為顆粒i和j的質(zhì)量.

切向彈性力fct為

(14)

式中,δt為切向位移;μp為靜摩擦系數(shù);切向剛度St為

(15)

式中,等效剪切模量Geq為

(16)

式中,Gi和Gj分別為顆粒i和j的剪切模量.

(17)

1.2 幾何模型

本文構(gòu)建了基于Weibel-23級(jí)呼吸道結(jié)構(gòu)[29]的G3～G5 3級(jí)肺部呼吸道模型,幾何尺寸如表1所示,分叉角為60°,模型結(jié)構(gòu)及分區(qū)如圖1所示.

表1 呼吸道幾何尺寸 mm

圖1 G3～G5呼吸道結(jié)構(gòu)圖

圖中數(shù)字為分區(qū)號(hào).為了分析顆粒物在呼吸道內(nèi)的運(yùn)動(dòng)特性,對(duì)模型進(jìn)行了劃分并編號(hào).控制體網(wǎng)格由Gambit生成,模型網(wǎng)格總數(shù)為2 172 414個(gè).

1.3 初始條件

流場(chǎng)計(jì)算的時(shí)間步長(zhǎng)為0.1 ms,顆粒物運(yùn)動(dòng)的時(shí)間步長(zhǎng)設(shè)置為Rayleigh時(shí)間TR[30]的50%,TR定義為

(18)

式中,Rp,ρp,Gp和υp分別為顆粒物的半徑、密度、剪切模量和泊松比.

文獻(xiàn)[24-25]表明,非穩(wěn)態(tài)吸氣所形成的沉積形式可以用穩(wěn)態(tài)吸氣替代,故本文采用了穩(wěn)態(tài)吸氣進(jìn)行模擬以減少計(jì)算時(shí)長(zhǎng).呼吸道入口速度場(chǎng)為拋物面型,平均速度為4 m/s,雷諾數(shù)Re約為1 600,為層流流動(dòng).首先單獨(dú)計(jì)算流場(chǎng)直到收斂(殘差小于10-4),之后在入口處隨即生成1.0×104個(gè)顆粒物,顆粒物隨氣流運(yùn)動(dòng),直至全部顆粒物沉積于模型表面或從出口處流出.相關(guān)參數(shù)的取值如表2所示.本文的CFD-DEM模擬利用FLUENT軟件求解流場(chǎng)分布,采用EDEM軟件計(jì)算顆粒相運(yùn)動(dòng)并進(jìn)行耦合,模擬完成后使用C++編程方法統(tǒng)計(jì)和分析顆粒物數(shù)據(jù).

表2 參數(shù)取值

1.4 顆粒物統(tǒng)計(jì)方法

由于DEM方法計(jì)算量遠(yuǎn)大于DPM方法,如采取傳統(tǒng)大量投入顆粒物取截面的方法將大大增加計(jì)算時(shí)間和存儲(chǔ)空間.本模擬中典型工況計(jì)算時(shí)長(zhǎng)約為7 d,產(chǎn)生數(shù)據(jù)量在300～500 GB.因此本文采取一次投入顆粒,分時(shí)間段區(qū)域取樣疊加方法進(jìn)行顆粒物位置與速度的統(tǒng)計(jì),如圖2所示.

2 結(jié)果與討論

2.1 局部速度分布

圖3為呼吸道中心平面G4和G5空氣速度云圖,在拋物面形速度進(jìn)口條件下,G3中部高速氣流沖向G4級(jí)分叉，平均分為2支,靠近分叉內(nèi)側(cè)為高速區(qū),而A′側(cè)速度較低,發(fā)生了邊界層分離.A側(cè)高速氣流帶動(dòng)A′側(cè)向下游運(yùn)動(dòng),其速度中心仍然靠近A側(cè),所以在進(jìn)入G5時(shí)更多的空氣進(jìn)入C-C′方向的支氣管.由于進(jìn)入C-C′方向的氣流速度依然較高,在C′側(cè)同樣產(chǎn)生了分離.B-B′支氣管速度分布較C-C′均勻,但受上游A側(cè)流速較高的影響,B側(cè)氣流速度略高于B′側(cè).

圖2 統(tǒng)計(jì)方法示意圖

圖3 呼吸道中心平面速度云圖

2.2 顆粒物位置與徑向速度分布

在拋物面形空氣進(jìn)口條件作用下,入口中部的顆粒物速度較快,大多在0～0.045 s內(nèi)通過(guò)呼吸道,靠近入口壁面的顆粒物速度較慢,大多在0.22～0.28 s內(nèi)通過(guò)呼吸道,x-z平面的顆粒物位置隨時(shí)間的變化可見(jiàn)文獻(xiàn)[31].圖4(a)和(b)為入口中部顆粒和靠近壁面顆粒物通過(guò)分區(qū)1時(shí)的位置和徑向速度分布.由圖4(a)可以發(fā)現(xiàn),顆粒物在通過(guò)入口中部時(shí)在二次流的作用下形成了近似二次流對(duì)稱(chēng)渦的速度分布.截面方向上顆粒物在呼吸道中心處速度最快,壁面附近速度較慢,且壁面附近呼吸道分叉A側(cè)的速度高于A′側(cè)速度.從濃度上看,壁面附近顆粒分布較多,中部濃度較壁面低,分布均勻；A側(cè)方向壁面附近幾乎沒(méi)有顆粒分布,這可能是由于此方向上顆粒物大多沉積于上游分叉處或受到分叉處逆壓力梯度的影響被排擠向呼吸道中部.圖4(b)中靠近入口壁面的顆粒物雖然也受到二次流的影響,但垂直于管流的速度分量遠(yuǎn)小于圖4(a)入口中部的顆粒,同時(shí)該顆粒物的分布也與圖4(a)中有明顯區(qū)別,集中在呼吸道中央偏分叉外側(cè)分區(qū)，而在壁面附近沒(méi)有分布.Comer等[32]的模擬也發(fā)現(xiàn),在G4壁面附近有顆粒物沉積,由此可以推斷該沉積源于入口靠近中部的顆粒物.

圖4 分區(qū)1內(nèi)顆粒物位置與徑向速度分布

圖5(a)和(b)分別為入口中部顆粒和靠近壁面顆粒物通過(guò)分區(qū)3時(shí)的位置和徑向速度分布.由圖5(a)可以發(fā)現(xiàn)，該分區(qū)內(nèi)二次流受上游影響形成了復(fù)雜的渦結(jié)構(gòu).在B′側(cè)虛線逆時(shí)針二次流形成的渦流中顆粒物濃度較高,右上方順時(shí)針運(yùn)動(dòng)的渦流中濃度次之,占據(jù)管流最大面積的左側(cè)逆時(shí)針渦流中顆粒物濃度最低.圖5(b)中靠近壁面的顆粒物僅參與了呼吸道中部?jī)?nèi)側(cè)渦流的運(yùn)動(dòng),其垂直于管流的速度分量大小與圖5(a)中參與該位置渦流運(yùn)動(dòng)的入口中部顆粒相近.

圖5 分區(qū)3內(nèi)顆粒物位置與徑向速度分布

圖6(a)和(b)分別為入口中部顆粒和靠近壁面顆粒物通過(guò)分區(qū)5時(shí)的位置和徑向速度分布.圖6中由于氣流的速度中心偏向呼吸道分叉的C′側(cè),二次流向上下排擠空氣,在靠近壁面處分向兩側(cè),各形成2個(gè)漩渦.圖6(a)入口中部的顆粒物主要集中于壁面和中部渦的邊緣,而圖6(b)靠近入口壁面顆粒只分布于呼吸道中心平面附近渦的內(nèi)部.

圖6 分區(qū)5內(nèi)顆粒物位置與徑向速度分布

3 結(jié)論

1) 模型入口中部的顆粒物在G4級(jí)壁面附近濃度較高,這也是G4級(jí)呼吸道壁面產(chǎn)生沉積的原因,入口壁面附近的顆粒物主要分布在G4級(jí)中央偏分叉外側(cè)分區(qū),在壁面附近沒(méi)有分布.

2) 入口中部的顆粒物在G5級(jí)外側(cè)分支分區(qū)3中呼吸道外側(cè)濃度分布較高,入口壁面附近的顆粒物位于呼吸道中部.

3) 入口中部的顆粒物在G5級(jí)內(nèi)側(cè)分支分區(qū)5中主要分布于壁面附近和渦的邊緣,而入口壁面附近顆粒物僅出現(xiàn)在渦的內(nèi)部.

)

[1] 中國(guó)氣象局. 中國(guó)氣象局2013年4月新聞發(fā)布會(huì)[EB/OL]. (2013-04-03)[2013-04-10]. http://www.scio.gov.cn/xwfbh/gbwxwfbh/fbh/201304/t1307058.htm.

[2] Ostro B D,Hurley S,Lipsett M J. Air pollution and daily mortality in the coachella valley,California: a study of pm10 dominated by coarse particles[J].EnvironmentalResearch,1999,81(3): 231-238.

[3] Schwartz J,Norris G,Larson T,et al. Episodes of high coarse particle concentrations are not associated with increased mortality[J].EnvironHealthPerspect,1999,107(5): 339-342.

[4] 郝捷,陳萬(wàn)青. 2012中國(guó)腫瘤登記年報(bào) [M]. 北京: 軍事醫(yī)學(xué)科學(xué)出版社,2012: 15-33.

[5] 曾敏捷,胡桂林,樊建人. 微顆粒在人體上呼吸道中運(yùn)動(dòng)沉積的數(shù)值模擬[J]. 浙江大學(xué)學(xué)報(bào):工學(xué)版,2006,40(7):1164-1167,1256.

Zeng Minjie,Hu Guilin,Fan Jianren. Numerical simulation of micro-particle movement and deposition in human upper respiratory tract[J].JournalofZhejiangUniversity:EngineeringScience,2006,40(7): 1164-1167,1256. (in Chinese)

[6] 林江,胡桂林,樊建人. 氣管支氣管樹(shù)內(nèi)氣流和顆粒運(yùn)動(dòng)的大渦模擬[J]. 工程熱物理學(xué)報(bào),2007,28(5):805-807.

Lin Jiang,Hu Guilin,Fan Jianren. Large eddy simulation of airflow and particle motion in tracheobronchial tree[J].JournalofEngineeringThermophysics,2007,28(5): 805-807. (in Chinese)

[7] 由長(zhǎng)福,李光輝,祁海鷹,等. 可吸入顆粒物近壁運(yùn)動(dòng)的直接數(shù)值模擬[J]. 工程熱物理學(xué)報(bào),2004,25(2):265-267.

You Changfu,Li Guanghui,Qi Haiying,et al. DNS of inhalable particle motion in channel flow[J].JournalofEngineeringThermophysics,2004,25(2): 265-267. (in Chinese)

[8] Luo H Y,Liu Y,Yang X L. Particle deposition in obstructed airways[J].JournalofBiomechanics,2007,40(14): 3096-3104.

[9] Wang S M,Inthavong K,Wen J,et al. Comparison of micron and nanoparticle deposition patterns in a realistic human nasal cavity[J].RespiratoryPhysiology&Neurobiology,2009,166(3): 142-151.

[10] 趙秀國(guó). 人體上呼吸道內(nèi)氣流運(yùn)動(dòng)特性與氣溶膠沉積規(guī)律的研究[D]. 天津:軍事醫(yī)學(xué)科學(xué)院衛(wèi)生裝備研究所,2008.

[11] Kleinstreuer C,Zhang Z. Airflow and particle transport in the human respiratory system [J].AnnualReviewofFluidMechanic,2010,42: 301-334.

[12] Ren Bing,Zhong Wengqi,Jin Baosheng,et al. Modeling of gas-particle turbulent flow in spout fluid bed by computational fluid dynamics with discrete element method[J].ChemicalEngineering&Technology,2011,34(12): 2059-2068.

[13] Ren Bing,Shao Yingjuan,Zhong Wenqi,et al. Investigation of mixing behaviors in a spouted bed with different density particles using discrete element method [J].PowderTechnology,2012,222: 85-94.

[14] Favier J F,Abbaspour-Fard M H,Kremmer M. Modeling nonspherical particles using multisphere discrete elements[J].JournalofEngineeringMechanics,2001,127(10): 971-977.

[15] Vu-Quoc L,Zhang X,Walton O R. A 3-D discrete-element method for dry granular flows of ellipsoidal particles[J].ComputerMethodsinAppliedMechanicsandEngineering,2000,187(3): 483-528.

[16] Abbaspour-Fard M H. Theoretical validation of a multi-sphere,discrete element model suitable for biomaterials handling simulation[J].BiosystemsEngineering,2004,88(2): 153-161.

[17] Abou-Chakra H,Baxter J,Tuzun U. Three-dimensional particle shape descriptors for computer simulation of non-spherical particulate assemblies[J].AdvancedPowderTechnology,2004,15(1): 63-77.

[18] Kodam M,Curtis J,Hancock B,et al. Discrete element method modeling of bi-convex pharmaceutical tablets: contact detection algorithms and validation[J].ChemicalEngineeringScience,2012,69(1): 587-601.

[19] Chen Xiaole,Zhong Wenqi,Sun Baobin,et al. Study on gas/solid flow in an obstructed pulmonary airway with transient flow based on CFD-DPM approach[J].PowderTechnology,2012,217: 252-260.

[20] Zhang Z,Kleinstreuer C,Kim C S. Effects of asymmetric branch flow rates on aerosol deposition in bifurcating airways[J].JournalofMedicalEngineering&Technology,2000,24(5): 192-202.

[21] Zhang Z,Kleinstreuer C. Effect of particle inlet distributions on deposition in a triple bifurcation lung airway model[J].JournalofAerosolMedicine,2001,14(1): 13-29.

[22] Zhang Z,Kleinstreuer C,Kim C S. Effects of curved inlet tubes on air flow and particle deposition in bifurcating lung models[J].JournalofBiomechanics,2001,34(5): 659-669.

[23] Zhang Z,Kleinstreuer C,Kim C S. Gas-solid two-phase flow in a triple bifurcation lung airway model[J].InternationalJournalofMultiphaseFlow,2002,28(6): 1021-1046.

[24] Zhang Z,Kleinstreuer C,Kim C S. Cyclic micron-size particle inhalation and deposition in a triple bifurcation lung airway model[J].JournalofAerosolScience,2002,33(2): 257-281.

[25] Li Z,Kleinstreuer C,Zhang Z. Particle deposition in the human tracheobronchial airways due to transient inspiratory flow patterns[J].JournalofAerosolScience,2007,38(6): 625-644.

[26] Morsi S A,Alexander A J. An investigation of particle trajectories in two-phase flow systems[J].JFluidMech,1972,55(2): 193-208.

[27] Mindlin R D,Deresiewicz H. Elastic spheres in contact under varying oblique forces [J].JournalofAppliedMechanics,1953,20(3): 327-344.

[28] Raji A O. Discrete element modelling of the deformation of bulk agricultural particulates [D]. Newcastle upon Tyne,UK: Newcastle University,1999.

[29] Weibel E R.Morphometryofthehumanlung[M]. Salt Lake City,USA: Academic Press,1963.

[30] Ning Z. Elasto-plastic impact of fine particles and fragmentation of small agglomerates [D]. Birmingham,UK: Aston University,1999.

[31] Chen Xiaole,Zhong Wenqi,Zhou Xiaoguang,et al. CFD-DEM simulation of particle transport and deposition in pulmonary airway[J].PowderTechnology,2012,228: 309-318.

[32] Comer J K,Kleinstreuer C,Kim C S. Flow structures and particle deposition patterns in double-bifurcation airway models. part 2. Aerosol transport and deposition[J].JournalofFluidMechanics,2001,435(1): 55-80.