嚴(yán)微微，孫玉勝，劉星利

(中國計量大學(xué) 計量測試工程學(xué)院，浙江 杭州 310018)

上呼吸道流場及顆粒物沉積規(guī)律的模擬

嚴(yán)微微，孫玉勝，劉星利

(中國計量大學(xué) 計量測試工程學(xué)院，浙江 杭州 310018)

基于醫(yī)學(xué)CT掃描圖像得到的真實人體上呼吸道模型，采用大渦模擬方法計算了上呼吸道內(nèi)部的氣流流動，分析了氣流流速和渦旋的特征.在拉格朗日框架下計算了可吸入顆粒物的沉積，對顆粒的局部沉積位置進行了直觀顯示.研究結(jié)果表明，舌咽部位渦旋和回流最多，狹窄結(jié)構(gòu)處流速最大，且產(chǎn)生了氣流噴射現(xiàn)象；大部分的顆粒沉積在鼻腔和喉部，顆粒粒徑對沉積率有很大影響.

上呼吸道；大渦模擬；氣流特征；顆粒沉積

霧霾天氣對人們的日常生活影響惡劣[1]，可吸入顆粒物會對人們的呼吸系統(tǒng)產(chǎn)生嚴(yán)重影響，導(dǎo)致上呼吸道感染、慢性阻塞性肺病、支氣管哮喘等疾病的發(fā)病率明顯升高[2].深入了解人體上呼吸道內(nèi)部的氣流特征及可吸入顆粒物在其中的沉積規(guī)律，不僅可以幫助判斷病人的病情，還可以為呼吸系統(tǒng)疾病的致病機理提供可靠的依據(jù).

目前，已有一些關(guān)于呼吸道內(nèi)氣流流場和顆粒運動的計算模擬研究，實驗方面：Cheng等根據(jù)某志愿者的呼吸道尺寸建立了上呼吸道三維模型[3]，并選用9種不同粒徑的顆粒進行了呼吸實驗，發(fā)現(xiàn)顆粒大小與呼吸強度對顆粒的沉積率有很大影響.Grgic等運用CT圖像生成技術(shù)，對人體呼吸道進行數(shù)據(jù)采集[4]，建立了三維呼吸道模型，實驗統(tǒng)計了顆粒在口腔、咽部、喉部和氣管的沉積百分比，發(fā)現(xiàn)顆粒大量沉積在喉部與氣管處，顆粒慣性是影響顆粒沉積率的關(guān)鍵因素.數(shù)值模擬方面：Stapleton等基于Grgic的呼吸道模型[5]，采用k-ε湍流模型模擬上呼吸道中的氣流，并將Lagrange跟蹤模型應(yīng)用于顆粒運動，結(jié)果顯示湍流模型能夠很好地預(yù)測低呼吸強度時上呼吸道中的氣流運動.Zhang等基于Cheng的呼吸道模型[6]，采用k-ε湍流模型計算流場，分別用歐拉-歐拉法和歐拉-拉格朗日方法模擬追蹤了呼吸道內(nèi)顆粒的運動情況，在不同入口氣流邊界條件下分別比較了從口腔到前三級支氣管內(nèi)，微顆粒和納米顆粒的沉積模式和局部濃度.縱觀國內(nèi)文獻，他們的幾何模型大多參照了Stapleton等的簡化模型[7-8]；所使用的計算模型多為低雷諾數(shù)k-ω湍流模型[9]和RNGk-ε模型[10]，發(fā)現(xiàn)顆粒沉積特性與其質(zhì)量相關(guān)[11].本文在此基礎(chǔ)上，基于醫(yī)學(xué)CT掃描的原始圖像，重構(gòu)得到真實人體上呼吸道的三維幾何模型；并且采用更為精確的可得到瞬時流場特性的大渦模擬方法求解氣流運動特性，利用與實際情況相符合的擬合函數(shù)仿真隨時間變化的呼吸強度；而且還在拉格朗日框架下跟蹤顆粒物的運動情況，討論了顆粒物在人體呼吸道內(nèi)部的沉積規(guī)律.

1 模型描述

1.1 幾何模型

原始圖像數(shù)據(jù)是由CT掃描一位中國男性患者的胸部得到的.圖像由分辨率為0.7×0.7 mm2的軸向平面制成，層厚為0.625 mm.上呼吸道模型三維點云數(shù)據(jù)由圖像處理軟件Mimics重建，其幾何模型如圖1.

