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

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

(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 dynamicsdiscrete 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.225 kg/m3；ε 為顆粒所在網格的空隙率；u為空氣速度矢量；p為空氣壓力；g為重力；(▽u)tr為▽u的轉置矩陣；F為流體與顆粒物的相互作用力，定義為

式中，fD，i為顆粒i所受到的曳力；n為該網格內顆粒物的數量；ΔV為該網格的體積.

單個顆粒的運動方程為

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

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

式中，dp為顆粒物直徑；CD為曳力系數，計算式為

式中，ReN為顆粒物的雷諾數(ReN=·Dp/μ)；a1，a2，a3是由 ReN所確定的系數［26］.

根據Hertz-Mindlin接觸理論［27-28］計算顆粒物間的碰撞，法向彈性力fcn為

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

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

法向阻尼力fdcn為

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

切向彈性力fct為

式中，δt為切向位移；μp為靜摩擦系數；切向剛度St為

式中，等效剪切模量Geq為

式中，Gi和Gj分別為顆粒i和j的剪切模量.切向阻尼力fdct為

1.2 幾何模型

本文構建了基于Weibel-23級呼吸道結構［29］的G3～G53級肺部呼吸道模型，幾何尺寸如表1所示，分叉角為60°，模型結構及分區如圖1所示.圖中數字為分區號.為了分析顆粒物在呼吸道內的運動特性，對模型進行了劃分并編號.控制體網格由Gambit生成，模型網格總數為2172414個.

表1 呼吸道幾何尺寸 mm

圖1 G3～G5呼吸道結構圖

1.3 初始條件

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

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

文獻［24-25］表明，非穩態吸氣所形成的沉積形式可以用穩態吸氣替代，故本文采用了穩態吸氣進行模擬以減少計算時長.呼吸道入口速度場為拋物面型，平均速度為4 m/s，雷諾數Re約為1600，為層流流動.首先單獨計算流場直到收斂(殘差小于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 局部速度分布

圖2 統計方法示意圖

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

圖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中主要分布于壁面附近和渦的邊緣，而入口壁面附近顆粒物僅出現在渦的內部.

References)

［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］.Environmental Research，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］.Environ Health Perspect，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］.Journal of Zhejiang U-niversity:Engineering Science，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］.Journal of Engineering Thermophysics，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］.Journal of Engineering Thermophysics，2004，25(2):265-267.(in Chinese)

［8］Luo H Y，Liu Y，Yang X L.Particle deposition in obstructed airways［J］.Journal of Biomechanics，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］.Respiratory Physiology ＆ Neurobiology，2009，166(3):142-151.

［10］趙秀國.人體上呼吸道內氣流運動特性與氣溶膠沉積規律的研究［D］.天津:軍事醫學科學院衛生裝備研究所，2008.

［11］Kleinstreuer C，Zhang Z.Airflow and particle transport in the human respiratory system ［J］.Annual Review of Fluid Mechanic，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］.Chemical Engineering ＆ 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］.Powder Technology，2012，222:85-94.

［14］Favier J F，Abbaspour-Fard M H，Kremmer M.Modeling nonspherical particles using multisphere discrete elements［J］. Journal of Engineering Mechanics，2001，127(10):971-977.

［15］Vu-Quoc L，Zhang X，Walton O R.A 3-D discreteelement method for dry granular flows of ellipsoidal particles［J］.Computer Methods in Applied Mechanics and Engineering，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］.Biosystems Engineering，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］.Advanced Powder Technology，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］.Chemical Engineering Science，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］.Powder Technology，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］.Journal of Medical Engineering ＆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］.Journal of Aerosol Medicine，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］.Journal of Biomechanics，2001，34(5):659-669.

［23］Zhang Z，Kleinstreuer C，Kim C S.Gas-solid twophase flow in a triple bifurcation lung airway model［J］.International Journal of Multiphase Flow，2002，28(6):1021-1046.

［24］Zhang Z，Kleinstreuer C，Kim C S.Cyclic micronsize particle inhalation and deposition in a triple bifurcation lung airway model［J］.Journal of Aerosol Science，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］.Journal of Aerosol Science，2007，38(6):625-644.

［26］Morsi S A，Alexander A J.An investigation of particle trajectories in two-phase flow systems［J］.J Fluid Mech，1972，55(2):193-208.

［27］Mindlin R D，Deresiewicz H.Elastic spheres in contact under varying oblique forces［J］.Journal of Applied Mechanics，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.Morphometry of the human lung［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］.Powder Technology，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］.Journal of Fluid Mechanics，2001，435(1):55-80.