上呼吸道流场及颗粒物沉积规律的模拟

严微微，孙玉胜，刘星利

(中国计量大学 计量测试工程学院，浙江 杭州 310018)

上呼吸道流场及颗粒物沉积规律的模拟

严微微，孙玉胜，刘星利

(中国计量大学 计量测试工程学院，浙江 杭州 310018)

基于医学CT扫描图像得到的真实人体上呼吸道模型，采用大涡模拟方法计算了上呼吸道内部的气流流动，分析了气流流速和涡旋的特征.在拉格朗日框架下计算了可吸入颗粒物的沉积，对颗粒的局部沉积位置进行了直观显示.研究结果表明，舌咽部位涡旋和回流最多，狭窄结构处流速最大，且产生了气流喷射现象；大部分的颗粒沉积在鼻腔和喉部，颗粒粒径对沉积率有很大影响.

上呼吸道；大涡模拟；气流特征；颗粒沉积

雾霾天气对人们的日常生活影响恶劣[1]，可吸入颗粒物会对人们的呼吸系统产生严重影响，导致上呼吸道感染、慢性阻塞性肺病、支气管哮喘等疾病的发病率明显升高[2].深入了解人体上呼吸道内部的气流特征及可吸入颗粒物在其中的沉积规律，不仅可以帮助判断病人的病情，还可以为呼吸系统疾病的致病机理提供可靠的依据.

目前，已有一些关于呼吸道内气流流场和颗粒运动的计算模拟研究，实验方面：Cheng等根据某志愿者的呼吸道尺寸建立了上呼吸道三维模型[3]，并选用9种不同粒径的颗粒进行了呼吸实验，发现颗粒大小与呼吸强度对颗粒的沉积率有很大影响.Grgic等运用CT图像生成技术，对人体呼吸道进行数据采集[4]，建立了三维呼吸道模型，实验统计了颗粒在口腔、咽部、喉部和气管的沉积百分比，发现颗粒大量沉积在喉部与气管处，颗粒惯性是影响颗粒沉积率的关键因素.数值模拟方面：Stapleton等基于Grgic的呼吸道模型[5]，采用k-ε湍流模型模拟上呼吸道中的气流，并将Lagrange跟踪模型应用于颗粒运动，结果显示湍流模型能够很好地预测低呼吸强度时上呼吸道中的气流运动.Zhang等基于Cheng的呼吸道模型[6]，采用k-ε湍流模型计算流场，分别用欧拉-欧拉法和欧拉-拉格朗日方法模拟追踪了呼吸道内颗粒的运动情况，在不同入口气流边界条件下分别比较了从口腔到前三级支气管内，微颗粒和纳米颗粒的沉积模式和局部浓度.纵观国内文献，他们的几何模型大多参照了Stapleton等的简化模型[7-8]；所使用的计算模型多为低雷诺数k-ω湍流模型[9]和RNGk-ε模型[10]，发现颗粒沉积特性与其质量相关[11].本文在此基础上，基于医学CT扫描的原始图像，重构得到真实人体上呼吸道的三维几何模型；并且采用更为精确的可得到瞬时流场特性的大涡模拟方法求解气流运动特性，利用与实际情况相符合的拟合函数仿真随时间变化的呼吸强度；而且还在拉格朗日框架下跟踪颗粒物的运动情况，讨论了颗粒物在人体呼吸道内部的沉积规律.

1 模型描述

1.1 几何模型

原始图像数据是由CT扫描一位中国男性患者的胸部得到的.图像由分辨率为0.7×0.7 mm2的轴向平面制成，层厚为0.625 mm.上呼吸道模型三维点云数据由图像处理软件Mimics重建，其几何模型如图1.

图1 上呼吸道三维结构模型Figure 1 3D model of upper respiratory tract

对此模型可以进行网格划分和计算流体力学分析.其中，鼻腔入口及喉部的延长段是为了更好地处理边界条件，方便计算.从图中可以看出，该呼吸道几何结较为复杂，鼻腔和咽喉部位突缩结构明显，因此与之前文献中所使用的简化模型相比，此物理模型更具准确性.

1.2 控制方程

由于呼吸道内结构复杂，存在急剧收缩和扩展的面，本文采用大涡模拟的方法描述气流场，用以更加准确地得到内部层流和复杂湍流的信息.在该方法中，对变量x进行滤波后的方程为

(1)

其中，V是计算单元的体积，滤波方程定义为

(2)

这一滤波过程有效地过滤了小于过滤宽度或网格间距的涡旋，因此经过滤波处理的Navier-Stokes方程为：

连续性方程

(3)

动量方程

对于颗粒相，认定其为稀相，且其密度远大于气相，不考虑颗粒之间的相互作用，只考虑颗粒受到的重力及Stokes力和Saffman升力，因此单颗粒的运动方程可以表达为

(5)

式中mp是颗粒质量，vi是颗粒的速度，dp是颗粒的直径，CD是颗粒的阻力系数，g是重力加速度.右边三项依次是单颗粒受到的Stokes阻力、重力和Saffman升力.

1.3 数值仿真方法