圖1 上呼吸道三維結(jié)構(gòu)模型Figure 1 3D model of upper respiratory tract

對此模型可以進行網(wǎng)格劃分和計算流體力學(xué)分析.其中，鼻腔入口及喉部的延長段是為了更好地處理邊界條件，方便計算.從圖中可以看出，該呼吸道幾何結(jié)較為復(fù)雜，鼻腔和咽喉部位突縮結(jié)構(gòu)明顯，因此與之前文獻中所使用的簡化模型相比，此物理模型更具準(zhǔn)確性.

1.2 控制方程

由于呼吸道內(nèi)結(jié)構(gòu)復(fù)雜，存在急劇收縮和擴展的面，本文采用大渦模擬的方法描述氣流場，用以更加準(zhǔn)確地得到內(nèi)部層流和復(fù)雜湍流的信息.在該方法中，對變量x進行濾波后的方程為

(1)

其中，V是計算單元的體積，濾波方程定義為

(2)

這一濾波過程有效地過濾了小于過濾寬度或網(wǎng)格間距的渦旋，因此經(jīng)過濾波處理的Navier-Stokes方程為：

連續(xù)性方程

(3)

動量方程

對于顆粒相，認(rèn)定其為稀相，且其密度遠(yuǎn)大于氣相，不考慮顆粒之間的相互作用，只考慮顆粒受到的重力及Stokes力和Saffman升力，因此單顆粒的運動方程可以表達為

(5)

式中mp是顆粒質(zhì)量，vi是顆粒的速度，dp是顆粒的直徑，CD是顆粒的阻力系數(shù)，g是重力加速度.右邊三項依次是單顆粒受到的Stokes阻力、重力和Saffman升力.

1.3 數(shù)值仿真方法

在進行數(shù)值仿真計算時，假設(shè)呼吸道壁面為剛性光滑面，計算域為鼻腔入口到會厭底部，計算軟件采用商業(yè)軟件Fluent，使用非結(jié)構(gòu)化網(wǎng)格，共劃分網(wǎng)格120萬左右，在顎咽即結(jié)構(gòu)最狹窄處進行了網(wǎng)格加密處理.使用大渦模擬(Large Eddy Simulation, LES)湍流模型，速度壓力耦合的求解采用PISO算法，瞬態(tài)方程采用二階隱式的方法求解.鼻腔入口采用速度入口邊界條件，出口采用壓力出口邊界條件，壓力為零.入口氣流流量隨時間周期變化，計算時間步長0.001 s，自行編制C語言函數(shù)模擬呼吸量Q=15 L/min、Q=30 L/min和Q=60 L/min時的呼吸量，并轉(zhuǎn)化為隨時間瞬時變化的鼻腔入口呼吸速率：

V(t)=Q/A(sin(2πt/T)).

(6)

式(6)中A為鼻孔面積，T為呼吸周期.

2 結(jié)果與討論

模型截面的選取自鼻腔至?xí)捁步厝×?0個特征面(圖2)：截面A-A為靠近鼻腔入口處；截面B-B和C-C為鼻腔中部，包含左右兩側(cè)的中鼻道和下鼻道；截面D-D為左右鼻道匯合處；截面E-E為鼻咽下界；截面F-F為腭咽下界；截面G-G為舌咽上界；截面H-H為舌咽下界.截面I-I為會厭處；截面J-J為下界.另外，在YZ軸組成的平面上選取了呼吸道正中間部分作為特征面.

圖2 CFD模型選取的截面Figure 2 Cross sections of CFD model

2.1 上呼吸道特征截面氣流運動特性

計算了人靜坐(Q=15 L/min)、輕微運動(Q=30 L/min)和劇烈運動(Q=60 L/min)呼吸量時的吸氣氣流場.圖3為該時刻正中氣流速度等值面圖.由于咽部非常狹窄，下游產(chǎn)生了湍流噴射現(xiàn)象，同一特征面內(nèi)氣流速度分布非常不均勻；而且可以看出，三種呼吸量的最大速度氣流，都在靠近呼吸道后壁面附近.

圖3 正中氣流速度等值面圖Figure 3 Middle airflow velocity contours of each respiratory intensity

圖4為人靜坐時吸氣中間時刻平面A-A到J-J的氣流速度等值面圖.可以看到，隨著呼吸道結(jié)構(gòu)的變化，呼吸道內(nèi)的氣流經(jīng)歷了從層流到湍流的變化過程.由于咽部和喉部的截面相較于其他部位變小了很多，致使這些部位的氣流速度較大，從而令下游各部分區(qū)域氣流速度變化劇烈；在F-F截面氣流速度達到最大值約5.3 m/s，且在舌咽上界(圖2G-G平面)靠近后壁面處，氣流逐漸產(chǎn)生了流動分離，形成了環(huán)流.

