人體呼吸道內可吸入顆粒徑向平面運動的CFD-DEM模擬

陳曉樂 鐘文琪 孫寶賓 金保昇 周獻光

(1東南大學能源熱轉換及其過程測控教育部重點實驗室,南京 210096)(2東南大學醫(yī)學院,南京 210096)

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

對可吸入顆粒物在呼吸道內運動機理的研究有利于掌握由此引發(fā)的呼吸道疾病的病理,能夠為可吸入藥劑研發(fā)提供重要理論依據(jù)和基礎數(shù)據(jù).然而我國在該領域起步較晚,開展研究的單位也較少[5-10],因此我國人體可吸入顆粒物的深入研究迫在眉睫.

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

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

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

1 模擬方法

1.1 控制方程

假設流體為不可壓縮層流流動,控制方程為

(1)

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

(2)

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

(3)

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

單個顆粒的運動方程為

(4)

(5)

式中,mp,i為顆粒i的質量;up,i和ωp,i為顆粒i的平移和轉動速度矢量;fc,i為顆粒i所受到的碰撞力;Ii為顆粒的轉動慣量;∑Ti為顆粒所受轉矩的矢量和.

曳力fD作為可吸入顆粒物在呼吸道內運動和沉積的主導作用力，可由下式確定:

(6)

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

(7)

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

(8)

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

(9)

(10)

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

(11)

(12)

(13)

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

切向彈性力fct為

(14)

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

(15)

式中,等效剪切模量Geq為

(16)

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

(17)

1.2 幾何模型

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

表1 呼吸道幾何尺寸 mm

圖1 G3～G5呼吸道結構圖

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

1.3 初始條件

流場計算的時間步長為0.1 ms,顆粒物運動的時間步長設置為Rayleigh時間TR[30]的50%,TR定義為

(18)

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

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

表2 參數(shù)取值

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

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

2 結果與討論

2.1 局部速度分布

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

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

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

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

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

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

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

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

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

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

3 結論

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

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

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

)

[1] 中國氣象局. 中國氣象局2013年4月新聞發(fā)布會[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] 郝捷,陳萬青. 2012中國腫瘤登記年報 [M]. 北京: 軍事醫(yī)學科學出版社,2012: 15-33.

[5] 曾敏捷,胡桂林,樊建人. 微顆粒在人體上呼吸道中運動沉積的數(shù)值模擬[J]. 浙江大學學報:工學版,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] 林江,胡桂林,樊建人. 氣管支氣管樹內氣流和顆粒運動的大渦模擬[J]. 工程熱物理學報,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] 由長福,李光輝,祁海鷹,等. 可吸入顆粒物近壁運動的直接數(shù)值模擬[J]. 工程熱物理學報,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ī)律的研究[D]. 天津:軍事醫(yī)學科學院衛(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.