CN104050321A - 肺泡内颗粒运动轨迹的检测方法 - Google Patents

Abstract

一种肺泡内颗粒运动轨迹的检测方法，包括：根据肺泡的几何数据对肺泡进行网格建模，建立肺泡动网格模型；对所述肺泡动网格模型中进行空气流场模拟，获得空气流场的速度；根据肺泡的空气流场的速度对肺泡的颗粒进行模拟，获得颗粒数值模拟模型；根据颗粒数值模拟模型确定稳态和非稳态气流下不同受力条件下不同颗粒在肺泡内的运动轨迹。本发明方案提高了获取颗粒在肺泡内的运动轨迹的准确率，进而可以根据不同受力条件不同颗粒大小下颗粒的运动轨迹确定颗粒沉积规律，从而可以提高气溶胶的药性。

Description

肺泡内颗粒运动轨迹的检测方法

技术领域

本发明涉及检测技术领域，特别是涉及一种肺泡内颗粒运动轨迹的检测方法。

背景技术

人体呼吸的主要功能是为身体的各个组织提供氧气和排除二氧化碳废气，人体的呼吸过程可以分为两个阶段：由外界环境向血液输送气体；气体经由血液进入到各个组织。随着社会、经济的不断进步，人类对生活环境的质量要求不断提高，对生存环境的保护意识逐渐增强。工业生产和环境恶化所带来的颗粒物污染已经成为评价生活质量和大气质量的一个重要指标之一。大气气溶胶颗粒物污染中，其中的部分微小的气溶胶颗粒，尤其是可吸入颗粒对人类的健康的影响更是深远，在进入人体呼吸道后，没有沉积在呼吸道的传导气管上，而是深入到人体呼吸道终末处的气体交换区沉积，很多研究表明这些颗粒对人体健康的危害最大。研究指出，人类许多疾病都和吸入颗粒物污染有着直接或间接的联系。因此，研究可吸入颗粒物在呼吸道中的运动特性，对于帮助了解可吸入颗粒的致病机理以及气溶胶治疗有着非常重要的意义，对于保护人类的健康，保护环境，提高生活质量具有积极意义。

传统技术中通过建立简化的气管支气管模型，借助计算流体动力学对其中的空气流动进行数值模拟，获得颗粒运动轨迹，根据运动轨迹确定颗粒沉积规律。然而，目前技术中的模型与人体相关性较小，获得的运动轨迹与实际情况相差较大。

发明内容

基于此，有必要针对获取颗粒运动轨迹准确率低的问题，提供一种肺泡内颗粒运动轨迹的检测方法。

一种肺泡内颗粒运动轨迹的检测方法，包括：

根据肺泡的几何数据对肺泡进行网格建模，建立肺泡动网格模型；

对所述肺泡动网格模型中进行空气流场模拟，获得空气流场的速度；

根据肺泡的空气流场的速度对肺泡的颗粒进行模拟，获得颗粒数值模拟模型；

根据颗粒数值模拟模型确定稳态和非稳态气流下不同受力条件下不同颗粒在肺泡内的运动轨迹。

上述肺泡内颗粒运动轨迹的检测方法，通过根据实际的肺泡几何数据对肺泡进行网格建模，大大提高了模型与人体肺泡的相似度。通过对所述肺泡动网格模型中进行空气流场模拟，获得空气流场的速度；根据肺泡的空气流场的速度对肺泡的颗粒进行模拟，获得颗粒数值模拟模型；从而可以根据颗粒数值模拟模型确定稳态和非稳态气流下、不同受力条件下、不同颗粒直径下颗粒在肺泡内的运动轨迹。提高了获取颗粒在肺泡内的运动轨迹的准确率，进而可以根据不同受力条件不同颗粒直径下颗粒的运动轨迹确定颗粒沉积规律，从而提高气溶胶的药性。

附图说明

图1为本发明肺泡内颗粒运动轨迹的检测方法实施例的流程示意图；

图2为本发明实施例中肺泡动网格模型的正面示意图；

图3为本发明实施例中肺泡动网格模型的状态示意图；

图4为呼吸过程中非稳态时气流流量示意图；

图5为肺泡囊的二维几何模型示意图；

图6为本发明应用实例中不同受力条件下0.1μm颗粒在肺泡模型内稳态气流下的运动轨迹示意图；

图7为本发明应用实例中不同受力条件下1.0μm颗粒在肺泡模型内稳态气流下的运动轨迹示意图；

图8为本发明应用实例中不同受力条件下10μm颗粒在肺泡模型内稳态气流下的运动轨迹示意图；

图9为本发明应用实例中在完整力作用下各个时间间隔中0.1μm颗粒所在位置示意图；

图10为本发明应用实例中在重力作用下各个时间间隔中0.1μm颗粒所在位置示意图；

图11为本发明应用实例中在失重作用下各个时间间隔中0.1μm颗粒所在位置示意图；

图12为本发明应用实例中在完整力作用下各个时间间隔中1.0μm颗粒所在位置示意图；

图13为本发明应用实例中在重力作用下各个时间间隔中1.0μm颗粒所在位置示意图；

图14为本发明应用实例中在失重作用下各个时间间隔中1.0μm颗粒所在位置示意图；

图15为本发明应用实例中在完整力作用下各个时间间隔中10μm颗粒所在位置示意图；

图16为本发明应用实例中在重力作用下各个时间间隔中10μm颗粒所在位置示意图；

图17为本发明应用实例中在失重作用下各个时间间隔中10μm颗粒所在位置示意图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细说明，但本发明的实施方式不限于此。

如图1所示，为本发明肺泡内颗粒运动轨迹的检测方法实施例的流程示意图，包括步骤：

步骤S101：根据肺泡的几何数据对肺泡进行网格建模，建立肺泡动网格模型；

肺泡的几何数据可以是指某一节肺泡的几何数据。所述几何数据可以包括：管长参数，内腔直径，肺泡直径，两肺泡间的中心距离，肺泡开口角，呼吸参数。比如，可以是第18节肺泡。几何数据可以包括：管长为740μm，内腔直径为320μm，肺泡直径为320μm，两肺泡间的中心距离为340μm，肺泡开口角为120°，呼吸参数为0.283。

