研究颗粒物在人体肺腺泡区的沉积规律对了解肺部疾病致病机制、促进肺部疾病的预防和治疗都有重要意义。现行研究大多关注颗粒物在人体肺腺泡区的最终沉积率，而对其沉积动态过程鲜有涉及。本文建立了 G3～G7 共五个多肺泡模型，引入颗粒物沉积速度评价参数 1/4 沉积时间，通过数值模拟研究了模型级数和结构、颗粒物粒径、呼吸模式等因素对 0.1～5 μm 粒径范围内颗粒物沉积规律的影响，归纳总结了不同粒径颗粒物的沉积特点，为进一步了解人体肺腺泡区颗粒物的运动规律提供了新视角。结果表明：模型级数和结构是影响各粒径颗粒物沉积的重要因素。0.1 μm 颗粒物沉积受布朗力主导，沉积率高，沉积速度快，沉积曲线呈对数型分布；5 μm 颗粒物沉积受重力主导，沉积率高，沉积较快，沉积曲线呈“S”型分布；0.3～1 μm 颗粒物沉积则受惯性冲击影响较大，随着呼吸模式的改变沉积规律变化明显。本文的研究方法和结果能为进一步探究肺部疾病致病机制和肺部疾病的预防与治疗提供理论依据和数据支持。
引用本文： 李鹏辉, 徐新喜, 李蓉, 乔扬. 多因素影响下的人体肺腺泡区吸入颗粒物沉积规律数值模拟研究. 生物医学工程学杂志, 2020, 37(5): 793-801. doi: 10.7507/1001-5515.201909004 复制
引言
人体肺腺泡是指终末细支气管以下的肺组织，包括呼吸性细支气管、肺泡管、肺泡囊和肺泡，是呼吸作用的主要区域。细颗粒物的吸入和沉积会影响肺部健康，引发肺部疾病，因此研究细颗粒物在人体肺腺泡区的沉积规律对肺部疾病致病机制的探究以及肺部疾病的预防和治疗都有重要意义。
人体肺腺泡区结构复杂、尺寸微小，尚未获得真实的计算机模型。现行数值模拟研究多采用简化的计算机模型，主要分为单肺泡模型和多肺泡模型。单肺泡研究仅对单个肺泡建模，计算量较小，着重关注单个肺泡的内部特性。Sznitman 等[1]对处于不同级数的单肺泡进行了数值模拟研究，发现肺泡内流场形态随着级数增加由循环状向辐射状变化。单肺泡研究只能关注单个肺泡内流场形态和颗粒物沉积状况，无法探究肺腺泡区的整体特性，因此国内外学者纷纷建立多肺泡模型进行进一步研究。目前，多肺泡模型主要分为四类[2]：“柱状气管+球状肺泡”的多肺泡模型、基于十四面体的多肺泡模型、泰森多边形法构建的多肺泡模型、计算机成像和重建技术构建的多肺泡模型。Ma 等[3]建立了“柱状气管+球状肺泡”的 5 级多肺泡模型，发现呼吸周期、潮气量和重力等对颗粒物沉积有较大影响。Longest 等[4]建立了 3 个基于十四面体的多肺泡模型，发现屏气能大幅度增加颗粒物的沉积率。Hofemeier 等[5]采用泰森多边形法构建了人体 1/8（第 3～8 级）肺腺泡模型，研究了 0.005～5.0 μm 粒径范围内不同大小颗粒物的沉积状况，发现颗粒物沉积率随粒径增加呈先下降后上升的趋势。Sera 等[6]采用微计算机断层扫描（micro focal computed tomography system，micro-CT）技术构建了大鼠肺腺泡模型，发现 0.8～1.0 μm 粒径颗粒物多沉积于肺腺泡入口处，而 0.1～0.2 μm 粒径颗粒物则主要沉积于模型中部和末端。
上述数值模拟研究大多关注颗粒物在人体肺腺泡区局部或整体的最终沉积率，而对其沉积动态过程鲜有涉及，无法进一步了解人体肺腺泡区颗粒物的运动规律。本文建立了五个“柱状气管+球状肺泡”多肺泡模型，通过数值模拟研究了模型级数和结构、颗粒物粒径、呼吸模式等因素对 0.1～5 μm 粒径范围内颗粒物沉积规律的影响，总结了颗粒物沉积率和沉积速度的变化规律，归纳了沉积曲线特征，以期为进一步研究肺部疾病致病机制和相关疾病的预防与治疗提供理论依据和数据支持。
1 模型与方法
1.1 人体肺腺泡模型
在经典的全肺模型（Weibel A 模型[7]）中，肺部分为 0～23 级。其中，17～23 级部分含有肺泡，是呼吸作用的主要区域，为肺腺泡区。人体肺腺泡区尺寸数据如表 1 所示。
表1
人体肺腺泡区尺寸数据
Table1.
Geometric data of human pulmonary acinus
由表 1 数据可知，肺腺泡区尺寸数据随级数变化显著。从 17 级至 23 级，肺腺管长度、管径数据逐级递减，最大缩减率分别为 64.54%、24.07%；肺泡直径和单节肺腺管肺泡数量逐级递增，最大增长率分别为 79.25%、300%。Sznitman 等[1]研究发现模型尺寸是影响数值模拟结果的重要因素。十四面体多肺泡模型中肺泡形状、大小一致[8]，仅适用于肺腺泡区局部特性探究或尺寸影响可忽略不计的情况，因而本文采用“柱状气管+球状肺泡”架构，将肺腺管简化为圆柱状，肺泡简化为附着在肺腺管表面的球，建立了如图 1 所示的多肺泡模型。
图1
“柱状气管+球状肺泡”架构多肺泡模型
Figure1.
Multi-alveolar model of “columnar trachea + globular alveolar” structure
然后，在该模型的基础上依据表 1 尺寸数据建立了图 2 所示的 G3（17～19 级）、G4（17～20 级）、G5（17～21 级）、G6（17～22 级）、G7（17～23 级）共 5 个多肺泡模型。各模型所含级数、肺泡数、长度、体积等参数如表 2 所示。
图2
多肺泡模型
Figure2.
Multi-alveolar models
表2
多肺泡模型参数
Table2.
Multi-alveolar models parameters
1.2 数值模拟方法
1.2.1 数值模拟软件
采用 Inventor 2015 进行多肺泡模型建模，采用 ICEM CFD 15.0 进行模型网格划分，采用 FLUENT 15.0 进行数值模拟计算。
1.2.2 网格划分结果
采用非结构网格对 G3～G7 模型进行网格划分，各模型网格总数及最小网格质量如表 3 所示，最小网格质量均大于 0.3，可以用于数值模拟计算。
表3
G3～G7 模型网格划分结果
Table3.
Meshing results of G3―G7 models
1.2.3 控制方程
① 流体控制方程：人体肺腺泡区气流运动为等温、不可压缩定常流动。满足质量、动量守恒方程。
质量守恒方程：
 |