2.2 上呼吸道內(nèi)部顆粒沉積特性

圖5顯示的是呼吸量為15 L/min時粒徑分別為2.5 μm、5 μm、10 μm的可吸入顆粒在模型中的沉積分布，從圖中可以看出，顆粒主要沉積在鼻腔以及咽部和喉部的后壁面，并且隨著粒徑的增大，鼻腔顆粒沉積數(shù)變多而整個咽部沉積數(shù)減少.與圖5比較得知，湍流強度大的區(qū)域，顆粒沉積數(shù)目非常少.

圖6顯示的是人靜坐時，三種粒徑d=2.5 μm、5 μm和10 μm的顆粒，在上呼吸道鼻腔、鼻咽、腭咽、舌咽和喉部五個部位的沉積分?jǐn)?shù)曲線.從圖中可以看出，可吸入顆粒物在鼻腔和喉部沉積較多，鼻咽、腭咽和舌咽沉積相對較少，2.5 μm和5 μm的可吸入顆粒在鼻咽、腭咽、舌咽和喉部的沉積率幾乎一樣，只有鼻腔的沉積率有所差別.而10 μm的顆粒除在鼻咽的沉積率與前二者相近外，其他四個部位都有所差別，10 μm顆粒在腭咽和舌咽的沉積率小于2.5 μm和5 μm的，鼻腔和喉部則明顯大于它們二者.

圖4 各特征平面氣流速度等值面圖Figure 4 Airflow velocity contours of each cross section

圖5 Q=15 L/min顆粒物粒徑d=2.5 μm、 5 μm、10 μm沉積分布圖Figure 5 Particle deposition distribution when d=2.5 μm,5 μm,10 μm and Q=15 L/min

圖6 Q=15 L/min三種粒徑在不同部位沉積率曲線Figure 6 Deposition rate curve of each diameters at different position when Q=15 L/min

圖7 三種粒徑的總沉積率曲線Figure 7 Total deposition rate of each diameters

圖7是三種粒徑d=2.5 μm、5 μm和10 μm在上呼吸道的總沉積率隨呼吸強度變化的情況.從圖中可以看出，可吸入顆粒物在人體上呼吸道內(nèi)的沉積率隨著呼吸強度的增大而變大，這主要是由于隨著呼吸量的增加引起了湍流強度的增加，造成了顆粒湍流擴散的增加，從而造成了顆粒沉積的增加.其中，在人靜坐和輕微運動時，2.5 μm和5 μm顆粒的沉積率增加得最快，之后隨著呼吸量的增大，它們的沉積率趨于平緩；而10 μm顆粒的沉積率一直處于較高水平，在人劇烈運動時，這些大顆粒幾乎全部沉積在了上呼吸道內(nèi).

3 結(jié) 語

本文使用以醫(yī)學(xué)CT掃描為基礎(chǔ)得到的真實人體上呼吸道模型，采用大渦模擬的方法對靜坐、輕微運動和劇烈運動三種呼吸強度下的氣流運動特性進行了仿真模擬，得到了上呼吸道內(nèi)部的氣流分布情況.在吸氣中間時刻，氣流在呼吸道內(nèi)經(jīng)歷了從層流到湍流的變化過程；咽部的縮小結(jié)構(gòu)，使這一區(qū)域氣流速度猛烈上升，加劇了后面部位氣流的紊亂程度.跟蹤計算了可吸入顆粒物的沉積情況，直觀地顯示了顆粒的局部沉積位置.顆粒的沉積特性表明，大部分的顆粒沉積在鼻腔和喉部；粒徑小的顆粒與粒徑大的顆粒沉積位置有所不同；粒徑大小對總沉積率有很大影響，小顆粒有一部分會進入下呼吸道對人體健康造成進一步傷害，而大顆粒幾乎全部沉積在了上呼吸道.本文對上呼吸道氣流流場和顆粒沉積的分析，可以作為研究呼吸系統(tǒng)的致病機理和研制開發(fā)治療人體呼吸道疾病的噴霧器的參考.

[1] 劉曉月.我國霧霾天氣現(xiàn)狀與綜合治理分析[J].科技風(fēng),2014,13:275-275. LIU Xiaoyue. Present situation and comprehensive analysis of haze-fog in China[J]. Technology Wind,2014,13:275-275.