本实施例还采用动网格模型模拟肺泡，提高模拟的真实性和相关性。肺泡动网格模型可以为3D模型。

步骤S102：对所述肺泡动网格模型中进行空气流场模拟，获得空气流场的速度；

进行空气流场模拟的方法有很多种，比如可以采用Fluent软件进行空气流场模拟。可以采用SIMPLE(Semi-Implicit Method for Pressure-Linked Equations)方法、SIMPLEC方法或PISO(Pressure-Implicit with Splitting of Operators)等方法计算空气流场的速度。

步骤S103：根据肺泡的空气流场的速度对肺泡的颗粒进行模拟，获得颗粒数值模拟模型；

颗粒数值模拟模型是一种数学模型，是可以用于计算颗粒在肺泡内运动轨迹的模型。

步骤S104：根据颗粒数值模拟模型确定稳态和非稳态气流下不同受力条件下不同颗粒在肺泡内的运动轨迹。

不同受力条件可以包括完整作用、失重作用和重力作用等情况。不同颗粒可以是指直径不同的颗粒。比如，直径分别为0.1μm，1.0μm和10μm的颗粒。在稳态气流测试时，可以将这些颗粒均匀的从肺泡的入口随气流一起射入肺泡内。其中，颗粒的入口速度可以与空气流场速度一致。

本实施例通过根据实际的肺泡几何数据对肺泡进行网格建模，大大提高了模型与人体肺泡的相似度。通过对所述肺泡动网格模型中进行空气流场模拟，获得空气流场的速度；根据肺泡的空气流场的速度对肺泡的颗粒进行模拟，获得颗粒数值模拟模型；从而可以根据颗粒数值模拟模型确定稳态和非稳态气流下不同受力条件下不同颗粒在肺泡内的运动轨迹。提高了获取颗粒在肺泡内运动轨迹的准确率，进而可以根据不同受力条件不同大小颗粒的运动轨迹确定颗粒沉积规律，从而提高气溶胶的药性。

在其中一个实施例中，所述颗粒数值模拟模型包括：

$\frac{d{X}_p}{\mathrm{dt}}=u_P,$

$\frac{{\mathrm{du}}_p}{\mathrm{dt}}=F_D(u-u_p)+\frac{g_i({\mathrm{\ρ}}_p-\mathrm{\ρ})}{{\mathrm{\ρ}}_p}+F_B,$

其中， $F_D=\frac{18\mathrm{\μ}}{d_p^2{\mathrm{\ρ}}_p{C}_C},C_C=1+\frac{2\mathrm{\λ}}{d_p}[1.257+0.4\exp (-1.1d_p/2\mathrm{\λ})]$

$F_B={\mathrm{\ζ}}_i\sqrt{\frac{\mathrm{\π}{S}_0}{\mathrm{\Δt}}},S_0=\frac{216v{k}_{B}T}{{\mathrm{\π}}^2\mathrm{\ρ}{d}_p^5{\left(\frac{{\mathrm{\ρ}}_p}{\mathrm{\ρ}}\right)}^2{C}_C}$

其中，X_p表示t时刻颗粒在空气流场中的位置，u_p表示颗粒的速度，u表示空气流场的速度，F_D表示作用在颗粒上的拽力，ρ_p表示颗粒的密度，ρ表示空气流场的密度，g_i表示重力矢量，F_B表示作用于颗粒上的布朗力，μ表示空气流场的粘度，d_p表示颗粒的直径，C_C表示滑移修正系数，λ表示分子的平均自由程，ζ_i表示一个具有零平均值和单位方差的高斯随机矢量，S₀表示白噪声过程的幅度，v表示空气的动力粘度，k_B＝1.38×10^-16，T表示空气流场的温度。

本实施例中的F_D充分考虑了微小颗粒下，分子间的滑移作用。通过本实施例中的颗粒数值模拟模型，可以获得各个时刻、不同颗粒大小、不同受力条件下颗粒的运动轨迹。

在其中一个实施例中，具体介绍了一种获得空气流场的速度的方法。即步骤S102，包括：

A1：设置肺泡动网格模型初始化的空气流场压强和质量通量；

在第一次计算中，变量由初始化过程更新获得。

A2：根据所述初始化的空气流场压强和质量通量求解空气流场的动量模型，并根据所述动量模型确定空气流场的第一速度；

动量模型可以根据动量守恒定律获得。根据动量模型即可获得空气流场速度。比如动量模型可以为：

$\frac{\∂\left(\mathrm{\ρ}{u}_i\right)}{\∂t}+\frac{\∂\left(\mathrm{\ρ}{u}_i{u}_j\right)}{\∂x_j}=-\frac{\∂P}{\∂x_i}+\frac{\∂\left({\mathrm{\τ}}_{\mathrm{ij}}\right)}{\∂x_j}$

ρ表示流体密度，u_i表示流体i方向速度，u_j表示流体j方向速度，x_i表示流体i方向位置，x_j表示流体j方向位置，其中，i为1或2或3，分别表示不同轴，j为1或2或3，分别表示不同轴，P表示流体微元上的压力，τ_ij表示粘性应力，

${\mathrm{\τ}}_{\mathrm{ij}}=\mathrm{\μ}(\frac{\∂u_i}{\∂x_j}+\frac{\∂u_j}{\∂x_i})-\frac{2\mathrm{\μ}{\mathrm{\δ}}_{\mathrm{ij}}\∂u_i}{3\∂x_j}$

μ表示流体的粘度(又称流体的运动粘度)，δ_ij表示流体的动力粘度。

A3：根据空气流场的连续方程获得修正模型，将所述第一速度、初始化的空气流场压强和质量通量代入修正模型进行修正，获得压强修正值、速度修正值和质量通量修正值；

修正模型可以是由连续方程导出的泊松型方程，用来对压强等参数进行修正，以便使修正后的数据满足连续方程。其中，连续方程可以为

A4：根据预设的收敛条件，判断所述压强修正值、速度修正值和质量通量修正值是否满足收敛条件，若是，根据压强修正值和质量通量修正值求解空气流场的动量模型，得到空气流场的第二速度，将第二速度设为空气流场的速度；若否，则将所述压强修正值、质量通量修正值替换初始化的空气流场压强和质量通量，重新计算压强修正值、速度修正值和质量通量修正值，直至满足收敛条件。