式中，
表示流体密度，
= 1.185 kg/m3，
为流体速度，t为时间。故上式可表示为：
 |
动量守恒方程：
 |
式中，P为流体压力，
= 1.785e-5 表示运动粘度，
为流体体积力矢量。
② 颗粒物控制方程：本研究中，忽略 Magnus 力、Saffman 力、Basset 力等对颗粒物运动影响较小的力和颗粒物之间的相互影响，仅考虑重力
、布朗力
和流体曳力
的作用效果，因此单个颗粒的运动控制方程为：
 |
 |
 |
 |
其中，其中 
和 
分别为颗粒物的密度和直径，
为随机数，
为玻尔兹曼常数，
为空气热力学温度，
为颗粒半径，
为计算时间间隔。
1.2.4 边界条件设置
如图 3 所示，以 G3 模型区域划分为例，本研究中将模型分为两部分：顶端平面（黑色阴影部分）和其他部分。其中，顶端平面为气流和颗粒物入口，采用入口速度函数实现气流和颗粒物的进入和排出模型，以此模拟呼吸作用下人体肺腺泡区的气流和颗粒物的运动。其他部分为无滑移壁面，颗粒物在运动过程中接触壁面会被捕捉，即为沉积。
图3
G3 模型区域划分
Figure3.
Parts of G3 model
1.2.5 入口速度函数
Oakes 等[9]测量了 21 岁健康男性休息状态下的呼吸流量曲线，近似为正弦曲线，因此本研究中入口速度采用随时间变化的正弦函数，速度正方向为 − Z，入口速度 
可表示为:
 |
其中，
为入口气流最大速度，
为呼吸周期。
则数值模拟一周期模型内气体最大进入量
为：
 |
积分可得:
 |
其中，
为入口管径，即 17 级肺腺管管径。
以上述建立的多肺泡模型代表人体肺腺泡进行数值模拟研究，模型体积
对应功能残气量（functional residual capacity，FRC），最大进气量
对应潮气量
，则有：
 |
根据 Haefeli-Bleuer 等[10]的测量结果，取人体正常呼吸状态下功能残气量为 3 L，潮气量为 600 mL，呼吸周期为 4 s。依据式（1）（2）（3）计算可得，正常呼吸状态下（
= 600 mL，T = 4 s），G3～G7 模型入口气流最大速度
和入口速度如表 4 所示。
表4
正常呼吸状态下 G3～G7 模型入口气流最大速度和入口 速度
Table4.
The maximum velocity and inlet velocity of inlet airflow in G3−G7 models under normal respiratory mode
1.2.6 呼吸模式
当人体处于不同呼吸状态时，潮气量和呼吸周期存在明显差异。由式（1）（2）（3）可知，潮气量、呼吸周期能直接决定模型入口速度函数，进而影响颗粒物在人体肺腺泡区的沉积规律。从肺腺泡形态结构来说，G7 模型包含级数更多，肺泡数量也更多，更加符合人体实际，因此本研究基于 G7 模型进一步研究呼吸模式对人体肺腺泡区颗粒物运动规律的影响。取呼吸周期为 4 s，潮气量分别为 300、600、900 mL，代表人体处于浅呼吸、正常呼吸和深呼吸状态。取潮气量为 600 mL，呼吸周期分别为 2、4、6 s，代表人体处于急促呼吸、正常呼吸和缓慢呼吸状态，由式（1）（2）（3）计算可得，不同呼吸模式下 G7 模型入口速度如表 5 所示。
表5
不同呼吸模式下 G7 模型入口速度
Table5.
The velocity of inlet airflow in G7 model under diffe rent respiratory mode
1.2.7 颗粒物设置
数值模拟计算中，将颗粒物简化为球状，密度为 1 kg/m3，选取直径为 0.1、0.3、0.5、0.8、1、3、5 μm 等能够进入人体肺腺泡区的颗粒物[11]作为研究对象。采用 FLUENT DPM（Discrete Phase Model）模型，激活重力、布朗力、曳力选项。数值模拟开始时，在模型入口平面基于网格释放颗粒物，使颗粒物均匀地分布在气流中并随气流进入模型。利用 DPM 模型中颗粒物沉积数量统计功能，每 
时间间隔记录一次颗粒物沉积数量，计算总时长为一个周期。定义 
为t时刻颗粒物沉积率，
为 t 时刻颗粒物沉积数量，
为释放的颗粒物总量，则：
 |
1.2.8 颗粒物沉积速度评价参数
通常而言，颗粒物沉积速度评价参数是首次沉积发生的时间，即颗粒物首次沉积发生的时间越短，沉积速度越快。但实际上，研究中多采用间隔采样法对颗粒物沉积率进行记录，无法得到颗粒物首次沉积发生的具体时间。而且，颗粒物沉积为非线性，颗粒物首次沉积时间并不能真实反映该粒径颗粒物沉积的快慢。为了更好地评价颗粒物在人体肺腺泡区沉积速度差异，引入颗粒物沉积速度评价参数：1/4 沉积时间，即某粒径颗粒物沉积率达到 1/4 最终沉积率的时间，记为 
。本研究中数值模拟时间为一周期，T 时刻颗粒物沉积率即为最终沉积率 
，则 1/4 最终沉积率为 
。假定 
、
（
）为两相邻记录时刻，即 
，且这两时刻沉积率满足：
 |
那么，1/4 沉积时间位于 
～
之间，假定 
～
内颗粒物沉积速度恒定，为线性沉积，则可利用插值计算估算得：
 |