[2] 張月琴,何艾濃,閆春良.可吸入顆粒物對呼吸系統(tǒng)疾病的影響[J].職業(yè)與健康,2014,30(24):3583-3585. ZHANG Yueqin, HE Ainong, YAN Chunliang. Effect of inhaled particulate matter on respiratory system disease[J]. Occup and Health,2014,30(24):3583-3585.

[3] CHENG Y S, ZHOU Y, CHEN B T. Particle deposition in a cast of human oral airways[J]. Aerosol Science and Technology,1999,31:286-300.

[4] GRGIC B, FINLAY W H, HEENAN A F. Regional aerosol deposition and flow measurements in an idealized mouth and throat[J]. Journal of Aerosol Science,2004,35:21-32.

[5] STAPLETON K W, G E, HOSKINSON M K. On the suitability of κ-ε turbulence modeling for aerosol deposition in the mouth and throat: A comparison with experiment[J]. Journal of Aerosol Science,2000,31:739-749.

[6] ZHANG Z, KLEINSTREUER C, DONOHUE J F. Comparison of micro- and nano-size particle depositions in a human upper airway model[J]. Journal of Aerosol Science,2005,36:211-233.

[7] 文謀,胡徐趣,盧志明.顆粒物在人體上呼吸道模型中沉降規(guī)律的數(shù)值研究[J].力學(xué)季刊,2011,32(4):515-524. WEN Mou, HU Xuqu, LU Zhiming.Numerical study on the deposition characteristics of particles in a human upper respiratory tract model[J]. Chinese Quarterly of Mechanics,2011,32(4):515-524.

[8] 趙秀國,徐新喜,孫棟,等.人體上呼吸道內(nèi)氣溶膠沉積的數(shù)值仿真研究[J].醫(yī)用生物力學(xué),2012,27(1):90-95. ZHAO Xiuguo, XU Xinxi, SUN Dong, et al. Numerical simulation of aerosol deposition in human upper respiratory tract[J]. Journal of Medical Biomechanics,2012,27(1):90-95.

[9] 徐新喜,趙秀國,譚樹林,等.人體上呼吸道內(nèi)氣流運動特性的數(shù)值模擬分析[J].計算力學(xué)學(xué)報,2010,27(5):881-886. XU Xinxi, ZHAO Xiuguo, TAN Shulin, etal. Numerical simulation of airflow characteristics in human upper respiratory tract[J]. Chinese Journal of Computational Mechanics,2010,27(5):881-886.

[10] 林江,胡桂林.人體呼吸時微細(xì)顆粒在呼吸道中運動特性的研究[J].工程熱物理學(xué)報,2006,27(1):225-228. LIN Jiang, HU Guilin. Research on micro-particle motion characteristics in the human airway during respiration[J]. Journal of Engineering Thermophysics,2006,27(1):225-228.

[11] 蔡志良,孫在,陳秋方,等.廚房油煙超細(xì)顆粒排放特征[J].中國計量學(xué)院學(xué)報,2015,26(1):75-79. CAI Zhiliang, SUN Zai, CHEN Qiufang, etal.The emission characteristics of ultrafine particles produced by cooking[J]. Journal of China University of Metrology,2015,26(1):75-79.

Simulation of airflow and inhalable particles deposition in human upper respiratory tract

YAN Weiwei, SUN Yusheng, LIU Xingli

(College of Metrology and Measurement Engineering, China Jiliang University, Hangzhou 310018, China)

The airflow in human upper respiratory tract was numerically studied by using the large eddy simulation method. Airflow velocity and properties of vortexes were studied. The motion and deposition of particles were investigated with Lagrangian method, and the deposition positions of the particles were displayed. The results show that vortexes and backflows occur mostly at glossopharyngeum and the velocity at the narrow part is the largest with flow-jet. Most of the particles deposit at the nasal cavity and the laryngopharynx. The deposition rate is influenced heavily by the particle diameter.

upper respiratory traet; large eddy simulation; flow characteristics; particle deposition.

2096-2835(2016)04-0372-05

10.3969/j.issn.2096-2835.2016.04.003

2016-07-16 《中國計量大學(xué)學(xué)報》網(wǎng)址：zgjl.cbpt.cnki.net

國家自然科學(xué)基金資助項目(No.11202203)，浙江省自然科學(xué)基金資助項目(No.LY15A020005).

TP391.9

A