收敛条件可以是分别将当前压强修正值、速度修正值和质量通量修正值和修正前的压强、第一速度、质量通量进行误差比较，误差同时小于对应的设定值时，则表示满足收敛条件，否则，表示不满足。当收敛条件满足时，则可以获得空气流程的速度。当收敛条件不满足时，则将获得的压强修正值、质量通量修正值作为当前数据代替初始化的空气流场压强和质量通量，重新返回步骤A1进行计算，直到收敛条件满足。

建立肺泡动网格模型，还可以对动网格模型进行更新。更新的方法有很多种，比如弹性光滑法(Spring-based Smoothing)和局部网格重划法(LocalRemeshing Method)。弹性光滑法中又可以采用棱边弹簧。

在棱边弹簧中，平衡长度等于边的原长，网格点在初始状态所受到的合力为零。根据胡克定律，当弹簧系统中网格点移动后，作用在其中某个网格节点i上的合力为

${\stackrel{\→}{F}}_i=\underset{j}{\overset{n_i}{\mathrm{\Σ}}}{k}_{\mathrm{ij}}(\mathrm{\Δ}{\stackrel{\→}{x}}_j-\mathrm{\Δ}{\stackrel{\→}{x}}_i)$

其中和分别为节点i以及和i的相邻节点j的矢量位移。n_i为与节点i相连的节点数量，k_ij是节点i、j之间弹簧的弹性系数。

在弹簧系统达到平衡状态时，每个节点上面受到的合力为零。通过对系统节点的平衡方程组进行迭代求解：

$\mathrm{\Δ}{\stackrel{\→}{x}}_i^{m+1}=\frac{{\mathrm{\Σ}}_i^{n_i}{k}_{\mathrm{ij}}\mathrm{\Δ}{\stackrel{\→}{x}}_j^m}{{\mathrm{\Σ}}_j^{n_i}{k}_{\mathrm{ij}}}$

其中边界上节点的位移已知，在边界上的节点更新后，可以对上式进行迭代求解。上式中，是节点j在第m次迭代时的位移值，计算收敛后就能得到网格系统中每个节点运动后的位移矢量。节点i在网格更新后的位置：

${\stackrel{\→}{x}}_i^{n+1}={\stackrel{\→}{x}}_i^n+\mathrm{\Δ}{\stackrel{\→}{x}}_i^{\mathrm{converged}}$

在顶点弹簧描述中弹簧的平衡长度为零。为了使网格在变形以后仍能保持合理的疏密分布，弹簧的倔强系数可以有多种选择，其中可以假设弹簧的倔强系数和两节点之间边长的某次幂成反比。

本实施例中棱边的弹性系数取值为：

$k_{\mathrm{ij}}=\frac{1}{\sqrt{{\stackrel{\→}{x}}_j-{\stackrel{\→}{x}}_i}}$

二维三角形网格和三维四面体网格可以使用弹性光顺法。

局部网格重划法，在网格系统用三角形或四面体网格组成时，如果边界的移动和变形远远大于网格尺寸，可能导致局部网格发生严重畸变，甚至出现体积为负的情况，或者网格的畸变过大使计算无法收敛。处理方法是将畸变率过大或尺寸变化过于剧烈的网格集中在一起进行局部网格的重新划分，如果重新划分后的网格可以满足畸变率和尺寸的要求，则用新的网格代替原来的网格，如果新的网格仍然无法满足要求，则放弃重新划分的结果。在局部网格重划法只能应用于二维三角形和三维四面体网格单元。上述两种方法可以在Fluent中实现。

本发明方案还提供一种优选的更新方法，该方法适用于六面体网格、楔形网格等可以在边界上分层的网格系统。具体包括：

根据肺泡的几何数据对肺泡进行网格建模后，检测相邻边界网格的高度与预设优化高度值的差值，当差值大于第一设定值时，在该相邻边界网格和边界层之间增加一层网格，当差值小于第二设定值时，删除该相邻边界网格，其中，所述相邻边界网格为与边界层相邻的一层网格。

本方法是在边界上假定一个优化的网格层高度值，在边界发生移动、变形时，如果临近边界的一层网格的高度同优化高度值相比大到一定程度时，就在边界面与相邻网格层之间增加一层网格。反之，如果临近边界的网格被压缩到一定程度时，临近一层网格又会被删除，用这种办法使边界上的网格保持在一定的密度。

在其中一个实施例中，本发明方案还可以采用DEFINE_GRID_MOTION宏控制动网格的自动更新。通过DEFINE_GRID_MOTION宏来定义呼吸过程中肺泡壁面随时间的扩张和收缩运动，其基本格式可以为：

DEFINE_GRID_MOTION(name，d，dt，time，dtime)，

name为用户指定的名称，d和dt为两个分别指向计算区域和动网格区域的指针，time和dtime分别为流场的当前时间和时间步长。

所述肺泡动网格模型为3D模型。在其中一个实施例中，所述肺泡动网格模型中网格可以是四面体、六面体、楔形体等中的一种或多种。其中，肺泡动网格模型的网格类型可以为混合网格，网格数可以为31万左右。如图2和图3所示，为本发明实施例中肺泡动网格模型的正面示意图和状态示意图。

数值模拟采用三维稳态和非稳态计算，对气流的计算采用SIMPLE方法求解N-S方程。计算的气流采用层流(laminar)计算模型，同时可以利用Fluent中的DPM模型对颗粒运动轨迹进行追踪。气流的入口速度按成人正常呼吸时，比如，第18节肺泡的流量给定。成人正常坐立状态下的第18节的呼吸参数：Re＝0.283,空气的参数为：ρ＝1.225kg/m³，μ＝1.7894e-05kg/m.s,ν＝μ/ρ＝1.46e-05m²/s，Re＝UD/ν，U＝0.0129m/s。非稳态时，肺泡内的非稳态气流流量随时间成正弦曲线分布，呼吸周期T为4秒，计算为时间步长0.02秒，每个时间间隔迭代次数以残差值趋向稳定为限， $Q_{\max }=\stackrel{\‾}{u}\×A\×\mathrm{\π}/2,\stackrel{\‾}{u}=0.129m/s,$ 求得 $u_{\max }=\stackrel{\‾}{u}\×\mathrm{\π}/2.$ 如图4所示，为呼吸过程中非稳态时气流流量示意图。