化简得：
 |
1/4 沉积时间以颗粒物达到 1/4 最终沉积率的时间作为沉积快慢评价参数，能对颗粒物沉积快慢进行较好地评价。
2 结果和讨论
2.1 模型对人体肺腺泡区吸入颗粒物沉积特性的影响
2.1.1 模型级数影响
模型的建立是数值模拟研究的前提和基础，也是影响数值模拟结果的重要因素。由表 2 可知，G3～G7 模型包含肺腺泡级数不同，模型总长度和体积也存在差异，因此研究 G3～G7 模型级数差异对颗粒物沉积特性的影响十分必要。正常呼吸模式下，G3～G7 模型中各粒径颗粒物最终沉积率如图 4 所示。由图 4 可知，由 G3 模型至 G7 模型，随着模型级数的增加，0.1 μm 粒径颗粒物最终沉积率呈微弱上升趋势，整体差异不大；0.3、0.5、0.8、1 μm 粒径颗粒物最终沉积率呈明显上升趋势；3、5 μm 粒径颗粒物最终沉积率则呈现出下降趋势。
图4
G3～G7 模型中各粒径颗粒物最终沉积率
Figure4.
Final deposition fraction of particles in G3―G7 models
研究表明，颗粒物在人体肺腺泡区沉积机制主要为布朗扩散、重力沉降和惯性冲击，而惯性冲击是影响 0.3～1 μm 粒径颗粒物沉积的重要因素[12]。由表 4 可知，正常呼吸状态下，模型所含级数越多，体积越大，则模型入口速度越大，惯性冲击对颗粒物沉积的影响也相应越大。因此随着模型级数的增多，0.3～1 μm 粒径颗粒物最终沉积率呈现出明显上升趋势。0.1 μm 粒径颗粒物沉积由布朗扩散主导，惯性冲击影响不大，因此沉积率仅微弱上升。由表 2 可知，模型级数的增加不仅体现为体积的增加，还体现在模型总长度的增长。3、5 μm 粒径颗粒物沉积由重力沉降主导，布朗扩散和重力沉积影响较小，随着模型总长度的增加，颗粒物沿重力方向沉积所需运动路程相应增加，沉积所需运动时间更长，最终沉积率也相应更低，因而最终沉积率随模型级数增加呈现出明显的下降趋势。Katan 等[13]在研究年龄对肺腺泡区颗粒物沉积状态的影响时发现，成人肺腺泡模型尺寸较大，颗粒物沉积路径增加，不利于颗粒物的沉积，因而沉积率较低，与本文研究结果一致。
2.1.2 模型结构影响
由图 4 可知，上一小节关于模型级数对颗粒物沉积影响规律的描述，只是整体趋势的归纳和总结。G6 模型中，0.1、0.3、0.5 μm 粒径颗粒物沉积率较 G5 模型中沉积率稍低，分别相差 2.23%、2.13%、2.90%，鉴于沉积率差值较小，可以认为基本不影响上一小节趋势的归纳。而 G3 模型中，3 μm 粒径颗粒物的最终沉积率则与理论值相差甚远，仅为 10.53%（理论值应大于 G4 模型中 3 μm 粒径颗粒物最终沉积率，即 > 46.49%），这是模型几何结构对颗粒物沉积影响的典型体现。如图 5 所示为
= 2 s 时刻 G3、G4 模型中 3 μm 粒径颗粒物的分布情况。正常呼吸模式下，2 s 时刻颗粒物运动路径最长，最容易接触壁面并沉积。而 G3 模型中，2 s 时刻颗粒物距接触壁面仍有一定距离，因此沉积率较小。G4 模型中，入口速度增大，运动路径增加，能够接触壁面并沉积，因而沉积率较高。同理，缩减模型中 17 级肺腺管长度（
）或改变重力方向[14]也能大幅度提升 G3 模型中 3 μm 颗粒物的沉积率。
图5
t = 2 s 时刻 G3、G4 模型中 3 μm 粒径颗粒物分布情况
Figure5.
Distribution of 3 μm particles in G3 and G4 model at t = 2 s
2.2 颗粒物粒径对人体肺腺泡区吸入颗粒物沉积特性的影响
正常呼吸模式下，G3～G7 模型中颗粒物最终沉积率随粒径变化趋势如图 6 所示，不同模型间各粒径颗粒物沉积曲线略有差异，但基本形状均为“U”型，即在颗粒物粒径为 0.1 μm 时，颗粒物沉积率较高，随着颗粒物粒径的增加，沉积率先下降、后上升，并在颗粒物粒径为 1 μm 左右时沉积率最低，呈现出两端高、中间低的趋势，这一趋势与文献数据[5, 15-16]保持一致。这是因为，布朗力和重力是影响颗粒物沉积的主要因素。布朗力对微小颗粒物（< 0.1 μm）沉积起主导作用，重力对大颗粒物（> 2 μm）沉积起主导作用。因此，当颗粒物粒径较小时，沉积受布朗力主导，沉积率较高。颗粒物粒径较大时，沉积受重力主导，沉积率较高。随着颗粒物粒径增加，布朗力作用减弱、重力作用加强，而 1 μm 左右粒径颗粒物受两者影响均较小，因而沉积率最低，故沉积曲线呈现为“U”型。
图6
G3～G7 模型中颗粒物最终沉积率随粒径变化趋势
Figure6.
Final deposition fraction of particles as a function of diameter in G3—G7 models
2.3 人体肺腺泡区吸入颗粒物动态沉积特性
2.1 小节探讨了模型级数对人体肺腺泡区颗粒物沉积特性的影响，从肺腺泡形态结构来说，G7 模型包含级数更多，肺泡数量也更多，更加符合人体实际，因此本研究基于 G7 模型进一步研究人体肺腺泡区吸入颗粒物的动态沉积特性和呼吸模式对颗粒物运动规律的影响。正常呼吸模式下（
= 600 mL，T = 4 s），G7 模型中各粒径颗粒物沉积率随时间变化曲线如图 7 所示，G7 模型中各粒径颗粒物沉积率均随时间增加，至 4.0 s 时刻达到最大值，即最终沉积率。其中，0.1 μm 颗粒物前期沉积速度快，中后期沉积缓慢并趋于稳定，沉积率变化曲线呈对数型分布。5 μm 颗粒物前期沉积速度慢，中期沉积快速并趋于稳定，后期沉积缓慢，沉积曲线呈“S”型。0.3 μm 颗粒物沉积趋势与 0.1 μm 基本相似，沉积率相对较低，对数型沉积曲线相对较矮。3 μm 颗粒物沉积曲线与 5 μm 基本相似，沉积率较低，呈较矮的“S”型曲线。其他粒径颗粒物前期沉积缓慢，中后期沉积率缓慢增长，沉积曲线基本呈折线型。
图7
G7 模型各粒径颗粒物沉积率变化曲线
Figure7.
Variation curves of deposition fraction of particles in G7 model
由图 7 沉积曲线可以看出，0.0～1.0 s 间除 0.1 μm 颗粒物沉积率增长较快外，其他粒径颗粒物沉积率均较小，很难判断各粒径颗粒物沉积速度差异，可见引入合理参数对沉积速度进行评价的必要性。G7 模型中各粒径颗粒物 1/4 沉积时间如图 8 所示，从 1/4 沉积时间数值来看，0.1、0.3、5 μm 粒径颗粒物 1/4 沉积时间相对较短，分别为 0.60、1.03、1.33 s。由此可见，微小颗粒物沉积速度最快，大颗粒物沉积速度次之，而 0.5～3 μm 颗粒物沉积速度则相对较慢。这是因为，数值模拟初期，气流速度较小，颗粒物速度较小，惯性冲击影响较小，重力沉降则需要颗粒物沿肺腺管运动一定时间至模型分叉处才能产生沉积。因此，数值模拟初期颗粒物沉积受布朗扩散主导，微小颗粒物沉积速度最快；大颗粒物沉积受重力主导，运动一定时间接触壁面后开始沉积，沉积速度也相应较快；而其他粒径颗粒物因数值模拟初期气流速度小，惯性冲击影响作用小，沉积较慢。
图8
G7 模型中各粒径颗粒物 1/4 沉积时间
Figure8.
1/4 deposition time of particles in G7 model
2.4 呼吸模式对人体肺腺泡区吸入颗粒物沉积特性的影响
2.4.1 潮气量影响
潮气量分别为 300、600、900 mL 时，G7 模型中各粒径颗粒物最终沉积率和 1/4 沉积时间如图 9 所示。
图9
不同潮气量下颗粒物的最终沉积率和 1/4 沉积时间
Figure9.
Final deposition fraction and 1/4 deposition time of particles under different tidal volume
由图 9 可知，三种潮气量下，颗粒物最终沉积率随粒径变化趋势均呈“U”型，1/4 沉积时间的变化趋势也基本保持一致，说明在人体肺腺泡区潮气量的变化不会改变流场状态[17]，仍为充分发展的层流。同一粒径颗粒物在潮气量为 300 mL 时最终沉积率最低，1/4 沉积时间最长，潮气量为 600 mL 时最终沉积率和 1/4 沉积时间次之，潮气量为 900 mL 时最终沉积率最高，1/4 沉积时间最短。潮气量对颗粒物沉积规律的影响规律可概括为：潮气量越大，颗粒物沉积率越高，沉积速度越快。Hofemeier 等[5]分别研究了潮气量为 500 mL 和潮气量为 2 500 mL 下的颗粒物沉积状况，发现潮气量为 2 500 mL 时颗粒物在肺腺泡内总体沉积率更高，与本文研究结果保持一致。如表 5 所示，潮气量的增加会导致模型入口速度的增加，受流场曳力影响，颗粒物运动速度随之增加，惯性冲击影响作用加大，因而颗粒物沉积数增多，沉积速度加快。
值得注意的是，不同粒径颗粒物沉积规律率受潮气量影响程度不同，0.1 和 5 μm 粒径颗粒物受潮气量影响较小，沉积率极差量分别为 5.47% 和 6.28%，1/4 沉积时间极差量分别为 0.26 s 和 0.42 s。0.3～1 μm 粒径颗粒物则受潮气量影响较大，沉积率极差量分别为 12.12%、18.17%、23.25%、20.85%，1/4 沉积时间极差量分别为 0.68、0.80、0.83、0.98 s。这一结果说明惯性冲击是影响 0.3～1 μm 颗粒物沉积的关键因素，对微小颗粒物和大颗粒物的沉积则影响较小。
2.4.2 呼吸周期影响
呼吸周期分别为 2、4、6 s 时，G7 模型中各粒径颗粒物的最终沉积率和 1/4 沉积时间如图 10 所示，呼吸周期对颗粒物沉积规律的影响与潮气量影响较为一致，其影响规律可概括为：呼吸周期越长，颗粒物最终沉积率越大，沉积速度越慢，这一结果与文献[18-19]保持一致。这是因为随着呼吸周期的增加，数值模拟时间加长，颗粒物沉积更加充分，沉积率相应更高。由表 5 可知，呼吸周期的增加还会导致模型入口气流速度的减小，惯性冲击作用减小，因而颗粒物沉积速度减缓，1/4 沉积时间更长。
图10
不同呼吸周期下颗粒物的最终沉积率和 1/4 沉积时间
Figure10.
Final deposition fraction and 1/4 deposition time of particles under different respiratory cycle
3 结语
本文一共建立了 G3～G7 五个多肺泡模型，采用数值模拟方法研究了模型级数和结构、颗粒物粒径、呼吸模式等因素对 0.1～5 μm 粒径范围内 7 种颗粒物沉积率和沉积速度的影响。结果表明：模型级数和几何结构等是造成模型间最终沉积率差异的重要因素。颗粒物沉积率随粒径变化曲线呈“U”型，其中，细微颗粒物（0.1 μm）沉积由布朗扩散主导，沉积率高，沉积速度快，沉积曲线呈对数型分布；大颗粒物（5 μm）沉积由重力沉降主导，沉积率高，沉积速度较快，沉积曲线呈“S”型分布；呼吸模式对颗粒物沉积规律的影响体现为：潮气量越大，颗粒物沉积率越高，沉积越快；呼吸周期越长，颗粒物沉积率越高，沉积越慢，这一结果与文献数据保持一致。本文建立的“柱状气管+球状肺泡”多肺泡模型和研究方法能为人体肺腺泡区吸入颗粒物运动规律的进一步探究提供参考和借鉴，研究结果和数据能为肺部疾病致病机制的探究、肺部疾病的预防和治疗提供理论依据和数据支持。
本文数值模拟过程中，模型壁面刚性固定，以模型入口气流速度的正弦变化模拟呼吸过程，忽略了真实呼吸过程中肺腺泡模型边界运动对结果的影响，后期还需进一步研究模型边界运动的影响。此外，颗粒物沉积的影响因素选取还需进一步细究，考虑更多其他因素，如重力[20]、摒气时间[21]、肺部病变[22]等的影响。
利益冲突声明：本文全体作者均声明不存在利益冲突。
引言
人体肺腺泡是指终末细支气管以下的肺组织，包括呼吸性细支气管、肺泡管、肺泡囊和肺泡，是呼吸作用的主要区域。细颗粒物的吸入和沉积会影响肺部健康，引发肺部疾病，因此研究细颗粒物在人体肺腺泡区的沉积规律对肺部疾病致病机制的探究以及肺部疾病的预防和治疗都有重要意义。
人体肺腺泡区结构复杂、尺寸微小，尚未获得真实的计算机模型。现行数值模拟研究多采用简化的计算机模型，主要分为单肺泡模型和多肺泡模型。单肺泡研究仅对单个肺泡建模，计算量较小，着重关注单个肺泡的内部特性。Sznitman 等[1]对处于不同级数的单肺泡进行了数值模拟研究，发现肺泡内流场形态随着级数增加由循环状向辐射状变化。单肺泡研究只能关注单个肺泡内流场形态和颗粒物沉积状况，无法探究肺腺泡区的整体特性，因此国内外学者纷纷建立多肺泡模型进行进一步研究。目前，多肺泡模型主要分为四类[2]：“柱状气管+球状肺泡”的多肺泡模型、基于十四面体的多肺泡模型、泰森多边形法构建的多肺泡模型、计算机成像和重建技术构建的多肺泡模型。Ma 等[3]建立了“柱状气管+球状肺泡”的 5 级多肺泡模型，发现呼吸周期、潮气量和重力等对颗粒物沉积有较大影响。Longest 等[4]建立了 3 个基于十四面体的多肺泡模型，发现屏气能大幅度增加颗粒物的沉积率。Hofemeier 等[5]采用泰森多边形法构建了人体 1/8（第 3～8 级）肺腺泡模型，研究了 0.005～5.0 μm 粒径范围内不同大小颗粒物的沉积状况，发现颗粒物沉积率随粒径增加呈先下降后上升的趋势。Sera 等[6]采用微计算机断层扫描（micro focal computed tomography system，micro-CT）技术构建了大鼠肺腺泡模型，发现 0.8～1.0 μm 粒径颗粒物多沉积于肺腺泡入口处，而 0.1～0.2 μm 粒径颗粒物则主要沉积于模型中部和末端。
上述数值模拟研究大多关注颗粒物在人体肺腺泡区局部或整体的最终沉积率，而对其沉积动态过程鲜有涉及，无法进一步了解人体肺腺泡区颗粒物的运动规律。本文建立了五个“柱状气管+球状肺泡”多肺泡模型，通过数值模拟研究了模型级数和结构、颗粒物粒径、呼吸模式等因素对 0.1～5 μm 粒径范围内颗粒物沉积规律的影响，总结了颗粒物沉积率和沉积速度的变化规律，归纳了沉积曲线特征，以期为进一步研究肺部疾病致病机制和相关疾病的预防与治疗提供理论依据和数据支持。
1 模型与方法
1.1 人体肺腺泡模型
在经典的全肺模型（Weibel A 模型[7]）中，肺部分为 0～23 级。其中，17～23 级部分含有肺泡，是呼吸作用的主要区域，为肺腺泡区。人体肺腺泡区尺寸数据如表 1 所示。
表1
人体肺腺泡区尺寸数据
Table1.
Geometric data of human pulmonary acinus
由表 1 数据可知，肺腺泡区尺寸数据随级数变化显著。从 17 级至 23 级，肺腺管长度、管径数据逐级递减，最大缩减率分别为 64.54%、24.07%；肺泡直径和单节肺腺管肺泡数量逐级递增，最大增长率分别为 79.25%、300%。Sznitman 等[1]研究发现模型尺寸是影响数值模拟结果的重要因素。十四面体多肺泡模型中肺泡形状、大小一致[8]，仅适用于肺腺泡区局部特性探究或尺寸影响可忽略不计的情况，因而本文采用“柱状气管+球状肺泡”架构，将肺腺管简化为圆柱状，肺泡简化为附着在肺腺管表面的球，建立了如图 1 所示的多肺泡模型。
图1
“柱状气管+球状肺泡”架构多肺泡模型
Figure1.
Multi-alveolar model of “columnar trachea + globular alveolar” structure
然后，在该模型的基础上依据表 1 尺寸数据建立了图 2 所示的 G3（17～19 级）、G4（17～20 级）、G5（17～21 级）、G6（17～22 级）、G7（17～23 级）共 5 个多肺泡模型。各模型所含级数、肺泡数、长度、体积等参数如表 2 所示。
图2
多肺泡模型
Figure2.
Multi-alveolar models
表2
多肺泡模型参数
Table2.
Multi-alveolar models parameters
1.2 数值模拟方法
1.2.1 数值模拟软件
采用 Inventor 2015 进行多肺泡模型建模，采用 ICEM CFD 15.0 进行模型网格划分，采用 FLUENT 15.0 进行数值模拟计算。
1.2.2 网格划分结果
采用非结构网格对 G3～G7 模型进行网格划分，各模型网格总数及最小网格质量如表 3 所示，最小网格质量均大于 0.3，可以用于数值模拟计算。
表3
G3～G7 模型网格划分结果
Table3.
Meshing results of G3―G7 models
1.2.3 控制方程
① 流体控制方程：人体肺腺泡区气流运动为等温、不可压缩定常流动。满足质量、动量守恒方程。
质量守恒方程：
|
式中，表示流体密度， = 1.185 kg/m3，为流体速度，t为时间。故上式可表示为：
|
动量守恒方程：
|
式中，P为流体压力， = 1.785e-5 表示运动粘度，为流体体积力矢量。
② 颗粒物控制方程：本研究中，忽略 Magnus 力、Saffman 力、Basset 力等对颗粒物运动影响较小的力和颗粒物之间的相互影响，仅考虑重力、布朗力和流体曳力的作用效果，因此单个颗粒的运动控制方程为：
|
|
|
|
其中，其中 和 分别为颗粒物的密度和直径， 为随机数， 为玻尔兹曼常数， 为空气热力学温度， 为颗粒半径， 为计算时间间隔。
1.2.4 边界条件设置
如图 3 所示，以 G3 模型区域划分为例，本研究中将模型分为两部分：顶端平面（黑色阴影部分）和其他部分。其中，顶端平面为气流和颗粒物入口，采用入口速度函数实现气流和颗粒物的进入和排出模型，以此模拟呼吸作用下人体肺腺泡区的气流和颗粒物的运动。其他部分为无滑移壁面，颗粒物在运动过程中接触壁面会被捕捉，即为沉积。
图3
G3 模型区域划分
Figure3.
Parts of G3 model
1.2.5 入口速度函数
Oakes 等[9]测量了 21 岁健康男性休息状态下的呼吸流量曲线，近似为正弦曲线，因此本研究中入口速度采用随时间变化的正弦函数，速度正方向为 − Z，入口速度 可表示为:
|
其中， 为入口气流最大速度， 为呼吸周期。
则数值模拟一周期模型内气体最大进入量为：
|
积分可得:
|
其中，为入口管径，即 17 级肺腺管管径。
以上述建立的多肺泡模型代表人体肺腺泡进行数值模拟研究，模型体积对应功能残气量（functional residual capacity，FRC），最大进气量对应潮气量，则有：
|
根据 Haefeli-Bleuer 等[10]的测量结果，取人体正常呼吸状态下功能残气量为 3 L，潮气量为 600 mL，呼吸周期为 4 s。依据式（1）（2）（3）计算可得，正常呼吸状态下（ = 600 mL，T = 4 s），G3～G7 模型入口气流最大速度和入口速度如表 4 所示。
表4
正常呼吸状态下 G3～G7 模型入口气流最大速度和入口 速度
Table4.
The maximum velocity and inlet velocity of inlet airflow in G3−G7 models under normal respiratory mode
1.2.6 呼吸模式
当人体处于不同呼吸状态时，潮气量和呼吸周期存在明显差异。由式（1）（2）（3）可知，潮气量、呼吸周期能直接决定模型入口速度函数，进而影响颗粒物在人体肺腺泡区的沉积规律。从肺腺泡形态结构来说，G7 模型包含级数更多，肺泡数量也更多，更加符合人体实际，因此本研究基于 G7 模型进一步研究呼吸模式对人体肺腺泡区颗粒物运动规律的影响。取呼吸周期为 4 s，潮气量分别为 300、600、900 mL，代表人体处于浅呼吸、正常呼吸和深呼吸状态。取潮气量为 600 mL，呼吸周期分别为 2、4、6 s，代表人体处于急促呼吸、正常呼吸和缓慢呼吸状态，由式（1）（2）（3）计算可得，不同呼吸模式下 G7 模型入口速度如表 5 所示。
表5
不同呼吸模式下 G7 模型入口速度
Table5.
The velocity of inlet airflow in G7 model under diffe rent respiratory mode
1.2.7 颗粒物设置
数值模拟计算中，将颗粒物简化为球状，密度为 1 kg/m3，选取直径为 0.1、0.3、0.5、0.8、1、3、5 μm 等能够进入人体肺腺泡区的颗粒物[11]作为研究对象。采用 FLUENT DPM（Discrete Phase Model）模型，激活重力、布朗力、曳力选项。数值模拟开始时，在模型入口平面基于网格释放颗粒物，使颗粒物均匀地分布在气流中并随气流进入模型。利用 DPM 模型中颗粒物沉积数量统计功能，每 时间间隔记录一次颗粒物沉积数量，计算总时长为一个周期。定义 为t时刻颗粒物沉积率， 为 t 时刻颗粒物沉积数量， 为释放的颗粒物总量，则：
|
1.2.8 颗粒物沉积速度评价参数
通常而言，颗粒物沉积速度评价参数是首次沉积发生的时间，即颗粒物首次沉积发生的时间越短，沉积速度越快。但实际上，研究中多采用间隔采样法对颗粒物沉积率进行记录，无法得到颗粒物首次沉积发生的具体时间。而且，颗粒物沉积为非线性，颗粒物首次沉积时间并不能真实反映该粒径颗粒物沉积的快慢。为了更好地评价颗粒物在人体肺腺泡区沉积速度差异，引入颗粒物沉积速度评价参数：1/4 沉积时间，即某粒径颗粒物沉积率达到 1/4 最终沉积率的时间，记为 。本研究中数值模拟时间为一周期，T 时刻颗粒物沉积率即为最终沉积率 ，则 1/4 最终沉积率为 。假定 、（）为两相邻记录时刻，即 ，且这两时刻沉积率满足：
|
那么，1/4 沉积时间位于 ～ 之间，假定 ～ 内颗粒物沉积速度恒定，为线性沉积，则可利用插值计算估算得：
|
化简得：
|
1/4 沉积时间以颗粒物达到 1/4 最终沉积率的时间作为沉积快慢评价参数，能对颗粒物沉积快慢进行较好地评价。
2 结果和讨论
2.1 模型对人体肺腺泡区吸入颗粒物沉积特性的影响
2.1.1 模型级数影响
模型的建立是数值模拟研究的前提和基础，也是影响数值模拟结果的重要因素。由表 2 可知，G3～G7 模型包含肺腺泡级数不同，模型总长度和体积也存在差异，因此研究 G3～G7 模型级数差异对颗粒物沉积特性的影响十分必要。正常呼吸模式下，G3～G7 模型中各粒径颗粒物最终沉积率如图 4 所示。由图 4 可知，由 G3 模型至 G7 模型，随着模型级数的增加，0.1 μm 粒径颗粒物最终沉积率呈微弱上升趋势，整体差异不大；0.3、0.5、0.8、1 μm 粒径颗粒物最终沉积率呈明显上升趋势；3、5 μm 粒径颗粒物最终沉积率则呈现出下降趋势。
图4
G3～G7 模型中各粒径颗粒物最终沉积率
Figure4.
Final deposition fraction of particles in G3―G7 models
研究表明，颗粒物在人体肺腺泡区沉积机制主要为布朗扩散、重力沉降和惯性冲击，而惯性冲击是影响 0.3～1 μm 粒径颗粒物沉积的重要因素[12]。由表 4 可知，正常呼吸状态下，模型所含级数越多，体积越大，则模型入口速度越大，惯性冲击对颗粒物沉积的影响也相应越大。因此随着模型级数的增多，0.3～1 μm 粒径颗粒物最终沉积率呈现出明显上升趋势。0.1 μm 粒径颗粒物沉积由布朗扩散主导，惯性冲击影响不大，因此沉积率仅微弱上升。由表 2 可知，模型级数的增加不仅体现为体积的增加，还体现在模型总长度的增长。3、5 μm 粒径颗粒物沉积由重力沉降主导，布朗扩散和重力沉积影响较小，随着模型总长度的增加，颗粒物沿重力方向沉积所需运动路程相应增加，沉积所需运动时间更长，最终沉积率也相应更低，因而最终沉积率随模型级数增加呈现出明显的下降趋势。Katan 等[13]在研究年龄对肺腺泡区颗粒物沉积状态的影响时发现，成人肺腺泡模型尺寸较大，颗粒物沉积路径增加，不利于颗粒物的沉积，因而沉积率较低，与本文研究结果一致。
2.1.2 模型结构影响
由图 4 可知，上一小节关于模型级数对颗粒物沉积影响规律的描述，只是整体趋势的归纳和总结。G6 模型中，0.1、0.3、0.5 μm 粒径颗粒物沉积率较 G5 模型中沉积率稍低，分别相差 2.23%、2.13%、2.90%，鉴于沉积率差值较小，可以认为基本不影响上一小节趋势的归纳。而 G3 模型中，3 μm 粒径颗粒物的最终沉积率则与理论值相差甚远，仅为 10.53%（理论值应大于 G4 模型中 3 μm 粒径颗粒物最终沉积率，即 > 46.49%），这是模型几何结构对颗粒物沉积影响的典型体现。如图 5 所示为 = 2 s 时刻 G3、G4 模型中 3 μm 粒径颗粒物的分布情况。正常呼吸模式下，2 s 时刻颗粒物运动路径最长，最容易接触壁面并沉积。而 G3 模型中，2 s 时刻颗粒物距接触壁面仍有一定距离，因此沉积率较小。G4 模型中，入口速度增大，运动路径增加，能够接触壁面并沉积，因而沉积率较高。同理，缩减模型中 17 级肺腺管长度（）或改变重力方向[14]也能大幅度提升 G3 模型中 3 μm 颗粒物的沉积率。
图5
t = 2 s 时刻 G3、G4 模型中 3 μm 粒径颗粒物分布情况
Figure5.
Distribution of 3 μm particles in G3 and G4 model at t = 2 s
2.2 颗粒物粒径对人体肺腺泡区吸入颗粒物沉积特性的影响
正常呼吸模式下，G3～G7 模型中颗粒物最终沉积率随粒径变化趋势如图 6 所示，不同模型间各粒径颗粒物沉积曲线略有差异，但基本形状均为“U”型，即在颗粒物粒径为 0.1 μm 时，颗粒物沉积率较高，随着颗粒物粒径的增加，沉积率先下降、后上升，并在颗粒物粒径为 1 μm 左右时沉积率最低，呈现出两端高、中间低的趋势，这一趋势与文献数据[5, 15-16]保持一致。这是因为，布朗力和重力是影响颗粒物沉积的主要因素。布朗力对微小颗粒物（< 0.1 μm）沉积起主导作用，重力对大颗粒物（> 2 μm）沉积起主导作用。因此，当颗粒物粒径较小时，沉积受布朗力主导，沉积率较高。颗粒物粒径较大时，沉积受重力主导，沉积率较高。随着颗粒物粒径增加，布朗力作用减弱、重力作用加强，而 1 μm 左右粒径颗粒物受两者影响均较小，因而沉积率最低，故沉积曲线呈现为“U”型。
图6
G3～G7 模型中颗粒物最终沉积率随粒径变化趋势
Figure6.
Final deposition fraction of particles as a function of diameter in G3—G7 models
2.3 人体肺腺泡区吸入颗粒物动态沉积特性
2.1 小节探讨了模型级数对人体肺腺泡区颗粒物沉积特性的影响，从肺腺泡形态结构来说，G7 模型包含级数更多，肺泡数量也更多，更加符合人体实际，因此本研究基于 G7 模型进一步研究人体肺腺泡区吸入颗粒物的动态沉积特性和呼吸模式对颗粒物运动规律的影响。正常呼吸模式下（ = 600 mL，T = 4 s），G7 模型中各粒径颗粒物沉积率随时间变化曲线如图 7 所示，G7 模型中各粒径颗粒物沉积率均随时间增加，至 4.0 s 时刻达到最大值，即最终沉积率。其中，0.1 μm 颗粒物前期沉积速度快，中后期沉积缓慢并趋于稳定，沉积率变化曲线呈对数型分布。5 μm 颗粒物前期沉积速度慢，中期沉积快速并趋于稳定，后期沉积缓慢，沉积曲线呈“S”型。0.3 μm 颗粒物沉积趋势与 0.1 μm 基本相似，沉积率相对较低，对数型沉积曲线相对较矮。3 μm 颗粒物沉积曲线与 5 μm 基本相似，沉积率较低，呈较矮的“S”型曲线。其他粒径颗粒物前期沉积缓慢，中后期沉积率缓慢增长，沉积曲线基本呈折线型。
图7
G7 模型各粒径颗粒物沉积率变化曲线
Figure7.
Variation curves of deposition fraction of particles in G7 model
由图 7 沉积曲线可以看出，0.0～1.0 s 间除 0.1 μm 颗粒物沉积率增长较快外，其他粒径颗粒物沉积率均较小，很难判断各粒径颗粒物沉积速度差异，可见引入合理参数对沉积速度进行评价的必要性。G7 模型中各粒径颗粒物 1/4 沉积时间如图 8 所示，从 1/4 沉积时间数值来看，0.1、0.3、5 μm 粒径颗粒物 1/4 沉积时间相对较短，分别为 0.60、1.03、1.33 s。由此可见，微小颗粒物沉积速度最快，大颗粒物沉积速度次之，而 0.5～3 μm 颗粒物沉积速度则相对较慢。这是因为，数值模拟初期，气流速度较小，颗粒物速度较小，惯性冲击影响较小，重力沉降则需要颗粒物沿肺腺管运动一定时间至模型分叉处才能产生沉积。因此，数值模拟初期颗粒物沉积受布朗扩散主导，微小颗粒物沉积速度最快；大颗粒物沉积受重力主导，运动一定时间接触壁面后开始沉积，沉积速度也相应较快；而其他粒径颗粒物因数值模拟初期气流速度小，惯性冲击影响作用小，沉积较慢。
图8
G7 模型中各粒径颗粒物 1/4 沉积时间
Figure8.
1/4 deposition time of particles in G7 model
2.4 呼吸模式对人体肺腺泡区吸入颗粒物沉积特性的影响
2.4.1 潮气量影响
潮气量分别为 300、600、900 mL 时，G7 模型中各粒径颗粒物最终沉积率和 1/4 沉积时间如图 9 所示。
图9
不同潮气量下颗粒物的最终沉积率和 1/4 沉积时间
Figure9.
Final deposition fraction and 1/4 deposition time of particles under different tidal volume
由图 9 可知，三种潮气量下，颗粒物最终沉积率随粒径变化趋势均呈“U”型，1/4 沉积时间的变化趋势也基本保持一致，说明在人体肺腺泡区潮气量的变化不会改变流场状态[17]，仍为充分发展的层流。同一粒径颗粒物在潮气量为 300 mL 时最终沉积率最低，1/4 沉积时间最长，潮气量为 600 mL 时最终沉积率和 1/4 沉积时间次之，潮气量为 900 mL 时最终沉积率最高，1/4 沉积时间最短。潮气量对颗粒物沉积规律的影响规律可概括为：潮气量越大，颗粒物沉积率越高，沉积速度越快。Hofemeier 等[5]分别研究了潮气量为 500 mL 和潮气量为 2 500 mL 下的颗粒物沉积状况，发现潮气量为 2 500 mL 时颗粒物在肺腺泡内总体沉积率更高，与本文研究结果保持一致。如表 5 所示，潮气量的增加会导致模型入口速度的增加，受流场曳力影响，颗粒物运动速度随之增加，惯性冲击影响作用加大，因而颗粒物沉积数增多，沉积速度加快。
值得注意的是，不同粒径颗粒物沉积规律率受潮气量影响程度不同，0.1 和 5 μm 粒径颗粒物受潮气量影响较小，沉积率极差量分别为 5.47% 和 6.28%，1/4 沉积时间极差量分别为 0.26 s 和 0.42 s。0.3～1 μm 粒径颗粒物则受潮气量影响较大，沉积率极差量分别为 12.12%、18.17%、23.25%、20.85%，1/4 沉积时间极差量分别为 0.68、0.80、0.83、0.98 s。这一结果说明惯性冲击是影响 0.3～1 μm 颗粒物沉积的关键因素，对微小颗粒物和大颗粒物的沉积则影响较小。
2.4.2 呼吸周期影响
呼吸周期分别为 2、4、6 s 时，G7 模型中各粒径颗粒物的最终沉积率和 1/4 沉积时间如图 10 所示，呼吸周期对颗粒物沉积规律的影响与潮气量影响较为一致，其影响规律可概括为：呼吸周期越长，颗粒物最终沉积率越大，沉积速度越慢，这一结果与文献[18-19]保持一致。这是因为随着呼吸周期的增加，数值模拟时间加长，颗粒物沉积更加充分，沉积率相应更高。由表 5 可知，呼吸周期的增加还会导致模型入口气流速度的减小，惯性冲击作用减小，因而颗粒物沉积速度减缓，1/4 沉积时间更长。
图10
不同呼吸周期下颗粒物的最终沉积率和 1/4 沉积时间
Figure10.
Final deposition fraction and 1/4 deposition time of particles under different respiratory cycle
3 结语
本文一共建立了 G3～G7 五个多肺泡模型，采用数值模拟方法研究了模型级数和结构、颗粒物粒径、呼吸模式等因素对 0.1～5 μm 粒径范围内 7 种颗粒物沉积率和沉积速度的影响。结果表明：模型级数和几何结构等是造成模型间最终沉积率差异的重要因素。颗粒物沉积率随粒径变化曲线呈“U”型，其中，细微颗粒物（0.1 μm）沉积由布朗扩散主导，沉积率高，沉积速度快，沉积曲线呈对数型分布；大颗粒物（5 μm）沉积由重力沉降主导，沉积率高，沉积速度较快，沉积曲线呈“S”型分布；呼吸模式对颗粒物沉积规律的影响体现为：潮气量越大，颗粒物沉积率越高，沉积越快；呼吸周期越长，颗粒物沉积率越高，沉积越慢，这一结果与文献数据保持一致。本文建立的“柱状气管+球状肺泡”多肺泡模型和研究方法能为人体肺腺泡区吸入颗粒物运动规律的进一步探究提供参考和借鉴，研究结果和数据能为肺部疾病致病机制的探究、肺部疾病的预防和治疗提供理论依据和数据支持。
本文数值模拟过程中，模型壁面刚性固定，以模型入口气流速度的正弦变化模拟呼吸过程，忽略了真实呼吸过程中肺腺泡模型边界运动对结果的影响，后期还需进一步研究模型边界运动的影响。此外，颗粒物沉积的影响因素选取还需进一步细究，考虑更多其他因素，如重力[20]、摒气时间[21]、肺部病变[22]等的影响。
利益冲突声明：本文全体作者均声明不存在利益冲突。
