1. 引言
肺癌的早期诊断和治疗仍然存在重大挑战。随着新兴产业的出现，伴随着的是大范围的雾霾天气以及呈现不同程度污染的空气质量，持续不断的雾霾天气导致呼吸道疾病患者爆发性增长 [1] 。从毒理学的角度来讲，所有直径小于10微米的颗粒都具有生物活性，可能导致敏感个体的过敏反应甚至癌症 [2] 。研究表明，环境条件的急剧变化、不健康的生活方式以及暴露于污染物和辐射 [3] 等因素导致预计未来20年肺癌病例数可能增加近50% [4] 。也就是说，吸入有毒和无毒的气溶胶粒子对人肺部的作用效果及其生理机制的影响存在显著差异。雾化治疗是最广泛使用的治疗和预防肺部疾病的方法。近年来，由于其优点，雾化治疗已越来越被认为是一种合适的用于治疗肺部和非肺部疾病的方法。医学对于治疗性吸入药物的传递和沉积有靶向区域的要求，因此了解和掌握粒子的在肺部的沉积情况有利于肺部疾病的有效治疗 [5] [6] [7] 。而且，为了成功地进行肺部给药，最小化药物剂量和低副作用，转移过程必须有针对性 [8] 。
当人处于不同的运动状态时，对应的呼吸频率和呼吸量是不同的，吸入的粒子浓度也是不同的。因此研究定量粒子在人体肺部的沉积效果，对正确认识到粒子沉积对人体肺部的影响具有重要意义 [9] ，在量化吸入颗粒的健康风险 [1] [10] 和药物治疗上有一定的指导意义。
国内外对粒子在人肺部的粒子沉积有多方面研究。Cheng等人通过尸体解剖制作的呼吸系统模具研究了不同呼吸模式下的粒子沉积 [11] ，Rahman等人研究了现实人体肺气道中的汽车尾气、烟雾和灰尘颗粒的输运情况 [2] ，另外，Dogan等人研究了人体肺腺泡内三种不同场景下的气溶胶行为 [12] 。上述研究均是利用单一的湍流模型在固定的人体肺气道中研究不同运动状态下不同粒子的输运和沉积情况，缺少对不同湍流模型下肺气道部分粒子沉积的比较研究，无法确定最合适实际情况的湍流模型。
本文以一个人体“左二右三”结构的肺部支气管模型为研究对象，考虑典型的湍流模型(k-ω、k-ε)和层流模型下对15 L/min、30 L/min和60 L/min的恒定吸气流速下颗粒在支气管模型中的沉积情况。本文研究所用的模型是半理想化口腔通道和非对称3级肺气管构成的更真实的人体上呼吸道模型 [13] ，该模型能够很好地将多相流理论应用于粒子在人体呼吸系统中的动力学行为研究。结合肺气道的对比模拟结果和本文研究的仿真情况，说明1) k-ω模型比其他模型的仿真情况更接近真实的肺部支气管粒子沉积情况。2) 粒子的沉积与斯托克斯数的大小有一定关系，“左二右三”支气管各部分的粒子沉积各不相同，随着吸入粒子速度的增大，粒子沉积在各隆凸脊的位置逐渐变为左肺叶的靠下位置。3) 粒子的沉积位置和粒子的运动状态有关。
2. 模型理论
人体的呼吸道主要有三个部分：胸外部分、气管和支气管区域和腺泡区，气管支气管区也称为传导区(第0~16代)，主要将吸入的空气传导至深肺位置。此外，称为过渡区的腺泡区域 [14] (第17~23代)和称为呼吸区的第20~23代 [15] ，呼吸区主要是发生气体交换。气管支气管(TB)和肺泡(AI)区域通常称为“胸内(IT)区域” [16] 。本文从1956年Weibel建立的一个完全对称的理想化肺部分叉树状结构出发，考虑人体所处环境的差异性，并且具体到肺部的每级的不对称支气管结构 [17] ，因此，选择恰当合适的大小、长度的支气管模型对实验的研究显得尤为重要。本文基于计算机流体力学和两相流理论的基础 [18] ，利用ANSYS CFX研究气溶胶粒子在人体内部运动动力学行为。
本文从人体实际情况出发 [19] ，利用CFX软件建立人体“左二右三”肺部支气管模型(肺部不对称支气管模型)，详见图1，该模型已在文献 [13] 中已经证明更接近真实的人体肺部模型，有关模型的具体参数详见文献 [13] 。本文的肺部结构忽略了与重力的夹角，即各级肺叶支气管的轴线位于同一平面内。
3. 数值方法
本文基于计算流体力学的数值研究 [21] ，利用商业软件ANSYS Workbench环境下的CFX 18.0计算程序进行连续性方程、动量方程以及k-ω和k-ε方程运算。假设稳定的吸入气流状态，入口处均匀入射等量颗粒，颗粒与流体的初始速度相等，气道壁上假设无滑动，空气流场和颗粒流场之间假设为单向耦合，数值模拟采用有限体积法，全隐式耦合算法求解湍流方程。计算在DELL-R740专用服务器上进行，该服务器规格为CPU4116/64G/6T。
当残差收敛曲线的所有残差值达到10e−4时，假设流场的稳态解收敛。使用收敛的流场解决方案，每个斯托克斯数对应的颗粒轨迹模拟运行时间为约1小时。运行典型湍流模型(k-ω模型和k-ε模型)的运行时间相差1~2小时，具体取决于入口雷诺数。
4. 模拟结果与讨论
本文通过改变粒子直径和改变粒子速度值的方式来改变斯托克斯数的大小，主要研究三种运动状态，分别是休息状态吸气流速(15 L/min)、轻微呼吸状态吸气流速(30 L/min)和运动状态吸气流速(60 L/min)的粒子沉积情况，此外，还分析了层流和湍流模型下支气管部分的流场和粒子沉积情况。
在1999 Cheng实验 [11] 中给出了上气道部分粒子沉积率随斯托克斯数变化而变化的曲线关系图。我们对上气道部分的粒子沉淀情况进行数值模拟，得到与实验较为吻合的结果，验证了模型的有效性。在对肺部支气管部分的模拟中，我们定义支气管部分颗粒沉积率为特定支气管部分中沉积颗粒的数量与进入该部分的颗粒的数量之比。
Stokes Number (斯托克斯数)的定义式 [11] ：
(1)
其中ρ为粒子密度，d为粒子的直径，U为口腔入口的平均速度，具体表达式为U = Qin/Smean，Smean为口咽通道的平均横截面积，μ为空气粘度系数，Dh为最小水力直径，Dh = 0.00817 m [13] 。
本文主要分析层流、k-ω和k-ε三种代表性模型对“左二右三”结构的肺部支气管部分的粒子沉积情况。
4.1. 流场特性
首先，我们对整体模型进行纵向切面，分析了“左二右三”结构的支气管肺叶部分在此切面上的整体流场情况。图2显示了在层流、k-ω和k-ε模型下入口流速为15 L/min、30 L/min和60 L/min的速度矢量图，图3为对应的流场线图。
Figure 2. Velocity vector diagrams of the bronchial section under laminar, k-ω, and k-ε models at inlet flow velocities of 0.8 m/s, 1.6 m/s, and 3.2 m/s
图2. 层流、k-ω和k-ε模型下入口流速为0.8 m/s、1.6 m/s、3.2 m/s时支气管部分速度矢量图
从休息状态条件(15 L/min对应0.8 m/s)下人体支气管部分中的速度矢量方向分布图来看，支气管入口处流量大，分流后各分管中流量各不相同。在层流模式速度矢量未到达壁面位置；从轻微呼吸状态吸气流速状态条件(30 L/min对应1.6 m/s)下的人体支气管部分中的速度矢量方向分布图来看，人体支气管部分在k-ε模型模拟下速度矢量最密集，层流和k-ω模型模拟速度矢量图相似；从运动状态吸气流速(60 L/min对应3.2 m/s)时的人体支气管部分中的速度矢量方向分布图来看，随着速度的增大，支气管入口处速度矢量密度也逐渐增加。
从速度矢量图的流线图更好的能看出粒子的速度运动情况，在支气管管壁位置k-ω模型模拟有明显的沿壁面粒子运动的情况；层流模型模拟时在右侧支气管上部位置存在流线图交叉的情况；k-ε型模拟结果上，各流线较为规整有序，随速度的增大，流场线的密度越大。
Figure 3. Streamline diagrams for the bronchial section under laminar, k-ω, and k-ε models at inlet velocities of 0.8 m/s, 1.6 m/s, and 3.2 m/s
图3. 层流、k-ω和k-ε模型下入口流速为0.8 m/s、1.6 m/s、3.2 m/s时支气管部分速度流线图
4.2. 分支肺部支气管粒子沉积
为描述不同斯托克斯数下的三种模型模拟在支气管部分的颗粒沉积情况，将支气管的分支气管粒子沉积进行对比(图4)。可见，在入口速度为0.8 m/s的情况下，三种模型模拟时支气管右下(L42)位置部分粒子沉积最大，左下(L33)部分k-ε和k-ω模型模拟粒子沉积率基本相同；在入口速度为1.6 m/s的情况下，层流和k-ω模型模拟在左下(L33)部分粒子沉积率最大，k-ε模型模拟在右下(L42)位置部分粒子沉积最大，左上(L34)部分三种模型模拟的粒子沉积率相同；在入口速度为3.2 m/s的情况下，层流模型模拟在右下(L42)部分粒子沉积率最大，其次为左下(L33)，k-ε和k-ω模型模拟在左下(L33)位置部分粒子沉积最大，右下(L42)和左上(L34)部分模拟的粒子沉积率相同。
图5对支气管层流和湍流模型下同一斯托克斯数(2.87e−1)颗粒在各气道模型中的沉积情况进行柱状对比，层流模型模拟下3.2 m/s时右中(L41)部分右粒子沉积率最大；k-ε模型对应的1.6 m/s时右上(L31)粒子沉积率最大，而在右下(L42)位置粒子沉积率最小，3.2 m/s时右下(L42)粒子沉积率比右中(L41)粒子沉积率大，但相差不大；k-ω模型模拟下1.6 m/s时右上(L31)部分有粒子沉积，右下(L42)和左下(L33)部分粒子在三种情况对应的流速下都有粒子沉积，且左下(L33)部分粒子沉积率最大。
Figure 4. Comparison of particle deposition rates at specific locations within the bronchial lobes under laminar, k-ω, and k-ε models at an inlet density of 1000 kg/m3 and velocities of 0.8 m/s, 1.6 m/s, and 3.2 m/s, corresponding to particle diameters of 1.0e−6 m and Stokes Numbers of 5.93e−4, 1.19e−3, and 2.37e−3, respectively
图4. 层流、k-ω和k-ε模型下入口处吸入密度为1000 kg/m3，速度分别为0.8 m/s、1.6 m/s、3.2 m/s时的粒子沉积在支气管肺叶的具体位置粒子沉积率对比(对应粒子直径为1.0e−6 m，对应斯托克斯数为5.93e−4、1.19e−3、2.37e−3)
Figure 5. Comparison of particle deposition rates at specific locations within the bronchial lobes for laminar, k-ω, and k-ε models at inlet inhalation velocities of 0.8 m/s, 1.6 m/s, and 3.2 m/s, corresponding to a Stokes Number of 2.87e−1
图5. 层流、k-ω和k-ε模型下入口吸入速度分别为0.8 m/s、1.6 m/s、3.2 m/s时的粒子沉积在支气管肺叶部分具体位置的粒子沉积率对比(对应斯托克斯数为2.87e−1)
通过这些模拟结果，我们观察到随着斯托克斯数的增大，粒子沉积率最高的区域从右下(L42)位置逐渐转移到左下(L33)位置。这一发现表明，不同的模拟模型对于预测支气管部分的粒子沉积效率存在显著差异。
4.3. 粒子沉积图
人体肺部支气管粒子的沉积情况受到多种因素的影响，包括流体动力学特性、气道的几何结构、颗粒的物理特性(如大小、形状和密度)以及呼吸模式等。为进一步研究在不同呼吸状态下人体肺部支气管粒子的明显沉积情况，我们利用CFX模拟时吸入10000个粒子，沉积结果见图6。
研究表明，在层流条件下，较低的流速(0.8 m/s)导致粒子在支气管上部直管区域沉积较多，这可能是因为在低流速下，粒子的惯性效应较小，更容易受到重力作用而在管道上部沉积。随着流速的增加(1.6m/s和3.2 m/s)，粒子在隆凸脊周围的沉积率增加，这可能是由于流速增加时，湍流效应增强，使得粒子更容易被输送到管道壁面并沉积。在k-ε模型下，粒子沉积随速度的增大而增大，但这种增长并不明显。这意味着在该湍流模型中，虽然流速增加有助于粒子的混合和沉积，但湍流的复杂性可能也导致了粒子再悬浮，从而抵消了部分沉积效果。在k-ω模型下，随着速度的增大，粒子沉积模式发生了变化，尤其是在垂直于重力方向的管壁上沉积的粒子明显增加。这表明在较高流速下，湍流和粒子的惯性效应共同作用，使得粒子更容易在管道的下半部分沉积，这是由于重力和湍流的联合作用。
从对比图来看，上游沉积和局部气流条件对粒子在支气管部分的沉积分布有显著影响。在支气管分叉处，由于流体动力学的复杂性，粒子会在这些区域积累，导致沉积率较高。数据表明，沉积率在支气管的不同区域存在局部变化，这与局部的流场特性、管道几何形状以及粒子的物理特性有关。
Figure 6. Particle deposition under laminar, k-ε, and k-ω models at inlet inhalation velocities of 0.8 m/s, 1.6 m/s, and 3.2 m/s, corresponding to a Stokes Number of 2.87e−1
图6. 层流、k-ε和k-ω模型下入口吸入速度分别为0.8 m/s、1.6 m/s、3.2 m/s时的粒子沉积情况(对应斯托克斯数为2.87e−1)
5. 结果
本文研究了稳态条件下人的休息、轻微运动和剧烈三种运动状态下人体“左二右三”结构的肺部支气管部分粒子的运动和沉积情况，主要利用CFX对层流和典型的湍流模型下不同的仿真结果进行对比分析，基于研究情况，可得出以下结论。
1) 从不同速度矢量图分析粒子运动情况，k-ω模型模拟时存在沿壁面粒子运动的情况；k-ε模型模拟显示，随着速度的增加，各运动轨迹的流线密度增大；对比实验研究结果，k-ω模型比其他模型的仿真情况更接近真实的肺部支气管粒子沉积情况。
2) 粒子的沉积与斯托克斯数的大小有一定关系，“左二右三”支气管各部分的粒子沉积各不相同，考虑到肺部支气管各分支的尺寸大小以及支气管部分独特的结构对粒子沉积有一定影响，各隆凸脊的周围成为支气管部分粒子沉积率较高部位。随着吸入粒子速度的增大，粒子沉积在各隆凸脊的位置逐渐变为左肺叶的靠下位置。
3) 粒子的沉积位置和粒子的运动状态有关，对于分析研究医学上特定的颗粒药物在肺部支气管的输送治疗效果能够提供一定的帮助，也可能对肺部支气管毒理学影响提供一定信息。