计算中，计算了一个完整周期的气流在肺泡内的流动的3D模型的稳态和非稳态过程，对颗粒的计算以颗粒完全逃离肺泡或者沉积为结束，最大模拟时间为T。颗粒的沉积率(DE)的计算方法为：DE＝Nd/Nt100％.Nd表示颗粒在壁面的沉积数，Nt表示总颗粒数。

在其中一个实施例中，在建立肺泡动网格模型过程中，还包括：建立肺泡囊的二维几何模型，用于模拟呼吸过程中的肺泡的扩张和收缩，其中，将前后开口的圆柱管形模拟呼吸道中的气管，将连接在圆柱壁面上的圆形模拟呼吸道中的肺泡，肺泡位于该段气管的正中位置。还可以设置相关参数，比如气管的半径为250um，肺泡的半径为200um，肺泡在主气管上的开口角设为60度，气管的长度为1000um。

如图5所示，为肺泡囊的二维几何模型。本实施例可以利用参数模型来模拟呼吸过程中肺泡的扩张和收缩。使用前后开口的圆柱管模拟呼吸道中的气管，使用连接在圆柱壁面上球形盖来描述呼吸道中的肺泡。本实施例中R_a和R_d分别为肺泡和气管的半径，γ为肺泡的开口角。结合参数模型中的数据，本实施例模型中设气管的半径R_d为250um，肺泡的半径为200um，肺泡在主气管上的开口角设为60度。气管的长度为1000um，并且在本文的计算模型中，肺泡位于该段气管的正中位置。

在其中一个实施例中，可以将人体肺泡和气管在呼吸过程中的变形过程视为一个保持几何形状相似的扩张和收缩运动。即在建立肺泡囊的二维几何模型过程中，还包括设置肺泡和气管的各个特征尺寸为随时间变化的正弦函数：

$L\left(t\right)=L_0[1+\frac{\mathrm{\β}}{2}+\frac{\mathrm{\β}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\π}}{2})]$

其中，L(t)表示肺泡或气管的特征尺寸的值，L₀表示几何模型中的各个特征尺寸在吸气开始时或呼气结束时的值，f表示呼吸的频率，f＝2π/T，T表示呼吸的周期时间，β表示各个特征尺寸的扩张幅度，β＝(C+1)^1/3-1，C＝(V_max-V_min)/V_min，V_max和V_min表示肺泡动网格模型的最大体积和最小体积，所述最大体积和最小体积分别对应在实际呼吸过程中呼气终了和吸气终了时的体积。

在t＝T/2时刻(吸气结束时或者呼气开始时)，模型中的特征长度达到最大值：

L_max＝L₀(1+β)

在呼吸过程中，每一级模型上的当地Re数和气管中的气流速度均为时间的函数：

在其中一个实施例中，在建立肺泡囊的二维几何模型过程中，还包括：根据呼吸过程中进出肺泡的气体流量与流过肺泡所处气管的气体流量的比值确定肺泡在呼吸道上所处位置，根据所述位置建立肺泡囊的二维几何模型，其中：

$Q_a\left(t\right)=\frac{27}{16}\mathrm{\π\βf}{R}_a^2{\mathrm{\λ}}^2\cos (\mathrm{ft}-\frac{\mathrm{\π}}{2}),\mathrm{\λ}=1+\frac{\mathrm{\β}}{2}+\frac{\mathrm{\β}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\π}}{2})$

$Q_d\left(t\right)=\frac{{\mathrm{dV}}_t}{\mathrm{dt}}$

Q_a(t)表示呼吸过程中进出肺泡的气体流量，R_a表示肺泡半径，Q_d(t)表示流过肺泡所处气管的气体流量，V_t表示将经过该段气管后的所有气管及气管上所有肺泡的体积求和获得的值。

在本实施例中通过在每个几何模型中取不同的Q_a/Q_d值，来表示该模型在整个肺泡区气管树上所处的位置。在人体的呼吸过程中，最大Re数发生在肺泡区第一级上、呼气或吸气的峰值时刻，计算表明这个最大Re的值大约在20左右，大多数情况下流动的Re数小于1。在这样的流场条件下，可以忽略进口处速度分布对整个流场的影响。对于低雷诺数的圆管内层流而言，从截面上的均匀速度分布发展成为抛物线型速度分布，需要的管长L_e可以由根据以下的表达式得到：

$\frac{L_e}{D}=0.06\mathrm{Re}$

其中，D为管道的直径。相对于本实施例模型中肺泡前的气管长度，这段长度是可以不予不计的。因此，在本实施例中设定一个在进口表面上均匀分布的速度是比较合理的，到达肺泡处时气管内的流动已经充分发展。

在本实施例的模型中，肺泡的边界的变形量大于位于边界上网格的高度，可以使用弹性光顺法和局部重划法来控制模型中网格的变形运动，因此对于本实施例的计算模型的网格划分可以使用三角形网格。

本方案基于上述检测方案，举其中一种具体运用实例进行说明。

本实施例中的空气流场作为层流流动。为了获得空气流场速度和颗粒运动轨迹，可以采用模拟方法实现。

B1：流场变量更新。在第一次计算时，变量由初始化过程更新。在随后的计算中，每迭代一次即得到一个更新的解。

B2：用当前压强和质量通量的值求解动量模型(动量方程)，以得到新的速度场。

B3：求解压强修正模型。压强修正模型是由连续方程导出的泊松型方程，求解这个方程可以得到对压强场、速度场和质量通量的修正，进而使连续方程得到满足。

B4：利用前面求出的解，求解动量模型、能量模型。

B5：在多相流(空气流场和颗粒)计算中，当考虑相间干扰时，则需要通过求解弥散相轨迹计算得到连续相方程中的源项解。

B6：检验收敛条件是否被满足。如果收敛条件被满足，则停止计算。如果计算没有收敛，则继续迭代过程。

在分离求解算法中，由于连续性方程和动量方程不需要耦合，由动量方程计算出的流场速度值可能不能同时满足连续性方程，则用一个由连续性方程和线性化后的动量方程推导得到的“泊松类”方程来对压力进行修正，以得到满足连续性方程的压力和速度。