在进行数值仿真计算时，假设呼吸道壁面为刚性光滑面，计算域为鼻腔入口到会厌底部，计算软件采用商业软件Fluent，使用非结构化网格，共划分网格120万左右，在颚咽即结构最狭窄处进行了网格加密处理.使用大涡模拟(Large Eddy Simulation, LES)湍流模型，速度压力耦合的求解采用PISO算法，瞬态方程采用二阶隐式的方法求解.鼻腔入口采用速度入口边界条件，出口采用压力出口边界条件，压力为零.入口气流流量随时间周期变化，计算时间步长0.001 s，自行编制C语言函数模拟呼吸量Q=15 L/min、Q=30 L/min和Q=60 L/min时的呼吸量，并转化为随时间瞬时变化的鼻腔入口呼吸速率：

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

(6)

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

2 结果与讨论

模型截面的选取自鼻腔至会厌共截取了10个特征面(图2)：截面A-A为靠近鼻腔入口处；截面B-B和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为该时刻正中气流速度等值面图.由于咽部非常狭窄，下游产生了湍流喷射现象，同一特征面内气流速度分布非常不均匀；而且可以看出，三种呼吸量的最大速度气流，都在靠近呼吸道后壁面附近.

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

图4为人静坐时吸气中间时刻平面A-A到J-J的气流速度等值面图.可以看到，随着呼吸道结构的变化，呼吸道内的气流经历了从层流到湍流的变化过程.由于咽部和喉部的截面相较于其他部位变小了很多，致使这些部位的气流速度较大，从而令下游各部分区域气流速度变化剧烈；在F-F截面气流速度达到最大值约5.3 m/s，且在舌咽上界(图2G-G平面)靠近后壁面处，气流逐渐产生了流动分离，形成了环流.

2.2 上呼吸道内部颗粒沉积特性

图5显示的是呼吸量为15 L/min时粒径分别为2.5 μm、5 μm、10 μm的可吸入颗粒在模型中的沉积分布，从图中可以看出，颗粒主要沉积在鼻腔以及咽部和喉部的后壁面，并且随着粒径的增大，鼻腔颗粒沉积数变多而整个咽部沉积数减少.与图5比较得知，湍流强度大的区域，颗粒沉积数目非常少.

图6显示的是人静坐时，三种粒径d=2.5 μm、5 μm和10 μm的颗粒，在上呼吸道鼻腔、鼻咽、腭咽、舌咽和喉部五个部位的沉积分数曲线.从图中可以看出，可吸入颗粒物在鼻腔和喉部沉积较多，鼻咽、腭咽和舌咽沉积相对较少，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在上呼吸道的总沉积率随呼吸强度变化的情况.从图中可以看出，可吸入颗粒物在人体上呼吸道内的沉积率随着呼吸强度的增大而变大，这主要是由于随着呼吸量的增加引起了湍流强度的增加，造成了颗粒湍流扩散的增加，从而造成了颗粒沉积的增加.其中，在人静坐和轻微运动时，2.5 μm和5 μm颗粒的沉积率增加得最快，之后随着呼吸量的增大，它们的沉积率趋于平缓；而10 μm颗粒的沉积率一直处于较高水平，在人剧烈运动时，这些大颗粒几乎全部沉积在了上呼吸道内.

3 结 语

本文使用以医学CT扫描为基础得到的真实人体上呼吸道模型，采用大涡模拟的方法对静坐、轻微运动和剧烈运动三种呼吸强度下的气流运动特性进行了仿真模拟，得到了上呼吸道内部的气流分布情况.在吸气中间时刻，气流在呼吸道内经历了从层流到湍流的变化过程；咽部的缩小结构，使这一区域气流速度猛烈上升，加剧了后面部位气流的紊乱程度.跟踪计算了可吸入颗粒物的沉积情况，直观地显示了颗粒的局部沉积位置.颗粒的沉积特性表明，大部分的颗粒沉积在鼻腔和喉部；粒径小的颗粒与粒径大的颗粒沉积位置有所不同；粒径大小对总沉积率有很大影响，小颗粒有一部分会进入下呼吸道对人体健康造成进一步伤害，而大颗粒几乎全部沉积在了上呼吸道.本文对上呼吸道气流流场和颗粒沉积的分析，可以作为研究呼吸系统的致病机理和研制开发治疗人体呼吸道疾病的喷雾器的参考.

[1] 刘晓月.我国雾霾天气现状与综合治理分析[J].科技风,2014,13:275-275. LIU Xiaoyue. Present situation and comprehensive analysis of haze-fog in China[J]. Technology Wind,2014,13:275-275.

[2] 张月琴,何艾浓,闫春良.可吸入颗粒物对呼吸系统疾病的影响[J].职业与健康,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] 文谋,胡徐趣,卢志明.颗粒物在人体上呼吸道模型中沉降规律的数值研究[J].力学季刊,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] 赵秀国,徐新喜,孙栋,等.人体上呼吸道内气溶胶沉积的数值仿真研究[J].医用生物力学,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] 徐新喜,赵秀国,谭树林,等.人体上呼吸道内气流运动特性的数值模拟分析[J].计算力学学报,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] 林江,胡桂林.人体呼吸时微细颗粒在呼吸道中运动特性的研究[J].工程热物理学报,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] 蔡志良,孙在,陈秋方,等.厨房油烟超细颗粒排放特征[J].中国计量学院学报,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 《中国计量大学学报》网址：zgjl.cbpt.cnki.net

国家自然科学基金资助项目(No.11202203)，浙江省自然科学基金资助项目(No.LY15A020005).

TP391.9

A