陈晓乐 钟文琪 孙宝宾 金保昇 周献光

(1东南大学能源热转换及其过程测控教育部重点实验室,南京 210096)(2东南大学医学院,南京 210096)

可吸入颗粒物(PM10)作为大气主要污染物,对人体健康尤其是呼吸系统和心血管具有严重危害.近年来,我国雾霾天数猛增,据中国气象局数据,2013年3月以来全国雾霾天数创52年新高[1].流行病学研究结果表明，可吸入颗粒物尤其是细颗粒物(PM2.5)与心血管和呼吸系统疾病的死亡率有密切联系[2-3].我国首次发布的肿瘤发病率登记年报中亦指出,我国恶性肿瘤发病第一位的是肺癌,并且明显呈现年轻化趋势[4].

对可吸入颗粒物在呼吸道内运动机理的研究有利于掌握由此引发的呼吸道疾病的病理,能够为可吸入药剂研发提供重要理论依据和基础数据.然而我国在该领域起步较晚,开展研究的单位也较少[5-10],因此我国人体可吸入颗粒物的深入研究迫在眉睫.

随着呼吸道内颗粒物运动与沉积研究的不断深入,目前的研究重点已经由局部呼吸道内球形颗粒物转向对特异性呼吸道、肺泡区、非球形颗粒物以及气固交互作用等方面的分析[11].计算流体力学-离散元方法(computational fluid dynamics-discrete element method,CFD-DEM)与传统的计算流体力学-离散相模型(computational fluid dynamics-discrete phase model,CFD-DPM)相比在颗粒-颗粒碰撞、颗粒旋转、非球形颗粒物构建等方面有突出优势[12-18],因此CFD-DEM方法在肺泡区和非球形颗粒物模拟方面可能具有较大优势,文献[19]对CFD-DEM方法进行了验证,并针对颗粒物起始位置与颗粒物轨迹关系进行了分析.

在过去的研究中,学者较多地关注了各种因素对颗粒物沉积的影响,如非对称速度进口[20]、颗粒物入口分布[21]、弯曲进口[22]、非稳态呼吸[23-25]等,但对颗粒物在呼吸道内输运过程中的分布和运动研究较少.多数模拟研究都采取了一个或几个时间步长间隔在入口处投入大量颗粒物,取部分截面观察颗粒物在截面内的分布情况.该方法可以获得颗粒物在呼吸道内分布的浓度和运动趋势,但因为呼吸道内流场分布的不均匀使颗粒物速度并不一致,在某一时刻截面上出现的颗粒物是不同时间投入控制体内的,所以不利于分析颗粒投入位置与输运过程中位置和速度之间的关系.

本文为进一步积累可吸入药剂研发的基础理论依据，分析呼吸道内可吸入颗粒物初始位置与输运特性之间的关系,构建了基于Weibel-23级肺结构的G3～G5级支气管模型,利用CFD-DEM方法模拟呼吸道模型内流场及颗粒物的运动与沉积,分析呼吸道入口中部和壁面附近颗粒物在G4与G5级分区内的输运特性.

1 模拟方法

1.1 控制方程

假设流体为不可压缩层流流动,控制方程为

(1)

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

(2)

式中,ρ为空气密度，ρ=1.225 kg/m3;ε为颗粒所在网格的空隙率;u为空气速度矢量;p为空气压力;g为重力;(u)tr为u的转置矩阵;F为流体与颗粒物的相互作用力,定义为

(3)

式中,fD,i为颗粒i所受到的曳力;n为该网格内颗粒物的数量;ΔV为该网格的体积.

单个颗粒的运动方程为

(4)

(5)

式中,mp,i为颗粒i的质量;up,i和ωp,i为颗粒i的平移和转动速度矢量;fc,i为颗粒i所受到的碰撞力;Ii为颗粒的转动惯量;∑Ti为颗粒所受转矩的矢量和.

曳力fD作为可吸入颗粒物在呼吸道内运动和沉积的主导作用力，可由下式确定:

(6)

式中,dp为颗粒物直径;CD为曳力系数,计算式为

(7)

根据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为静摩擦系数;切向刚度St为

(15)

式中,等效剪切模量Geq为

(16)

式中,Gi和Gj分别为颗粒i和j的剪切模量.

(17)

1.2 几何模型

本文构建了基于Weibel-23级呼吸道结构[29]的G3～G5 3级肺部呼吸道模型,几何尺寸如表1所示,分叉角为60°,模型结构及分区如图1所示.

表1 呼吸道几何尺寸 mm

图1 G3～G5呼吸道结构图

图中数字为分区号.为了分析颗粒物在呼吸道内的运动特性,对模型进行了划分并编号.控制体网格由Gambit生成,模型网格总数为2 172 414个.

1.3 初始条件

流场计算的时间步长为0.1 ms,颗粒物运动的时间步长设置为Rayleigh时间TR[30]的50%,TR定义为

(18)

式中,Rp,ρp,Gp和υp分别为颗粒物的半径、密度、剪切模量和泊松比.