在求解低速不可压缩问题时，速度、压力、温度之间的耦合程度不是很高，采用上述方法具有很高的数值稳定性。

其中，流动的控制方程包括：

根据质量守恒，且流体作为不可压缩流体，可以将空气密度视为常数，即得到连续方程：

$\frac{\∂u_i}{\∂x_i}=0$

根据动量守恒定律得到动量模型：

$\frac{\∂\left(\mathrm{\ρ}{u}_i\right)}{\∂t}+\frac{\∂\left(\mathrm{\ρ}{u}_i{u}_j\right)}{\∂x_j}=-\frac{\∂P}{\∂x_i}+\frac{\∂\left({\mathrm{\τ}}_{\mathrm{ij}}\right)}{\∂x_j}$

ρ表示流体密度，u_i表示流体i方向速度，u_j表示流体j方向速度，x_i表示流体i方向位置，x_j表示流体j方向位置，其中，i为1或2或3，分别表示各轴，j为1或2或3，分别表示各轴，P表示流体微元上的压力，τ_ij表示粘性应力，

${\mathrm{\τ}}_{\mathrm{ij}}=\mathrm{\μ}(\frac{\∂u_i}{\∂x_j}+\frac{\∂u_j}{\∂x_i})-\frac{2\mathrm{\μ}{\mathrm{\δ}}_{\mathrm{ij}}\∂u_i}{3\∂x_j}$

μ表示流体的运动粘度，δ_ij表示流体的动力粘度。

根据能量守恒定律获得能量模型：

$\frac{\∂\left(\mathrm{\ρT}\right)}{\∂t}+\frac{\∂\left(\mathrm{\ρ}{u}_{j}T\right)}{\∂x_j}=\frac{\∂\left(\frac{k\∂T}{c_P\∂x_j}\right)}{\∂x_j}+S_T$

其中，c_P为比热容，T为温度，k为流体的传热系数，S_T为粘性耗散项。

把气相作为连续相考虑，通过Navier-Stokes方程进行求解。通过计算得到的连续相流场来跟踪流场中一定数量的颗粒、气泡或液滴的运动，计算时可以考虑颗粒相同气相之间发生的动量、质量和能量交换。该模型应用在气流场中离散颗粒相应足够稀疏的情况下，离散相只占据很小的体积分数，这样颗粒与颗粒之间的作用、颗粒的体积分数对连续相的作用都可以忽略不计。则颗粒数值模拟模型包括：

$\frac{d{X}_p}{\mathrm{dt}}=u_P,$

$\frac{{\mathrm{du}}_p}{\mathrm{dt}}=F_D(u-u_p)+\frac{g_i({\mathrm{\ρ}}_p-\mathrm{\ρ})}{{\mathrm{\ρ}}_p}+F_B.$

基于上述模拟方法，进行了稳态气流下可吸入颗粒在在肺泡内运动特性的3D数值模拟和非稳态气流下可吸入颗粒在在肺泡内运动特性的3D数值模拟。

首先，进行稳态气流下可吸入颗粒在在肺泡内运动特性的3D数值模拟。

在3D肺泡动网格模型内稳态气流下，分别选用了直径为0.1-，1.0-和10μm的颗粒均匀的从肺泡的入口随气流一起射入肺泡内，颗粒形状为球形，密度为1000kg/m³。颗粒的入口速度和气流一致，即0.0129m/s，颗粒的分布方式可以选用的是FLUENT里面DMP模型中的surface分布方式，颗粒数和网格数量有关，本计算中颗粒数量可以为852。

计算中，模拟计算在完整力作用，失重作用和重力作用三种情况下的可吸入颗粒运动沉积特性，通过对三种受力条件下颗粒的运动轨迹的变化以及沉积效率的改变来评估惯性，重力和布朗运动对颗粒运动特性的影响。

首先，进行稳态气流下0.1μm颗粒在肺泡内运动特性3D数值模拟结果与分析。如图6所示，为本发明应用实例中不同受力条件下0.1μm颗粒在肺泡模型内稳态气流下的运动轨迹示意图。其中，6a的受力条件为完整力作用，6b的受力条件为失重作用，6c的受力条件为重力作用。

检测0.1μm颗粒在不同受力条件下的运动轨迹可以得出，三种受力条件下，颗粒的运动轨迹几乎是一样的，三种颗粒完全随主气流运动，运动轨迹和主气流流线重合，颗粒并没有进入肺泡的回流区。由于采用的是层流计算模型，使得稳态下气流没有扰动，主气流和回流之间存在明显的分隔面，同时，重力、布朗运动力都不能对颗粒产生足够使颗粒发生扰动的作用力，最终导致了三种颗粒的沉积效率为0。综上可得到，对于0.1μm的颗粒在肺泡内的运动，重力和布朗运动力都不是颗粒收到的主要运动力，颗粒主要受惯性作用。颗粒的沉积主要是由扰动气流产生的在泡内的沉积，以及在各节之间的分叉处的碰撞沉积。由于缺少扰动和分叉，导致了沉积效率过低。

同时，在肺泡内，气流强度是非常小的，气流扰动也不可能很大，同时颗粒主要受惯性作用，因此对于0.1μm颗粒的沉积效率是很低的。在本模型中，0.1μm颗粒在肺泡内稳态气流下的运动主要受惯性作用，颗粒沉积效率几乎为0。

接着，进行稳态气流下1.0μm颗粒在肺泡内运动特性3D数值模拟结果与分析。如图7所示，为本发明应用实例中不同受力条件下1.0μm颗粒在肺泡模型内稳态气流下的运动轨迹示意图。其中，7a的受力条件为完整力作用，7b的受力条件为失重作用，7c的受力条件为重力作用。

从图7中可以看出，颗粒在肺泡的轨迹结果和0.1μm颗粒的结果一致，即三种受力条件下表现出来的结果是一致的，颗粒的沉积仍然是0。可以看出，对于1.0μm的颗粒在肺泡内的运动，惯性力仍然是主要作用力，重力和布朗运动仍然很不明显。

由于采用的是稳态层流流动计算模型，气流的扰动不存在，使得颗粒只能随主气流流动，无法进入肺腺泡的回流区。因此，重力和布朗运动不能够产生这种作用(主要是产生足够的扰动，导致颗粒能够越过分隔面，进行对流交换)，即这两种力对0.1-1.0μm的颗粒的作用是不明显的。

接着，进行稳态气流下10μm颗粒在肺泡内运动特性3D数值模拟结果与分析。如图8所示，为本发明应用实例中不同受力条件下10μm颗粒在肺泡模型内稳态气流下的运动轨迹示意图。其中，8a的受力条件为完整力作用，8b的受力条件为失重作用，8c的受力条件为重力作用。

从图8可以看出，在完整力和重力条件下，颗粒出现了大量沉积，沉积率为62.8％(535/852，536/852)。在失重条件下，颗粒轨迹和主气流重合。表明重力作用对10μm颗粒的运动产生了很重要的作用，而布朗运动则对10μm颗粒运动的作用不明显。从气流场中知道，回流区和主气流之间存在分隔面，两者之间不存在对流交换，而重力作用使得颗粒能够穿过分隔面，同时在腺泡内，颗粒也不随回流运动，而主要受重力作用，说明对于10μm颗粒，颗粒主要受到重力作用。同时，颗粒在失重条件下，主要随气流运动，说明惯性仍然是颗粒运动的主要作用，但对于颗粒的沉积，重力是主要作用力。

上述模拟中，利用肺泡模型模拟计算了稳态气流下肺泡内的气流特征以及稳态气流下不同粒径的颗粒(0.1-，1.0-，10μm)在不同的受力条件下的运动沉积特性。从模拟计算的结果可以看出，0.1-和1.0μm在肺泡内的运动和沉积特性相似度很高，颗粒运动主要受到惯性作用，颗粒沉积非常的少；10μm颗粒则受到重力作用很大，颗粒出现的较大量的沉积。同时，总体来说，稳态气流下可吸入颗粒在肺泡内的运动特性的3D数值模拟结果可以归纳为：

微细颗粒(比如0.1-，1.0μm)在肺泡内的颗粒沉积特性是相同的，主要受到惯性力作用，颗粒沉积非常的少。粗颗粒(比如10μm)在肺泡内的沉积主要受到重力作用，颗粒出现大量沉积；采用层流计算模型导致气流扰动的缺少，是导致微细颗粒布朗运动不明显以及颗粒沉积效率偏低的原因之一。

接着进行非稳态气流下可吸入颗粒在在肺泡内运动特性的3D数值模拟。在3D非稳态气流下模拟颗粒在肺泡内的运动沉积，在计算中，为了节约时间，提高效率，只在初始时刻以速度0射入不同直径的颗粒。颗粒在肺泡内的运动特性计算采用FLUENT中的DPM模型，计算中总共射入的颗粒数为852个，同时，计算中对每一个时间间隔(时间间隔可以设为0.02s)的颗粒特性都进行统计，直到肺泡内没有颗粒存在为止。在3D非稳态模拟计算中，同样选择不同粒径的颗粒在不同的受力条件下来模拟颗粒沉积，来评估计算颗粒直径以及不同的受力对颗粒的运动沉积的作用。

首先，进行非稳态气流下0.1μm颗粒在肺泡内运动特性3D数值模拟结果与分析。

在计算模型中，0.16秒后0.1μm就全部离开了肺泡出口，并且在三种不同的受力条件下，0.1μm颗粒均没有在肺泡内沉积，颗粒的沉积率均为0。如图9，为本发明应用实例中在完整力作用下各个时间间隔中0.1μm颗粒所在位置示意图，9a至9f分别为0.02s时0.1μm颗粒所在位置示意图、0.04s时0.1μm颗粒所在位置示意图、0.06s时0.1μm颗粒所在位置示意图、0.08s时0.1μm颗粒所在位置示意图、0.10s时0.1μm颗粒所在位置示意图、0.12s时0.1μm颗粒所在位置示意图。如图10，为本发明应用实例中在重力作用下各个时间间隔中0.1μm颗粒所在位置示意图，10a至10f分别为0.02s时0.1μm颗粒所在位置示意图、0.04s时0.1μm颗粒所在位置示意图、0.06s时0.1μm颗粒所在位置示意图、0.08s时0.1μm颗粒所在位置示意图、0.10s时0.1μm颗粒所在位置示意图、0.12s时0.1μm颗粒所在位置示意图。如图11，为本发明应用实例中在失重作用下各个时间间隔中0.1μm颗粒所在位置示意图，11a至11f分别为0.02s时0.1μm颗粒所在位置示意图、0.04s时0.1μm颗粒所在位置示意图、0.06s时0.1μm颗粒所在位置示意图、0.08s时0.1μm颗粒所在位置示意图、0.10s时0.1μm颗粒所在位置示意图、0.12s时颗粒所在位置示意图。

从图中可以看出，颗粒主要随主气流运动，在肺泡入口处，颗粒随主气流有部分的绕入到肺泡内，但是最终仍然是随主气流绕出肺泡，并没有进入腺泡的回流区。这表明，0.1μm颗粒的在肺泡内的运动没有在分隔面上发生物质的对流交换。同时，颗粒的运动主要随主气流运动，重力和布朗运动不明显，主要是惯性作用对颗粒产生影响。从颗粒进入到离开肺泡的时间可以知道只有0.16秒，肺泡内主气流速度过快，导致布朗运动不能够及时产生作用很可能是布朗运动表现不明显的一个重要原因。三种受力条件下，颗粒的运动轨迹几乎是一致的，说明重力对0.1颗粒的运动产生的作用不大，表现不明显。因此，可以指出，所采用的模型中，0.1μm颗粒在肺泡内的运动主要受惯性作用，布朗运动和重力作用不明显。

接着，进行非稳态气流下1.0μm颗粒在肺泡内运动特性3D数值模拟结果与分析。