文献[24-25]表明,非稳态吸气所形成的沉积形式可以用稳态吸气替代,故本文采用了稳态吸气进行模拟以减少计算时长.呼吸道入口速度场为抛物面型,平均速度为4 m/s,雷诺数Re约为1 600,为层流流动.首先单独计算流场直到收敛(残差小于10-4),之后在入口处随即生成1.0×104个颗粒物,颗粒物随气流运动,直至全部颗粒物沉积于模型表面或从出口处流出.相关参数的取值如表2所示.本文的CFD-DEM模拟利用FLUENT软件求解流场分布,采用EDEM软件计算颗粒相运动并进行耦合,模拟完成后使用C++编程方法统计和分析颗粒物数据.

表2 参数取值

1.4 颗粒物统计方法

由于DEM方法计算量远大于DPM方法,如采取传统大量投入颗粒物取截面的方法将大大增加计算时间和存储空间.本模拟中典型工况计算时长约为7 d,产生数据量在300～500 GB.因此本文采取一次投入颗粒,分时间段区域取样叠加方法进行颗粒物位置与速度的统计,如图2所示.

2 结果与讨论

2.1 局部速度分布

图3为呼吸道中心平面G4和G5空气速度云图,在抛物面形速度进口条件下,G3中部高速气流冲向G4级分叉，平均分为2支,靠近分叉内侧为高速区,而A′侧速度较低,发生了边界层分离.A侧高速气流带动A′侧向下游运动,其速度中心仍然靠近A侧,所以在进入G5时更多的空气进入C-C′方向的支气管.由于进入C-C′方向的气流速度依然较高,在C′侧同样产生了分离.B-B′支气管速度分布较C-C′均匀,但受上游A侧流速较高的影响,B侧气流速度略高于B′侧.

图2 统计方法示意图

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

2.2 颗粒物位置与径向速度分布

在抛物面形空气进口条件作用下,入口中部的颗粒物速度较快,大多在0～0.045 s内通过呼吸道,靠近入口壁面的颗粒物速度较慢,大多在0.22～0.28 s内通过呼吸道,x-z平面的颗粒物位置随时间的变化可见文献[31].图4(a)和(b)为入口中部颗粒和靠近壁面颗粒物通过分区1时的位置和径向速度分布.由图4(a)可以发现,颗粒物在通过入口中部时在二次流的作用下形成了近似二次流对称涡的速度分布.截面方向上颗粒物在呼吸道中心处速度最快,壁面附近速度较慢,且壁面附近呼吸道分叉A侧的速度高于A′侧速度.从浓度上看,壁面附近颗粒分布较多,中部浓度较壁面低,分布均匀；A侧方向壁面附近几乎没有颗粒分布,这可能是由于此方向上颗粒物大多沉积于上游分叉处或受到分叉处逆压力梯度的影响被排挤向呼吸道中部.图4(b)中靠近入口壁面的颗粒物虽然也受到二次流的影响,但垂直于管流的速度分量远小于图4(a)入口中部的颗粒,同时该颗粒物的分布也与图4(a)中有明显区别,集中在呼吸道中央偏分叉外侧分区，而在壁面附近没有分布.Comer等[32]的模拟也发现,在G4壁面附近有颗粒物沉积,由此可以推断该沉积源于入口靠近中部的颗粒物.

图4 分区1内颗粒物位置与径向速度分布

图5(a)和(b)分别为入口中部颗粒和靠近壁面颗粒物通过分区3时的位置和径向速度分布.由图5(a)可以发现，该分区内二次流受上游影响形成了复杂的涡结构.在B′侧虚线逆时针二次流形成的涡流中颗粒物浓度较高,右上方顺时针运动的涡流中浓度次之,占据管流最大面积的左侧逆时针涡流中颗粒物浓度最低.图5(b)中靠近壁面的颗粒物仅参与了呼吸道中部内侧涡流的运动,其垂直于管流的速度分量大小与图5(a)中参与该位置涡流运动的入口中部颗粒相近.

图5 分区3内颗粒物位置与径向速度分布

图6(a)和(b)分别为入口中部颗粒和靠近壁面颗粒物通过分区5时的位置和径向速度分布.图6中由于气流的速度中心偏向呼吸道分叉的C′侧,二次流向上下排挤空气,在靠近壁面处分向两侧,各形成2个漩涡.图6(a)入口中部的颗粒物主要集中于壁面和中部涡的边缘,而图6(b)靠近入口壁面颗粒只分布于呼吸道中心平面附近涡的内部.

图6 分区5内颗粒物位置与径向速度分布

3 结论

1) 模型入口中部的颗粒物在G4级壁面附近浓度较高,这也是G4级呼吸道壁面产生沉积的原因,入口壁面附近的颗粒物主要分布在G4级中央偏分叉外侧分区,在壁面附近没有分布.

2) 入口中部的颗粒物在G5级外侧分支分区3中呼吸道外侧浓度分布较高,入口壁面附近的颗粒物位于呼吸道中部.

3) 入口中部的颗粒物在G5级内侧分支分区5中主要分布于壁面附近和涡的边缘,而入口壁面附近颗粒物仅出现在涡的内部.

)

[1] 中国气象局. 中国气象局2013年4月新闻发布会[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]. 北京: 军事医学科学出版社,2012: 15-33.

[5] 曾敏捷,胡桂林,樊建人. 微颗粒在人体上呼吸道中运动沉积的数值模拟[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] 由长福,李光辉,祁海鹰,等. 可吸入颗粒物近壁运动的直接数值模拟[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] 赵秀国. 人体上呼吸道内气流运动特性与气溶胶沉积规律的研究[D]. 天津:军事医学科学院卫生装备研究所,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.