对1.0μm颗粒在肺泡模型内非稳态气流下的运动特性进行了模拟计算。如图12，为本发明应用实例中在完整力作用下各个时间间隔中1.0μm颗粒所在位置示意图，12a至12f分别为0.02s时1.0μm颗粒所在位置示意图、0.04s时1.0μm颗粒所在位置示意图、0.06s时1.0μm颗粒所在位置示意图、0.08s时1.0μm颗粒所在位置示意图、0.10s时1.0μm颗粒所在位置示意图、0.12s时1.0μm颗粒所在位置示意图。如图13，为本发明应用实例中在重力作用下各个时间间隔中1.0μm颗粒所在位置示意图，13a至13f分别为0.02s时1.0μm颗粒所在位置示意图、0.04s时1.0μm颗粒所在位置示意图、0.06s时1.0μm颗粒所在位置示意图、0.08s时1.0μm颗粒所在位置示意图、0.10s时1.0μm颗粒所在位置示意图、0.12s时1.0μm颗粒所在位置示意图。如图14，为本发明应用实例中在失重作用下各个时间间隔中1.0μm颗粒所在位置示意图，14a至14f分别为0.02s时1.0μm颗粒所在位置示意图、0.04s时1.0μm颗粒所在位置示意图、0.06s时1.0μm颗粒所在位置示意图、0.08s时1.0μm颗粒所在位置示意图、0.10s时1.0μm颗粒所在位置示意图、0.12s时1.0μm颗粒所在位置示意图。

从检测结果可以看出，1.0μm颗粒和0.1μm颗粒在肺泡内的运动轨迹几乎是一样的。1.0μm颗粒在肺泡的运动也主要是随主气流运动的，颗粒从入口到出口所用的时间为0.16秒，颗粒在三种受力情况下的沉积率均为0。通过三种受力下的各个时间间隔的颗粒位置进行比较，可以得出，三种受力下，颗粒的运动轨迹是一样的，颗粒主要随主气流运动，在腺泡附近随气流绕入并绕出，非稳态气流的扰动不足以使颗粒穿过回流去和主气流之间的分隔面使得颗粒从主气流中偏离到回流区中。同时，对于1.0μm颗粒，重力和布朗运动对颗粒的沉积产生的影响并不重要，1.0μm颗粒在肺泡内的运动主要是受惯性作用。可见，1.0μm颗粒在肺泡内的运动特性和0.1μm颗粒在肺泡内的运动特性是非常的相似。

接着，进行非稳态气流下10μm颗粒在肺泡内运动特性3D数值模拟结果与分析。

本检测试验中，10μm颗粒在肺泡内的运动沉积特性得到了模拟计算。如图15，为本发明应用实例中在完整力作用下各个时间间隔中10μm颗粒所在位置示意图，15a至15f分别为0.02s时10μm颗粒所在位置示意图、0.04s时10μm颗粒所在位置示意图、0.06s时10μm颗粒所在位置示意图、0.08s时10μm颗粒所在位置示意图、0.10s时10μm颗粒所在位置示意图、0.12s时10μm颗粒所在位置示意图。如图16，为本发明应用实例中在重力作用下各个时间间隔中10μm颗粒所在位置示意图，16a至16f分别为0.02s时10μm颗粒所在位置示意图、0.04s时10μm颗粒所在位置示意图、0.06s时10μm颗粒所在位置示意图、0.08s时10μm颗粒所在位置示意图、0.10s时10μm颗粒所在位置示意图、0.12s时10μm颗粒所在位置示意图。如图17，为本发明应用实例中在失重作用下各个时间间隔中10μm颗粒所在位置示意图，17a至17f分别为0.02s时10μm颗粒所在位置示意图、0.04s时10μm颗粒所在位置示意图、0.06s时10μm颗粒所在位置示意图、0.08s时10μm颗粒所在位置示意图、0.10s时10μm颗粒所在位置示意图、0.12s时10μm颗粒所在位置示意图。

从计算得到的结果可以看出，10μm颗粒从入口到出口总共用的时间为0.12秒，在重力和完整力作用下，出现了大量的沉积(68.2％)，失重条件下，颗粒主要随主气流运动，没有出现沉积现象。同时从颗粒沉积的位置可以看出，颗粒主要是在腺泡内部的沉积。同时，进入腺泡内的颗粒并没有随着回流运动，而主要是在重力作用下运动。这表明，对于10μm颗粒，在肺泡内的沉积主要受到重力作用，重力作用能够使颗粒穿越主气流和腺泡内回流区之间的分隔面，而进入腺泡，并沉积。因此，10μm颗粒在肺泡内的运动主要是受惯性和重力作用，同时，重力是颗粒沉积的主要原因。

上述检测中，模拟计算了非稳态气流下3D肺泡模型内的气流特征以及0.1-，1.0-和10μm非稳态气流下在肺泡内的运动特性以及颗粒的沉积效率。通过计算，可以得到了特征时刻下颗粒在肺泡内的位置和沉积状况。可以看出，颗粒在肺泡内的沉积和稳态气流下存在很大的相同性，0.1-和1.0μm在肺泡内的运动和颗粒沉积效率很低；10μm颗粒沉积效率较大。即可以得出肺泡内非稳态气流下颗粒运动的以下几个特点：

微细颗粒(比如0.1-1.0μm)在3D模型中的运动和沉积主要受惯性作用，粗颗粒(比如10μm)主要受重力作用。微细颗粒(0.1-1.0μm)在3D模型中效率最低(0)，粗颗粒在模型中的沉积最高，达到68.2％。从而可以根据这些规律研究生产气溶胶，提高气溶胶的药性。

以上实施方式中的各种技术特征可以任意进行组合，只要特征之间的组合不存在冲突或矛盾，但是限于篇幅，未进行一一描述，因此上述实施方式中的各种技术特征的任意进行组合也属于本说明书公开的范围。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种肺泡内颗粒运动轨迹的检测方法，其特征在于，包括：

根据肺泡的几何数据对肺泡进行网格建模，建立肺泡动网格模型；

对所述肺泡动网格模型中进行空气流场模拟，获得空气流场的速度；

根据肺泡的空气流场的速度对肺泡的颗粒进行模拟，获得颗粒数值模拟模型；

根据颗粒数值模拟模型确定稳态和非稳态气流下不同受力条件下不同颗粒在肺泡内的运动轨迹。

2.根据权利要求1所述的肺泡内颗粒运动轨迹的检测方法，其特征在于，所述颗粒数值模拟模型包括：

$\frac{d{X}_p}{\mathrm{dt}}=u_P,$

$\frac{{\mathrm{du}}_p}{\mathrm{dt}}=F_D(u-u_p)+\frac{g_i({\mathrm{\ρ}}_p-\mathrm{\ρ})}{{\mathrm{\ρ}}_p}+F_B,$

其中， $F_D=\frac{18\mathrm{\μ}}{d_p^2{\mathrm{\ρ}}_p{C}_C},C_C=1+\frac{2\mathrm{\λ}}{d_p}[1.257+0.4\exp (-1.1d_p/2\mathrm{\λ})]$

$F_B={\mathrm{\ζ}}_i\sqrt{\frac{\mathrm{\π}{S}_0}{\mathrm{\Δt}}},S_0=\frac{216v{k}_{B}T}{{\mathrm{\π}}^2\mathrm{\ρ}{d}_p^5{\left(\frac{{\mathrm{\ρ}}_p}{\mathrm{\ρ}}\right)}^2{C}_C}$

其中，X_p表示t时刻颗粒在空气流场中的位置，u_p表示颗粒的速度，u表示空气流场的速度，F_D表示作用在颗粒上的拽力，ρ_p表示颗粒的密度，ρ表示空气流场的密度，g_i表示重力矢量，F_B表示作用于颗粒上的布朗力，μ表示空气流场的粘度，d_p表示颗粒的直径，C_C表示滑移修正系数，λ表示分子的平均自由程，ζ_i表示一个具有零平均值和单位方差的高斯随机矢量，S₀表示白噪声过程的幅度，v表示空气的动力粘度，k_B＝1.38×10^-16，T表示空气流场的温度。

3.根据权利要求1所述的肺泡内颗粒运动轨迹的检测方法，其特征在于，所述对所述肺泡动网格模型中进行空气流场模拟，获得空气流场的速度步骤，包括：

设置肺泡动网格模型初始化的空气流场压强和质量通量；

根据所述初始化的空气流场压强和质量通量求解空气流场的动量模型，并根据所述动量模型确定空气流场的第一速度；

根据空气流场的连续方程获得修正模型，将所述第一速度、初始化的空气流场压强和质量通量代入修正模型进行修正，获得压强修正值、速度修正值和质量通量修正值；

根据预设的收敛条件，判断所述压强修正值、速度修正值和质量通量修正值是否满足收敛条件，若是，根据压强修正值和质量通量修正值求解空气流场的动量模型，得到空气流场的第二速度，将第二速度设为空气流场的速度；若否，则将所述压强修正值、质量通量修正值替换初始化的空气流场压强和质量通量，重新计算压强修正值、速度修正值和质量通量修正值，直至满足收敛条件。

4.根据权利要求1所述的肺泡内颗粒运动轨迹的检测方法，其特征在于，所述几何数据包括：管长参数，内腔直径，肺泡直径，两肺泡间的中心距离，肺泡开口角，呼吸参数。

5.根据权利要求1所述的肺泡内颗粒运动轨迹的检测方法，其特征在于，建立肺泡动网格模型时，还包括对动网格模型进行更新，包括：

根据肺泡的几何数据对肺泡进行网格建模后，检测相邻边界网格的高度与预设优化高度值的差值，当差值大于第一设定值时，在该相邻边界网格和边界层之间增加一层网格，当差值小于第二设定值时，删除该相邻边界网格，其中，所述相邻边界网格为与边界层相邻的一层网格。

6.根据权利要求1所述的肺泡内颗粒运动轨迹的检测方法，其特征在于，所述肺泡动网格模型中网格包括四面体、六面体、楔形体中的一种或多种。

7.根据权利要求1所述的肺泡内颗粒运动轨迹的检测方法，其特征在于，所述肺泡动网格模型的网格类型为混合网格，网格数为31万。

8.根据权利要求1所述的肺泡内颗粒运动轨迹的检测方法，其特征在于，在建立肺泡动网格模型过程中，还包括：

建立肺泡囊的二维几何模型，用于模拟呼吸过程中的肺泡的扩张和收缩，其中，将前后开口的圆柱管形模拟呼吸道中的气管，将连接在圆柱壁面上的圆形模拟呼吸道中的肺泡，肺泡位于该段气管的正中位置。

9.根据权利要求8所述的肺泡内颗粒运动轨迹的检测方法，其特征在于，在建立肺泡囊的二维几何模型过程中，还包括设置肺泡和气管的各个特征尺寸为随时间变化的正弦函数：

$L\left(t\right)=L_0[1+\frac{\mathrm{\β}}{2}+\frac{\mathrm{\β}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\π}}{2})]$

其中，L(t)表示肺泡或气管的特征尺寸的值，L₀表示几何模型中的各个特征尺寸在吸气开始时或呼气结束时的值，f表示呼吸的频率，f＝2π/T，T表示呼吸的周期时间，β表示各个特征尺寸的扩张幅度，β＝(C+1)^1/3-1，C＝(V_max-V_min)/V_min，V_max和V_min表示肺泡动网格模型的最大体积和最小体积，所述最大体积和最小体积分别对应在实际呼吸过程中呼气终了和吸气终了时的体积。

10.根据权利要求9所述的肺泡内颗粒运动轨迹的检测方法，其特征在于，在建立肺泡囊的二维几何模型过程中，还包括：根据呼吸过程中进出肺泡的气体流量与流过肺泡所处气管的气体流量的比值确定肺泡在呼吸道上所处位置，根据所述位置建立肺泡囊的二维几何模型，其中：

$Q_a\left(t\right)=\frac{27}{16}\mathrm{\π\βf}{R}_a^2{\mathrm{\λ}}^2\cos (\mathrm{ft}-\frac{\mathrm{\π}}{2}),\mathrm{\λ}=1+\frac{\mathrm{\β}}{2}+\frac{\mathrm{\β}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\π}}{2})$

$Q_d\left(t\right)=\frac{{\mathrm{dV}}_t}{\mathrm{dt}}$

Q_a(t)表示呼吸过程中进出肺泡的气体流量，R_a表示肺泡半径，Q_d(t)表示流过肺泡所处气管的气体流量，V_t表示将经过该段气管后的所有气管及气管上所有肺泡的体积求和获得的值。