CN104027114B - 肺泡收缩和扩展过程的流场数值模拟测量方法和系统 - Google Patents

Abstract

本发明提供一种肺泡收缩和扩展过程的流场数值模拟测量方法，包括如下步骤：根据预设的肺泡的几何参数建立肺泡的几何模型，对所述肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型；获取所述肺泡对应的呼吸参数及流体参数；根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型；根据所述流场模型数值模拟测量所述肺泡在收缩和扩展过程时的流场特征。本发明还提供对应的系统，能缩短测量周期，减少运算负担，准确地描述肺泡内的流场特征。

Description

肺泡收缩和扩展过程的流场数值模拟测量方法和系统

技术领域

本发明涉及肺泡模拟检测技术领域，特别是涉及一种肺泡收缩和扩展过程的流场数值模拟测量方法，以及一种肺泡收缩和扩展过程的流场数值模拟测量系统。

背景技术

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

发明内容

基于此，本发明提供一种肺泡收缩和扩展过程的流场数值模拟测量方法和系统，能缩短测量周期，减少运算负担，准确地描述肺泡内的流场特征。

一种肺泡收缩和扩展过程的流场数值模拟测量方法，包括如下步骤：

根据预设的肺泡的几何参数建立肺泡的几何模型，对所述肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型；

获取所述肺泡对应的呼吸参数及流体参数；

根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型；

根据所述流场模型数值模拟测量所述肺泡在收缩和扩展过程时的流场特征。

一种肺泡收缩和扩展过程的流场数值模拟测量系统，包括：

建立模块，用于根据预设的肺泡的几何参数建立肺泡的几何模型，对所述肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型；

获取模块，用于获取所述肺泡对应的呼吸参数及流体参数；

模拟模块，用于根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型；

测量模块，用于根据所述流场模型数值模拟测量所述肺泡在收缩和扩展过程时的流场特征。

本发明肺泡收缩和扩展过程的流场数值模拟测量方法和系统，通过对肺泡建立几何模型，再对肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型，用以模拟肺泡在呼吸时的运动；接着获取所述肺泡对应的几何参数、呼吸参数及流体参数，根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，可模拟出所述肺泡收缩和扩展时肺泡内的流场模型，根据生成的流场模型可测量出所述肺泡在收缩和扩展过程时的流场特征，能缩短测量周期，减少运算负担，准确地获取肺泡内的流场特征数据。

附图说明

图1为本发明肺泡收缩和扩展过程的流场数值模拟测量方法在一实施例中的流程示意图。

图2为本发明肺泡收缩和扩展过程的流场数值模拟测量方法在一实施例中肺泡的几何模型示意图。

图3为本发明肺泡收缩和扩展过程的流场数值模拟测量方法在一实施例中肺泡的动网格模型示意图。

图4为本发明肺泡收缩和扩展过程的流场数值模拟测量方法在一实施例中肺泡及气管的模型示意图。

图5为一个呼吸周期内不同时刻时肺泡内部的流场流线图。

图6为吸气峰值时刻各级肺泡内部的流场流线图。

图7为呼气峰值时刻不同级肺泡内部的流线示意图。

图8为不同的呼吸强度下第七级肺泡在吸气阶段不同时刻肺泡内部的流场示意图。

图9为本发明肺泡收缩和扩展过程的流场数值模拟测量系统在一实施例中的结构示意图。

具体实施方式

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

如图1所示，是本发明肺泡收缩和扩展过程的流场数值模拟测量方法在一实施例中的流程示意图，包括如下步骤：

S11、根据预设的肺泡几何参数建立肺泡的几何模型，对所述肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型；

S12、获取所述肺泡对应的呼吸参数及流体参数；

S13、根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型；

S14、根据所述流场模型测量所述肺泡在收缩和扩展过程时的流场特征；

本实施例中，通过对肺泡建立几何模型，再对肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型，用以模拟肺泡在呼吸时的运动；接着获取所述肺泡对应的几何参数、呼吸参数及流体参数，根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，可模拟出所述肺泡收缩和扩展时肺泡内的流场模型，根据生成的流场模型可测量出所述肺泡在收缩和扩展过程时的流场特征，能缩短测量周期，减少运算负担，准确地获取肺泡内的流场特征数据。

本实施例中，考虑到自然界的流动可以分为层流流动(laminar)和湍流流动(turbulence)。从实验的角度来看，在层流流动中流体的层与层之间没有相互干扰，层之间没有质量的传递也没有动量的传递；而在湍流流动过程中，层与层之间的流体质点具有不断互相混渗的现象，速度和压力等物理量在时间和空间上具有随机性质的脉动值。

湍流是自然界普遍存在的流动，而层流属于个别的情况。常见的层流有毛细管或多孔介质中的流动、扰流物体表面边界层中的流动，层流一般出现在在流体的密度、特征速度和物体的特征长度较小，或者是流体的粘度很大等情况。同层流流动比较，湍流流动的研究比较复杂，至今一些基本的问题尚未完全解决。

根据人体呼吸道的几何形态，可以将空气在其中的流动视为管内流动，对于圆管内的流动可以定义Reynold数(Re数)为：

$\mathrm{Re}=\frac{\mathrm{ud}}{v}$

其中，u为流体的流速，d为管径，v为流体的动力学粘度。

一般认为Re≤2300时，管内的流动为层流流动；当Re≥8000时，管内的流动为湍流流动；当2300<Re<8000时，流动处于层流到湍流的过渡区。按照正常人的呼吸状态计算，气体在上呼吸道流动时，由于咽部和喉部形状起到喷管的作用，当地Re数可能超过8000，流动由层流转变为湍流状态。而在本实施例中人体呼吸道肺泡区，由于呼吸道的不断分叉，气体流动的横截面不断增大，在呼吸道末端处流速一般为cm/s，气管的直径为毫米左右，当地Re数远远小于100，属于典型的层流流动。

根据流动过程中流体的密度是否为常数可将流动区分为可压缩流动和不可压流动流动。当密度为常数时，流体为不可压缩流体；密度为变量时，流体为可压缩流动。虽然空气为可压缩流体，在本实施例的层流低速流动中，可以作为不可压缩流体来处理。

稳态流动和非稳态流动是依据流动过程速中速度、压力、温度等物理量是否随时间变化来区分。当各个物理量在流动过程中不随时间而变化时，流动为稳态定常流动；物理量在流动过程中随时间发生变化时，则为非稳态流动。

由于人体呼吸道的肺泡区的流动特征尺度为小于微米的量级，因此需考虑微小的尺度对气体流动分析的影响。在宏观尺度下研究流动问题的一个重要的假设就是连续介质假设，即在所研究的问题的特征尺度远远大于流体分子的平均自由程时，可以将流体视为一种连续体对待，认为流体的分子结构和分子热运动只能通过影响流体本身的热力学特性来间接影响流体的宏观运动。但是随着流动问题的几何尺度不断缩小，达到微米甚至是纳米量级时，流动过程就不能继续用传统的宏观理论来进行描述。

在研究小尺寸通道与宏观尺寸通道中流动现象之间的不同时，首先需要确定的是通道尺寸的划分问题，一般将大于1mm的尺度称为宏观尺度，1um～1mm之间的尺度称为微尺度。具体的尺度划分也没有统一的标准，有的是按照几何尺寸来进行划分，有的是根据流动中不同作用力的范围来进行尺度划分。

在微尺度流动中，一些在宏观流动中常常忽略不计的因素可能不再是可以忽略的。较小的特征尺度对于流动的影响主要体现于以下方面：

在宏观条件下，流体的粘度只与流体本身的性质有关。而在微观的条件下，流体的粘性受管道尺寸、横截面形状、温度、压强等多方面因素的影响，而目前还不能使用量化的方式准确地表达粘度同这些因素之间的关系。在微尺度流动中，就可能不能简单地将流体粘度视为常量来计算。

随着几何尺度的减小，微尺度流动中表面力将起主要的作用。流体虽然整体上不会呈现极性，但是，含有极性离子的流体与非极性流体的在流动特征上存在显著的差异。随着流动尺度的减小，粗糙度相对于管道的直径的比例增加，将导致表面粗糙度的影响显著增加。

Knudson数(Knudsonnumber)定义为:

$\mathrm{Kn}=\frac{\mathrm{\&lambda;}}{L}$

其中，λ为分子的平均自由程，L是所研究问题的特征尺度。

分子的平均自由程是指分子两次连续碰撞之间，一个分子自由运动的平均路程。平均自由程的计算式为：

$\mathrm{\&lambda;}=\frac{\mathrm{kT}}{\sqrt{2}\mathrm{\&pi;}{d}^{2}P}$

根据实际流动过程中Kn数的大小，可以将气体在管道内的流动分区如下：

(1)Kn≤0.001时，流动过程遵守连续介质假设，气体分子间的碰撞频率远高于气体分子和物体表面之间的碰撞频率，在这种情况下，牛顿粘性应力定律和傅立叶导热定律能够适用，可以求解具有无滑移边界条件的控制方程来描述流动过程；

(3)0.001≤Kn≤0.1时，流动处于滑移区，气体分子之间的碰撞仍远高于气体分子和物体表面的碰撞，但是气体和固体表面的交界处存在速度滑移和温度跳跃，可以求解具有滑移边界条件的方程来描述流动过程；

(4)0.1≤Kn≤10时，流动处于过渡区，这种情况下，分子之间的碰撞和分子和物体表面之间的碰撞处于同一个量级，需要同时考虑位于这个区域的流动对计算模型的要求很高，可以采用包含积分项的Boltzmann方程或对流动进行直接模拟的DMSC(directsimulationMonteCarlo)、分子动力学方法(MD)；

(5)Kn＞10时，流动为分子自由流。同气体分子和物体表面之间的碰撞相比较，分子间的碰撞可以忽略不计，可以采用无碰撞积分项的Boltzmann方程来描述流动过程。

空气分子的有效直径约为d＝3.1×10^-10m，可以根据kn数的定义计算出在温度为310K、压强为101325Pa时，即本实施例中的空气流动条件下，分子平均自由程大约为8.7×10^-8m。

按照人体的肺泡区主流气管的平均直径为0.5mm计算，可以的到流场的Kn数大约为Kn≈1.74×10^-4。因此根据Kn数的划分，本实施例中肺泡内流场需建立设定有无滑移边界条件的控制方程来进行测量。

下呼吸道的气体交换区几何形状十分复杂，本实施例建立的肺泡的几何模型如图2所示，图3是将几何模型进行网格划分后生成的动网格模型；

模型的建立可根据肺泡的几何参数而建立，几何参数包括肺泡气管管长、内腔直径、肺泡直径、两个肺泡之间的中心距离和肺泡开口角；本实施例的几何参数可为：管长(ductlength)740μm，内腔直径(lumenDiameter)320μm，肺泡直径(alveolardiameter)320μm，两肺泡间的中心距离为340μm，肺泡开口角(openingangle)120°，Re＝ud/v＝0.283。

模型中肺泡入口为速度入口，出口为压力出口边界；气流的入口速度按成人正常呼吸时的平均肺泡流量给定，成人正常坐立状态下的呼吸参数空气密度、呼吸周期和压力系数等参数；

其中，稳态时，Re＝0.283，空气的参数为：ρ＝1.225kg/m³，Re＝UD/ν，U＝0.0129m/s，μ＝1.7894e-05kg/m.s，ν＝μ/ρ＝1.46e-05m²/s。非稳态时，肺泡内的非稳态气流流量随时间成正弦曲线分布，呼吸周期T为4秒，计算为时间步长0.02秒，每个时间间隔迭代次数以残差值趋向稳定为限， $\stackrel{\&OverBar;}{u}=0.129m/s,$ 求得 $u_{\max }=\stackrel{\&OverBar;}{u}\&times;\mathrm{\&pi;}/2.$

在一较佳实施例中，所述根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型的步骤包括：

根据下列公式在所述肺泡的动网格模型中模拟所述肺泡内的流场：

$\left(\begin{array}{c}\frac{\&PartialD;u_i}{\&PartialD;x_i}=0\\ \frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i{u}_j\right)}{\&PartialD;x_j}=-\frac{\&PartialD;P}{\&PartialD;x_i}+\frac{\&PartialD;\left({\mathrm{\&tau;}}_{\mathrm{ij}}\right)}{\&PartialD;x_j}\\ \frac{\&PartialD;\left(\mathrm{\&rho;T}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_{j}T\right)}{\&PartialD;x_j}=\frac{\&PartialD;\left(\frac{k\&PartialD;T}{c_P\&PartialD;x_j}\right)}{\&PartialD;x_j}+S_T\end{array}\right.$

其中， ${\mathrm{\&tau;}}_{\mathrm{ij}}=\mathrm{\&mu;}(\frac{\&PartialD;u_i}{\&PartialD;x_j}+\frac{\&PartialD;u_j}{\&PartialD;x_i})-\frac{2\mathrm{\&mu;}{\mathrm{\&delta;}}_{\mathrm{ij}}\&PartialD;u_i}{3\&PartialD;x_j};$

i，j＝1，2，3；u_i和u_j为流体速度，x_i为流场变量，t为时间，ρ为流体密度，P为流体微元上的压力，τ_ij为粘性应力，c_P为比热容，T为温度，k为流体的传热系数，S_T为预设的粘性耗散项，μ为动力粘度。

所述根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型的步骤还可包括：

根据下列公式在所述肺泡的动网格模型中模拟所述肺泡的收缩和扩展过程：

$\begin{array}{c}\frac{{\left(\mathrm{\&rho;\&phi;V}\right)}^{n+1}-{\left(\mathrm{\&rho;\&phi;V}\right)}^n}{\mathrm{\&Delta;t}}+{\&Integral;}_{\&PartialD;V}\mathrm{\&rho;\&phi;}(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&PartialD;V}\mathrm{\&Gamma;}\&dtri;\mathrm{\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}\mathrm{dV}\end{array}$

$\begin{array}{c}\frac{\&PartialD;}{\&PartialD;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}\mathrm{dV}=\underset{\&PartialD;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;\mathrm{ndS}\end{array}$

$\begin{array}{c}\frac{\mathrm{dV}}{\mathrm{dt}}={\&Integral;}_{\&PartialD;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;\mathrm{dA}=\underset{l}{\overset{n_l}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l\end{array}$

其中， $V^{n+1}=V^n+\frac{\mathrm{dV}}{\mathrm{dt}}\mathrm{\&Delta;t},{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l=\frac{\mathrm{\&delta;}{V}_l}{\mathrm{\&Delta;t}};$

ρ为流体密度，为流体速度，为所述动网格模型中运动网格节点的速度，Γ为扩散系数，为控制体的边界，为边界的面积矢量；V为运动网格的单位体积，Vⁿ为n时刻所述运动网格的体积，dV/dt为控制体积对时间的导数，n为控制体积面的个数，S_φ为预设的广义源项，A_l为第l个面的面积矢量，δV_l表示第l个面在时间Δt内扫过的体积，D为肺泡所在管道的直径，S为控制面积。

所述根据所述流场模型测量所述肺泡在收缩和扩展过程时的流场特征的步骤可包括：

根据下式在所述流场模型中测量进出所述肺泡时的气流流量及所述肺泡在扩张和收缩时的气体流量：

$Q_a\left(t\right)=\frac{27}{16}\mathrm{\&pi;\&beta;f}{R}_a^2{\mathrm{\&lambda;}}^2\cos (\mathrm{ft}-\frac{\mathrm{\&pi;}}{2})$ 及 $Q_d\left(t\right)=\frac{d{V}_t}{\mathrm{dt}};$

其中，β＝(C+1)^1/3-1，C＝(V_max-V_min)/V_min；

Q_a为进出所述肺泡时的气流流量，Q_d为所述肺泡在扩张和收缩时的气体流量，β为肺泡的扩张幅度、V_max和V_min为所述几何模型的最大体积和最小体积，f为呼吸频率、t为呼吸时刻、R_a为肺泡的半径、λ为分子的平均自由程，V_t为肺泡总体积。

任何流动问题都必须满足质量守恒定律，可以得到流体的质量守恒方程(massconservationequation):

$\frac{\&PartialD;\mathrm{\&rho;}}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i\right)}{\&PartialD;x_i}=0(i,j=\mathrm{1,2,3})$

本实施例中流体可以视为不可压缩流体，将空气的密度视为常数，因此可以将上式简化为

$\frac{\&PartialD;u_i}{\&PartialD;x_i}=0$

动量守恒定律和能量守恒定律也是任何流动系统都必须满足的基本定律，根据动量守恒定律得到能量守恒方程(momentumconservationequation)：

$\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i{u}_j\right)}{\&PartialD;x_j}=-\frac{\&PartialD;P}{\&PartialD;x_i}+\frac{\&PartialD;\left({\mathrm{\&tau;}}_{\mathrm{ij}}\right)}{\&PartialD;x_j}$

上式中P是流体微元上的压力，τ_ij为粘性应力，可以表示为：

${\mathrm{\&tau;}}_{\mathrm{ij}}=\mathrm{\&mu;}(\frac{\&PartialD;u_i}{\&PartialD;x_j}+\frac{\&PartialD;u_j}{\&PartialD;x_i})-\frac{2\mathrm{\&mu;}{\mathrm{\&delta;}}_{\mathrm{ij}}\&PartialD;u_i}{3\&PartialD;x_j}$

根据能量守恒定律可得到能量守恒方程(energyconservationequation)：

$\frac{\&PartialD;\left(\mathrm{\&rho;T}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_{j}T\right)}{\&PartialD;x_j}=-\frac{\&PartialD;\left(\frac{k\&PartialD;T}{c_P\&PartialD;x_j}\right)}{\&PartialD;x_j}+S_T$

其中c_P为比热容，T为温度，k为流体的传热系数，S_T为粘性耗散项。对于不可压流动，如果热交换量很小时可以不考虑能量方程，而在本实施例中为了计算布朗力的作用，也需要求解能量方程。

上述的方程是流动所必须满足的瞬时控制方程，可为：

$\frac{\&PartialD;\left(\mathrm{\&rho;\&phi;}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_j\mathrm{\&phi;}\right)}{\&PartialD;x_j}=\frac{\&PartialD;\left(\frac{{\mathrm{\&Gamma;}}_{\mathrm{\&phi;}}\&PartialD;\mathrm{\&phi;}}{\&PartialD;x_j}\right)}{\&PartialD;x_j}+S_{\mathrm{\&phi;}}$

其中Γ_φ为广义扩散系数，S_φ为广义源项，上式中的各项依次为时间变化项、对流项、扩散项和源项。

对于描述流体运动的偏微分方程，只在极少数的情况下才能够得到分析解。而对于那些无法得到精确分析解的流动控制方程，只能使用数值离散的方法，将偏微分方程在时间和空间上离散，转化为代数方程进行计算机求解，数值方法的精度则取决于离散方法的选择。

数值模拟求解方法可包括分离求解法和耦合求解法。在分离求解法中，连续方程、动量方程和能量方程分别进行求解；在耦合求解法中，连续方程、动量方程和能量方程同时进行求解，在上述控制方程被求解后，再求解湍流、辐射等方程。这两种求解算法的求解对象是相同的，即它们所求解的控制方程均为描述质量守恒、动量守恒和能量守恒的连续方程、动量方程和能量方程，并且它们都用有限体积法作为对计算对象进行离散求解的基础方法。

有限体积法(FiniteVolumeMethod)又称为控制体积法。其主要包括如下步骤：

(1)将计算区域划分为一系列不重复的离散控制体积，并使每个网格点周围有一个控制体积；

(2)将待求解的微分控制方程对每一个控制体积积分，得出一组离散方程。

(3)将离散方程组线化，然后通过求解线性化方程组来获得变量的迭代解。

本实施例中的数值模拟即为分离求解法。在求解低速不可压缩问题时，速度、压力、温度之间的耦合程度不是很高，分离求解法具有很高的数值稳定性。分离求解法的计算步骤可以概括为：

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

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

(3)因为(2)中得到的速度场的数值解无法完全满足连续方程，于是再求解压强修正方程。压强修正方程是由连续方程导出的泊松型方程，求解这个方程可以得到对压强场、速度场和质量通量的修正，进而使连续方程得到满足。

(4)利用前面求出的解，求解湍流方程、能量方程、组元方程和能量方程。

(5)在多相流计算中如果考虑相间干扰，则需要通过求解弥散相轨迹计算得到连续相方程中的源项解。

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

在分离求解算法中，由于连续性方程和动量方程不需要耦合，由动量方程计算出的流场速度值可能不能同时满足连续性方程，解决的办法之一是用一个由连续性方程和线性化后的动量方程推导得到的“泊松类”方程来对压力进行修正，以得到满足连续性方程的压力和速度。分离求解算法中采用的压强速度耦合方法，包括SIMPLE(Semi-ImplicitMethodforPressure-LinkedEquations)方法、SIMPLEC方法和PISO(Pressure-ImplicitwithSplittingofOperators)方法三种。

SIMPLE的基本思想是：先使用假定的压强场来求解动量方程得到边界点上的通量，如果在计算中假定的压强场不准确，那么所求得的通量就必然不能满足连续方程。于是在通量上面添加修正项，使经过修正后所得的通量能够满足连续性方程。而添加得通量修正项是压强修正项的函数，因此将修正过的通量带入连续性方程，就可以得到一个关于压强修正项的方程，通过求解这个方程得到压强修正项的解。在压强修正项的前面乘以亚松弛因子，再与原压强相加就得到一个新的压强场。以这个新的压强场为起点重复上述过程，就形成交替求解压强场、速度场的迭代过程，直到最后得到收敛解。

PISO算法的主要思想是将SIMPLE和SIMPLEC算法中在求解压力修正方程时需要的重复计算移走。经过一个或者多个额外增加的PISO迭代步后，修正后的速度更加接近地满足连续性方程和动量方程。在每个迭代步中PISO算法增加少量的CPU时间，但是可以很大程度上减少迭代次数。PISO算法非常适用于瞬态问题的计算。

SIMPLEC算法是SIMPLE算法的改进算法之一，其算法与SIMPLE算法的基本思路一致，仅在通量修正方法上有所改进，加快计算的收敛速度。SIMPLEC格式的稳定性较好，在计算过程中可以将亚松弛因子适当的放大。在层流计算中如果没有使用辐射模型等辅助方程，使用SIMPLEC算法可以加快计算收敛速度。由于本实施例研究的流动是典型的层流流动，因此使用SIMPLEC算法来计算空气流场，以加快求解速度。

拉格朗日方法和欧拉方法是描述流体运动的两种基本方法：在拉格朗日描述下，网格随流体一起运动，网格点的速度和当地流体微团的速度相同；在欧拉描述下，网格的空间位置固定，即网格点的速度为零，流体微团穿过网格单元构成的控制体。任意拉格朗日欧拉方法(arbitarylagrangian-eulerian，ALE)将这两种基本方法统一起来，允许网格以任意的速度运动。当网格速度为零时则为欧拉方法，当网格运动的速度等于流体的速度时则为拉格朗日方法。

对于任意的标量φ，可以得到ALE形式有限体积控制方程：

$\begin{array}{c}\frac{d}{\mathrm{dt}}{\&Integral;}_V\mathrm{\&rho;\&phi;dV}+{\&Integral;}_{\&PartialD;V}\mathrm{\&rho;\&phi;}(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&PartialD;V}\mathrm{\&Gamma;}\&dtri;\mathrm{\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}\mathrm{dV}\end{array}$

其中ρ为流体的密度，为流体的速度，运动网格节点的速度，Γ为扩散系数，表示控制体的边界，表示边界的面积矢量。上式左端第一项时间导数项用一阶后差格式离散：

$\frac{d}{\mathrm{dt}}{\&Integral;}_V\mathrm{\&rho;\&phi;dV}=\frac{{\left(\mathrm{\&rho;\&phi;V}\right)}^{n+1}-{\left(\mathrm{\&rho;\&phi;V}\right)}^n}{\mathrm{\&Delta;t}}$

网格单元的体积不再是常量而是随时间变化的变量，n+1时刻的体积Vⁿ⁺¹可以写成：

$V^{n+1}=V^n+\frac{\mathrm{dV}}{\mathrm{dt}}\mathrm{\&Delta;t}$

式中，dV/dt为控制体积对时间的导数。

在应用运动网格数值方法时，为了避免网格的运动引入额外的误差，几何守恒律(GCL)和质量守恒、动量守恒、能量守恒一样需要在数值上得到满足。根据几何守恒定律：

$\begin{array}{c}\frac{\&PartialD;}{\&PartialD;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}\mathrm{dV}=\underset{\&PartialD;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;\mathrm{ndS}\end{array}$

$\begin{array}{c}\frac{\mathrm{dV}}{\mathrm{dt}}={\&Integral;}_{\&PartialD;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;\mathrm{dA}=\underset{j}{\overset{n_j}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,j}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_j\end{array}$

其中，n_j为控制体积面的个数，A_j是第j个面的面积矢量。其中

${\stackrel{\&RightArrow;}{u}}_{g,j}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_j=\frac{\mathrm{\&delta;}{V}_j}{\mathrm{\&Delta;t}}$

其中，δV_j表示第j个面在时间Δt内扫过的体积。在网格更新后，根据网格体积的变化，先求解出网格的速度，再求解流体的速度和其它的变量，因此可在所述肺泡的动网格模型中模拟所述肺泡内的流场。

本实施例建立肺泡的动网格模型进行模拟，在计算流体力学中将网格分为结构和非结构网格两大类。

结构网格是过去普遍采用的网格划分方法。在结构网格中，节点排列有序，相邻节点之间的关系明确。结构网格常用的生成方法上有代数生成方法、保角变换方法、偏微分方程方法、变分原理方法等。结构化网格的优点是可以准确地处理计算中的边界条件，计算精度高，并且可以使用很多高效的隐式算法和多重网格法，计算效率高。但是对于具有复杂外形的计算区域，结构化网格生成通常难以实现，而即便是生成多块结构化网格，块与块之间的界面条件处理也十分复杂，因此结构化网格的使用受到很大的限制。

在非结构网格中，网格单元和节点彼此之间没有固定的规律可循，节点的分布完全是任意的。其基本思想基于这样的假设：任何空间区域可以被四面体或者是三角形单元所填满。对于具有复杂形状的计算区域，也可以很方便的进行处理。采用非结构网格的划分方式也便于采用自适应的方法细化和粗化网格，使网格的划分更加方便和灵活，减少了在计算过程中划分网格的时间，对于参数变化剧烈、梯度很大的流动，也便于捕捉流场中的物理特性。同结构网格相比，非结构网格也存在计算量较大、高精度解算器发展还不是很成熟等缺点。

目前，在处理包含运动边界的流动问题时，一般基于非结构网格来处理网格的变形。常用的动网格变形方法有：

弹簧近似方法(SpringAnalog)：弹簧近似方法最早用于三角形非结构网格来求解翼型强迫震动扰流。弹性光顺法将系统中任意两个网格节点之间的边理想化为一个弹簧，系统的初始网格被看成是一个保持平衡的弹簧网格系统。系统中任意一个网格节点的位移都会导致与之相连接的弹簧中产生弹性力，进而导致邻近网格点上的力平衡被打破。由此波及出去，当经过反复的迭代至整个网格系统达到新的平衡时，得到一个变形后的网格系统。标准的弹簧近似方法由于没有考虑剪切应力的作用，变形能力不是很强，在二维和三维网格系统中，不能避免网格边的互相交叉。当存在较大的运动位移时，网格单元严重扭曲，容易发生控制体为非负的情况。针对上述的缺陷，可使用改进的弹簧近似方法，使网格的变形能力大大的提高并且保持了较高的效率。

弹簧近似方法又可以分为棱边弹簧(segmentspring)和顶点弹簧(vertexspring)两种描述形式。在棱边弹簧描述中，平衡长度等于边的原长，网格点在初始状态所受到的合力为零。根据胡克定律，当弹簧系统中网格点移动后，作用在其中某个网格节点i上的合力为：

${\stackrel{\&RightArrow;}{F}}_i=\underset{j}{\overset{n_i}{\mathrm{\&Sigma;}}}{k}_{\mathrm{ij}}(\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{x}}_j-\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{x}}_i)$

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

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

$\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{x}}_i^{m+1}=\frac{{\mathrm{\&Sigma;}}_i^{n_i}{k}_{\mathrm{ij}}\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{x}}_j^m}{{\mathrm{\&Sigma;}}_j^{n_i}{k}_{\mathrm{ij}}}$

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

${\stackrel{\&RightArrow;}{x}}_i^{n+1}={\stackrel{\&RightArrow;}{x}}_i^n+\mathrm{\&Delta;}{\stackrel{\&RightArrow;}{x}}_i^{\mathrm{converged}}$

和棱边弹簧描述方法不同的是，在顶点弹簧描述中弹簧的平衡长度为零。事实上，这两种弹簧近似方法的描述是一致的，其具体的证明过程可以参考其他的文献。为了使网格在变形以后仍能保持合理的疏密分布，弹簧的倔强系数可以有多种选择，其中应用最广泛的方法是假设弹簧的倔强系数和两节点之间边长的某次幂成反比。

弹性体方法(ElasticityMethod)：在弹性体方法中，整个网格区域被看成一个弹性体，从而所有网格点的移动由弹性力学的一组基本方程控制，通过求解这组方程，就能得到网格系统内部点的位移。弹性体方法在任何网格类型和运动类型中都能适用，并且具有较强的变形能力，更适合于由于物体旋转而引起网格旋转的情况，但是弹性体方法的计算工作量大，计算效率不如弹簧近似方法。

局部网格重构方法(LocalRemeshing)：局部网格重构方法最先在二维外挂物的投放数值模拟中发展。在大位移的问题中，当网格发生较大的变形时，可以对局部网格的进行重新生成，然后通过插值的方法求得新生成网格的流场属性值。局部网格重构方法一般包括调用自动网格生成器和由旧网格向新网格进行插值两部分。局部网格重构方法有以下的缺点：一是局部网格的重构方法需要在每个时间步长中重新生成网格，大大地增加了计算的负担；二是网格的重构过程还给数值计算带来了插值误差。在三维非结构动网格问题中，局部网格重构是一个非常繁重的工作。

可采用三种动网格的控制方法：弹性光滑法(Spring-basedSmoothing)、局部网格重划法(LocalRemeshing)和动态层法(DynamicLayering)。

弹性光滑法即所述的弹性近似方法。在Fluent中使用的是棱边弹簧描述方法，并且棱边的弹性系数取值为：

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

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

动态层方法是根据边界的移动量动态地增加或减少边界上网格层的技术，因此动态层方法适用于六面体网格、楔形网格等可以在边界上分层的网格系统。动态层技术在边界上假定一个优化的网格层高度值，在边界发生移动、变形时，如果临近边界的一层网格的高度同优化高度值相比大到一定程度时，就在边界面与相邻网格层之间增加一层网格。反之，如果临近边界的网格被压缩到一定程度时，临近一层网格又会被删除，用这种办法使边界上的网格保持在一定的密度。

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

在欧拉-欧拉方法中，不同的相被视为可以相互渗透的连续相，在其中引入相体积分数的概念，各相的体积分数之和为1，各相的体积分数是一个在时间和空间上连续的函数，根据守恒方程对每个相能得到具有相似结构的方程，通过引入另外经验的信息可以使方程组封闭求解。

在本实施例中，颗粒只占据两相流动中很小的体积分数，采用欧拉-拉格朗日方法来计算颗粒相的运动。

颗粒在气固两相流动过程中的受力是相当复杂，和颗粒与流体、颗粒与颗粒、颗粒与壁面之间的相互作用相关。一般而言，根据颗粒与流体之间的作用机理，两相间的作用力包括气体阻力、Basset力、Saffman升力、Magnus升力、虚假质量力、压力梯度力、热泳力、布朗力。其中，Basset力是颗粒在流体中加速或减速运动而产生的力，Saffman力是由于流场中存在速度梯度而引起的力，Magnus力是颗粒在流场中旋转产生的力，布朗力是微小颗粒在流场中受到流体分子的无规律碰撞作用产生的力。

虽然颗粒在两相流动中的受力十分复杂，但是并不是所有的力都起同等重要的作用。对于本文研究的下呼吸道肺泡区，气体的流速很小，在本实施例的计算中忽略掉颗粒所受的Basset力、Saffman升力、Magnus升力、虚假质量力。对于能够进入呼吸道肺泡区的颗粒，其温度已经接近气管或肺泡壁面的温度，也可将热泳力、电泳力等力的作用予以忽略不计。层流流动中的微小颗粒，其受到的布朗力同重力和颗粒受到的阻力具有相同的量级，是不可忽略的。

表2.1中列出了标准的大气压下、温度为26度的静止空气中，不同粒径的颗粒在布朗扩散作用和重力沉降作用下1s内的位移量，可以看出在对于粒径大于1微米的颗粒，布朗力的作用相对于重力来说可以忽略不计，而对于小于0.5um的颗粒，布朗力的作用要大于重力。

表2.1布朗力和重力的对颗粒运动的相对影响

在拉格朗日框架下，颗粒的运动方程通过计算颗粒上的受力平衡，可通过下式计算：

$\frac{d{u}_p}{\mathrm{dt}}=F_D(u-u_p)+\frac{g_i({\mathrm{\&rho;}}_p-\mathrm{\&rho;})}{{\mathrm{\&rho;}}_p}+F_B$

其中，u为空气流场的速度，u_p为颗粒的速度，ρ_p为颗粒的密度，ρ为流体的密度，F_D为作用在颗粒上的拽力，g_i为重力矢量，F_B是作用于颗粒上的布朗力。

对于本实施例所研究的微小颗粒而言，分子间的滑移作用可能会起很大的作用。因此，本实施例通过下式计算其中的阻力FD：

$F_D=\frac{18\mathrm{\&mu;}}{d_p^2{\mathrm{\&rho;}}_p{C}_C}$

$C_C=1+\frac{2\mathrm{\&lambda;}}{d_p}[1.257+0.4\exp (-1.1d_p/2\mathrm{\&lambda;})]$

μ为流体的粘度，C_C为Cunningham滑移修正系数，λ为分子的平均自由程，d_p为颗粒的直径。

颗粒运动方程中的布朗力F_B，可以看成是一个随机的矢量，用颗粒和流体的物性和计算求解时间步长来定义：

$S_0=\frac{216v{k}_{B}T}{{\mathrm{\&pi;}}^2\mathrm{\&rho;}{d}_p^5{\left(\frac{{\mathrm{\&rho;}}_p}{\mathrm{\&rho;}}\right)}^2{C}_C}$

上式中，为一个具有零平均值和单位方差的高斯随机矢量，S₀为白噪声过程的幅度，v为空气的动力粘度，k_B＝1.38×10^-16为Boltzmann常数，T为空气流场的温度。

通过对颗粒的运动方程积分得到颗粒的速度矢量u_p，颗粒的在计算区域的位置根据下式的积分求得：

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

对上式进行积分，就可以得到各个时刻颗粒在流场中的位置。

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

如图4所示，是本实施例肺泡及气管的模型示意图，本实施例中使用前后开口的圆柱管模拟呼吸道中的气管，使用连接在圆柱壁面上球形盖来描述呼吸道中的肺泡。本实施例中R_a和R_d分别为肺泡和气管的半径，γ为肺泡的开口角。

本实施例中气管的半径R_d可为250um，肺泡的半径为200um，肺泡在主气管上的开口角设为60度。气管的长度为1000um，并且在本实施例的计算模型中，肺泡位于该段气管的正中位置。

可以将人体肺泡和气管在呼吸过程中的变形过程视为一个保持几何形状相似的扩张和收缩运动。肺泡和气管的各个特征尺寸为随时间变化的正弦函数：

$L\left(t\right)=L_0[1+\frac{\mathrm{\&beta;}}{2}+\frac{\mathrm{\&beta;}}{2}\sin (\mathrm{ft}-\frac{\mathrm{\&pi;}}{2})]=L_0\mathrm{\&lambda;}$

其中L₀为几何模型中的各个特征尺寸在时刻t＝0时(吸气开始时或呼气结束时)的值，f＝2π/T是呼吸的频率，T为呼吸的周期时间，β为各个特征尺度的扩张幅度。

可以定义肺泡的特征体积扩张系数(excursionfactor)为：

C＝(V_max-V_min)/V_min

其中V_max和V_min为几何模型的最大和最小体积，分别对应在实际呼吸过程中呼气终了和吸气终了时的体积。由此可以得到肺泡的特征尺度扩张系数：

β＝(C+1)^1/3-1

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

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

Re(t)＝U_d(t)D_d(t)/υ

$U_d\left(t\right)=4Q_d\left(t\right)/\left(\mathrm{\&pi;}{D}_d^2\left(t\right)\right)$

可通过一个无量纲特征参数来表征肺泡在呼吸道上所处的位置，定义为呼吸过程中进出肺泡的气体流量Q_a与流过肺泡所处气管的气体流量Q_d的比值。

在肺泡区的每一级气管上，根据气体不可压缩的假设，进出肺泡的气体流量为：

$Q_a\left(t\right)=\frac{d{V}_a}{\mathrm{dt}}---(3-7)$

V_a为气管上单位肺泡的体积。根据几何模型可以得到：

$Q_a\left(t\right)=\frac{27}{16}\mathrm{\&pi;\&beta;f}{R}_a^2{\mathrm{\&lambda;}}^2\cos (\mathrm{ft}-\frac{\mathrm{\&pi;}}{2})---(3-8)$

流过该级气管的气体流量可以下面的表达式求得：

$Q_d\left(t\right)=\frac{d{V}_t}{\mathrm{dt}}---(3-9)$

V_t为经过该段气管后的所有气管及这些气管上所有肺泡的体积之和，Q_d即这部分空间扩张和收缩变化产生的流量。

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

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

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

在本实施例的模型中，肺泡的边界的变形量大于位于边界上网格的高度，需要使用弹性光顺法和局部重划法来控制模型中网格的变形运动，因此对于本实施例的计算模型的网格划使用三角形网格。由于三维模型的计算中，为了使网格的边界不发生局部的变形，就需要在计算中使用较小的时间步长，而且在每一个时间步长中，为了使整个弹簧系统达到平衡，需要进行多次的迭代，这给计算带来了很大的负担。因此，在本实施例中，只对简化的二维模型进行数值模拟。

根据模拟出的所述肺泡收缩和扩展时肺泡内的流场模型，可得到出肺泡内稳态气流的3D速度等值线图和流线图，从等值线图和流线图可知，肺腺泡内的速度和主气道内的速度相比是可以忽略不计的，肺腺泡后的主气道内速度出现最大值，在圆柱主气道内，速度沿着肺泡延伸有增大趋势，速度值随着半径逐渐减小；同时，XZ截面和XY界面上的速度分布一样，可见气流流动不受重力的影响。

根据模拟出的所述肺泡收缩和扩展时肺泡内的流场模型，也可得到肺泡内的不同界面位置上的流线图，3D肺泡内的流线特别是肺腺泡内的流线结构复杂。肺泡内同样有回流区的存在，但是回流区不仅有非平面结构，而且也有立体结构的回流，肺泡内从外到内，回流区逐渐从平面结构转变为立体结构，回流也变得越来越复杂。在圆柱主气道内，气流流线就变得简单的多，气流基本上和X轴平行，在腺泡入口处，逐渐产生向泡内的绕流。同时，主气流和回流之间不存在对流交换，绕入的气流最终还是从开口处绕出，说明回流区和主气流之间存在一个分隔面，绕入气流进入回流。

利用模拟出的所述肺泡收缩和扩展时肺泡内的流场模型，也可输出呼吸过程的各个特征时刻时的气流的流线和速度等值线图，在本实施例中，除呼吸转换时刻(1/2T和T时刻)外，其他时刻，气流速度在肺泡内的分布是相似的，只是速度的大小不相同。需要特别指出的是，横向截面上的速度分布是一样的，纵截面上，无论是呼气还是吸气，速度最大值都是在压力出口处，速度入口处的速度在呼入和呼出时速度都是均匀的。同时，主气道内的气流和腺泡内的气流相比，腺泡内的气流使可以忽略不计的。从截面上的速度流线可以看出，腺泡内存在回流区，回流去的位置发生变换，同时回流区和主气流之间同样存在分隔线(面)。在呼气和吸气转换的时刻(1/2T和T时刻)，气流发生了明显的改变，原有气流场在此时刻发生了变化，这一时刻回流混合主气流之间发生了明显的物质对流交换，此时气流的扰动也是最大的。

1/4T时刻时各个肺泡内的流线是相同的，肺泡内的回流要复杂的多，同时肺泡回流和主气流之间存在分隔界面。各个时刻的肺泡内的回流中心是在发生稍许的变化的，这表明非稳态气流能够产生部分的气流扰动，使得肺泡回流和主气流之间发生轻微的对流交换，从而产生了呼气和吸气过程中的回流中心的变化。但是这种扰动效应要远比呼吸转换过程中所发生的气流扰动要小的多。同时，由于回流区的复杂，因此不能够明显的判断主气流绕入肺泡的弧度是否随着呼吸过程而发生了明显的改变。

利用模拟出的所述肺泡收缩和扩展时肺泡内的流场模型，也可测量出周期性壁面收缩扩张运动下肺泡内的流场特征。流动过程中，模型内气管和呼吸道没有扩张的情形，尽管这种呼吸状态在实际呼吸过程中并不存在，但是为了比较扩张收缩壁面运动对肺泡内部流场的影响，可先计算无壁面运动的非稳态呼吸过程。在这一过程中，气体的进口速度仍为随时间变化的正弦波分布。

如图5所示，是一个呼吸周期内不同时刻时肺泡内部的流场示意图，图5a和图5b分别为1/10T周期时肺泡壁面静止和肺泡壁面扩张收缩运动两种情况下第一级气管上肺泡内部的流场示意图；图5c和图5d分别为1/10T周期时肺泡壁面静止和肺泡壁面扩张收缩运动两种情况下第一级气管上肺泡内部的流场示意图；图5e和图5f分别为1/10T周期时肺泡壁面静止和肺泡壁面扩张收缩运动两种情况下第一级气管上肺泡内部的流场示意图；图5g和图5h分别为1/10T周期时肺泡壁面静止和肺泡壁面扩张收缩运动两种情况下第一级气管上肺泡内部的流场示意图.

本实施例中肺泡扩张系数为0.05，可以看出，在没有肺泡的壁面的情况下，气管内的流动在肺泡开口处形成了一个剪切流，在肺泡区内形成了一个回旋的流动区域。整个流场可以分为两个区域，分别为气管内的流动和肺泡内的流动，能近似地使用一条连接肺泡开口处前后两角的连线分开。对于实际三维肺泡结构，这条分隔线可以看成是一个曲面。气管内的流体和肺泡内的流体几乎没有对流混合的作用。气管内的气流朝着出口方向作对流运动，肺泡内部的流体以平行于肺泡壁面作回旋运动，回旋运动的区域占据整个肺泡空腔，这是外略空腔剪切流的主要特征。肺泡内的回旋流动的速度很小，远远小于气管内的流动速度，只有主流速度百分之一的量级。这一流场形态和剪切流外略轴对称空腔的研究中得到的结果是符合的。理论上，对于轴对称的空腔，在剪切流的速度很小的时候，回旋区域也应该是轴对称的。第一级气管上由于主流气管中的气体速度还比较大，因此回旋区域的中心略偏向上游的一侧。

尽管肺泡的气体流量同气管的气体流量(Q_a/Q_d＝0.004％)相比较，几乎是一个可以忽略不计的值。但同没有肺泡壁面运动情况下肺泡内部的流场结构相比较，可以看到肺泡内部的流场特征有了明显的差别。气管中的剪切流在肺泡内部仍然形成了一个很大的回流区域。但是壁面的运动对流场产生了作用：气管内的一小部分的流体同肺泡内的流体发生了对流混合，这部分流体来自于气管靠近壁面的位置。肺泡的扩张给了肺泡内作回旋运动的流体一个切向的速度，给气管内的气体一个径向速度，使气管靠近壁面的流体和气管内的主流脱离开来，这一部分流体随对流运动进入肺泡后，在肺泡开口靠近出口的位置方向发生转折，并沿着肺泡的壁面运动。肺泡内部的回旋流动区域的形态也发生了很大的变化。没有肺泡壁面运动的流场中，可以看到回旋流动是轴对称的。肺泡内回旋流动的中心也不断地发生偏移，这个偏移的方向和吸气过程中的偏移方向相反，和气管内的气体流动方向相同，回旋区域的面积随着偏移减小。

如图6所示，是吸气峰值时刻各级肺泡内的流场示意图，从图中可以看出，肺泡的壁面处于扩张运动，肺泡内回旋流动的中心偏移到了靠近肺泡开口靠近气管上游的位置。在吸气过程中，回旋流动的中心朝着肺泡的轴线移动，移动和主气管内的气体流动方向相同，回旋区域的面积增加，在吸气结束时，回旋已经比较接近肺泡的轴线。

如图7所示，是呼气峰值时刻不同级上肺泡内部的流线示意图；呼气的过程中，肺泡处于收缩运动，这时肺泡内的流线和吸气过程同一时刻的方向正好相反。

尽管气管内部的流体和肺泡内的流体发生了对流混合，但只发生在气管内接近管壁的一小部分流体上，主流上的大部分流体与肺泡内的流体没有发生对流混合，这部分流体朝着气管出口方向对流运动，进入气管的下一级。随着肺泡在呼吸道上的位置接近终端，肺泡和气管内气流的对流作用就越明显。

当Q_a/Q_d增大到某一值时，在整个呼吸过程中，肺泡内都不会出现回旋区域，可以观察到，肺泡的第七级上面仍然出现了回旋区域，因此人体呼吸道上的绝大多数肺泡内都存在回旋运动区域。

如图8所示，是不同的呼吸强度下(左为10L/min，右为20L/min)，第七级肺泡在吸气阶段不同时刻肺泡内部的流场示意图；在不同的肺泡扩张系数下，同一肺泡的内的流动也略有不同。扩张系数由0.05增加到0.1，肺泡的扩张幅度增加一倍，此时气管内的流量正好也增加了一倍，气管内主流气体速度增加产生的作用更明显一些，从图8左右两列图的对比中可以看出，气管主流速度产生作用增加更明显。

如图9所示，是本发明一种肺泡收缩和扩展过程的流场数值模拟测量系统在一实施例中的结构示意图，包括：

建立模块91，用于根据预设的肺泡的几何参数建立肺泡的几何模型，对所述肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型；

获取模块92，用于获取所述肺泡对应的呼吸参数及流体参数；

模拟模块93，用于根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型；

测量模块94，用于根据所述流场模型数值模拟测量所述肺泡在收缩和扩展过程时的流场特征。

在一较佳实施例中，所述模拟模块93还用于：

根据下列公式在所述肺泡的动网格模型中模拟所述肺泡内的流场：

$\left(\begin{array}{c}\frac{\&PartialD;u_i}{\&PartialD;x_i}=0\\ \frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_i{u}_j\right)}{\&PartialD;x_j}=-\frac{\&PartialD;P}{\&PartialD;x_i}+\frac{\&PartialD;\left({\mathrm{\&tau;}}_{\mathrm{ij}}\right)}{\&PartialD;x_j}\\ \frac{\&PartialD;\left(\mathrm{\&rho;T}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;}{u}_{j}T\right)}{\&PartialD;x_j}=\frac{\&PartialD;\left(\frac{k\&PartialD;T}{c_P\&PartialD;x_j}\right)}{\&PartialD;x_j}+S_T\end{array}\right.$

其中， ${\mathrm{\&tau;}}_{\mathrm{ij}}=\mathrm{\&mu;}(\frac{\&PartialD;u_i}{\&PartialD;x_j}+\frac{\&PartialD;u_j}{\&PartialD;x_i})-\frac{2\mathrm{\&mu;}{\mathrm{\&delta;}}_{\mathrm{ij}}\&PartialD;u_i}{3\&PartialD;x_j};$

i，j＝1，2，3；u_i和u_j为流体速度，x_i为流场变量，t为时间，ρ为流体密度，P为流体微元上的压力，τ_ij为粘性应力，c_P为比热容，T为温度，k为流体的传热系数，S_T为预设的粘性耗散项，μ为动力粘度；

和/或

根据下列公式在所述肺泡的动网格模型中模拟所述肺泡的收缩和扩展过程：

$\begin{array}{c}\frac{{\left(\mathrm{\&rho;\&phi;V}\right)}^{n+1}-{\left(\mathrm{\&rho;\&phi;V}\right)}^n}{\mathrm{\&Delta;t}}+{\&Integral;}_{\&PartialD;V}\mathrm{\&rho;\&phi;}(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&PartialD;V}\mathrm{\&Gamma;}\&dtri;\mathrm{\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}\mathrm{dV}\end{array}$

$\begin{array}{c}\frac{\&PartialD;}{\&PartialD;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}\mathrm{dV}=\underset{\&PartialD;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;\mathrm{ndS}\end{array}$

$\begin{array}{c}\frac{\mathrm{dV}}{\mathrm{dt}}={\&Integral;}_{\&PartialD;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;\mathrm{dA}=\underset{l}{\overset{n_l}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l\end{array}$

其中， $V^{n+1}=V^n+\frac{\mathrm{dV}}{\mathrm{dt}}\mathrm{\&Delta;t},{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l=\frac{\mathrm{\&delta;}{V}_l}{\mathrm{\&Delta;t}};$

ρ为流体密度，为流体速度，为所述动网格模型中运动网格节点的速度，Γ为扩散系数，为控制体的边界，为边界的面积矢量；V为运动网格的单位体积，Vⁿ为n时刻所述运动网格的体积，dV/dt为控制体积对时间的导数，n为控制体积面的个数，S_φ为预设的广义源项，A_l为第l个面的面积矢量，δV_l表示第l个面在时间Δt内扫过的体积，D为肺泡所在管道的直径，S为控制面积。

在一较佳实施例中，所述测量模块94还用于：

根据下式在所述流场模型中测量进出所述肺泡时的气流流量及所述肺泡在扩张和收缩时的气体流量：

$Q_a\left(t\right)=\frac{27}{16}\mathrm{\&pi;\&beta;f}{R}_a^2{\mathrm{\&lambda;}}^2\cos (\mathrm{ft}-\frac{\mathrm{\&pi;}}{2})$ 及 $Q_d\left(t\right)=\frac{d{V}_t}{\mathrm{dt}};$

其中，β＝(C+1)^1/3-1，C＝(V_max-V_min)/V_min；

Q_a为进出所述肺泡时的气流流量，Q_d为所述肺泡在扩张和收缩时的气体流量，β为肺泡的扩张幅度、V_max和V_min为所述几何模型的最大体积和最小体积，f为呼吸频率、t为呼吸时刻、R_a为肺泡的半径、λ为分子的平均自由程，V_t为肺泡总体积。

本发明肺泡收缩和扩展过程的流场数值模拟测量方法和系统，通过对肺泡建立几何模型，再对肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型，用以模拟肺泡在呼吸时的运动；接着获取所述肺泡对应的几何参数、呼吸参数及流体参数，根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，可模拟出所述肺泡收缩和扩展时肺泡内的流场模型，根据生成的流场模型可测量出所述肺泡在收缩和扩展过程时的流场特征，能缩短测量周期，减少运算负担，准确地获取肺泡内的流场特征数据。

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

Claims (9)

1.一种肺泡收缩和扩展过程的流场数值模拟测量方法，其特征在于，包括如下步骤：

根据预设的肺泡的几何参数建立肺泡的几何模型，对所述肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型；

获取所述肺泡对应的呼吸参数及流体参数；

根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型；

根据所述流场模型数值模拟测量所述肺泡在收缩和扩展过程时的流场特征；

其中，所述根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型的步骤包括：

根据下列公式在所述肺泡的动网格模型中数值模拟所述肺泡内的流场：

$\left(\begin{array}{c}\frac{\&part;u_i}{\&part;x_i}=0\\ \frac{\&part;\left({\mathrm{\&rho;u}}_i\right)}{\&part;t}+\frac{\&part;\left({\mathrm{\&rho;u}}_i{u}_j\right)}{\&part;x_j}=-\frac{\&part;P}{\&part;x_i}+\frac{\&part;\left({\mathrm{\&tau;}}_{ij}\right)}{x_j}\\ \frac{\&part;\left(\mathrm{\&rho;}T\right)}{\&part;t}+\frac{\&part;\left({\mathrm{\&rho;u}}_{j}T\right)}{\&part;x_j}=\frac{\&part;\left(\frac{k\&part;T}{c_P\&part;x_j}\right)}{\&part;x_j}+S_T\end{array}\right.$

其中， ${\mathrm{\&tau;}}_{ij}=\mathrm{\&mu;}(\frac{\&part;u_i}{\&part;x_j}+\frac{\&part;u_j}{\&part;x_i})-\frac{2{\mathrm{\&mu;\&delta;}}_{ij}\&part;u_i}{3\&part;x_j};$

i，j＝1，2，3；u_i和u_j为流体速度，x_i为流场变量，t为时间，ρ为流体密度，P为流体微元上的压力，τ_ij为粘性应力，c_P为比热容，T为温度，k为流体的传热系数，S_T为预设的粘性耗散项，μ为动力粘度，x_j为位置坐标，为导数符号。

2.根据权利要求1所述的肺泡收缩和扩展过程的流场数值模拟测量方法，其特征在于，所述肺泡的几何参数包括肺泡气管管长、内腔直径、肺泡直径、两个肺泡之间的中心距离和肺泡开口角。

3.根据权利要求1所述的肺泡收缩和扩展过程的流场数值模拟测量方法，其特征在于，所述呼吸参数包括空气密度、呼吸周期和压力系数。

4.根据权利要求1所述的肺泡收缩和扩展过程的流场数值模拟测量方法，其特征在于，流体参数包括气流的雷诺数、流体密度、流体速度和流体粘度。

5.根据权利要求1所述的肺泡收缩和扩展过程的流场数值模拟测量方法，其特征在于，所述根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型的步骤还包括：

根据下列公式在所述肺泡的动网格模型中模拟所述肺泡的收缩和扩展过程：

$\left(\begin{array}{c}\frac{{\left(\mathrm{\&rho;}\mathrm{\&phi;}V\right)}^{n+1}-{\left(\mathrm{\&rho;}\mathrm{\&phi;}V\right)}^n}{\mathrm{\&Delta;}t}+{\&Integral;}_{\&part;V}\mathrm{\&rho;}\mathrm{\&phi;}\left(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g\right)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&part;V}\mathrm{\&Gamma;}\&dtri;\mathrm{\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}dV\\ \frac{\&part;}{\&part;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}dV=\underset{\&part;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;ndS\\ \frac{dV}{dt}={\&Integral;}_{\&part;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;dA=\underset{l}{\overset{n_l}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l\end{array}\right.$

其中， $V^{n+1}=V^n+\frac{dV}{dt}\mathrm{\&Delta;}t,$ ${\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l=\frac{{\mathrm{\&delta;V}}_l}{\mathrm{\&Delta;}t};$

ρ为流体密度，Δt为时间步长，φ为变量符号，▽为偏导数符号，为流体速度，为所述动网格模型中运动网格节点的速度，Γ为扩散系数，V为控制体的边界，为边界的面积矢量；Ω(t)为曲面的体积，dS为曲面的面积；为微元面积，l表示X方向,n_l表示法向量，表示矢量；V为运动网格的单位体积，Vⁿ为n时刻所述运动网格的体积，dV/dt为控制体积对时间的导数，n为控制体积面的个数，S_φ为预设的广义源项，A_l为第l个面的面积矢量，δV_l表示第l个面在时间Δt内扫过的体积，D为肺泡所在管道的直径，S为控制面积。

6.根据权利要求1所述的肺泡收缩和扩展过程的流场数值模拟测量方法，其特征在于，所述根据所述流场模型测量所述肺泡在收缩和扩展过程时的流场特征的步骤包括：

根据下式在所述流场模型中测量进出所述肺泡时的气流流量及所述肺泡在扩张和收缩时的气体流量：

$Q_a\left(t\right)=\frac{27}{16}{\mathrm{\&pi;\&beta;fR}}_a^2{\mathrm{\&lambda;}}^{2}cos(ft-\frac{\mathrm{\&pi;}}{2})$ 及 $Q_d\left(t\right)=\frac{{\mathrm{dV}}_t}{dt};$

其中，β＝(C+1)^1/3-1，C＝(V_max-V_min)/V_min；

Q_a为进出所述肺泡时的气流流量，Q_d为所述肺泡在扩张和收缩时的气体流量，β为肺泡的扩张幅度，C为常数，t表示时间，cos表示余弦函数，V_max和V_min为所述几何模型的最大体积和最小体积，f为呼吸频率、t为呼吸时刻、R_a为肺泡的半径、λ为分子的平均自由程，V_t为肺泡总体积。

7.一种肺泡收缩和扩展过程的流场数值模拟测量系统，其特征在于，包括：

建立模块，用于根据预设的肺泡的几何参数建立肺泡的几何模型，对所述肺泡的几何模型进行网格划分，生成所述肺泡的动网格模型；

获取模块，用于获取所述肺泡对应的呼吸参数及流体参数；

模拟模块，用于根据所述肺泡的动网格模型及几何参数、呼吸参数及流体参数，模拟所述肺泡收缩和扩展时肺泡内的流场模型；

测量模块，用于根据所述流场模型数值模拟测量所述肺泡在收缩和扩展过程时的流场特征；

其中，所述模拟模块还用于：

根据下列公式在所述肺泡的动网格模型中模拟所述肺泡内的流场：

$\left(\begin{array}{c}\frac{\&part;u_i}{\&part;x_i}=0\\ \frac{\&part;\left({\mathrm{\&rho;u}}_i\right)}{\&part;t}+\frac{\&part;\left({\mathrm{\&rho;u}}_i{u}_j\right)}{\&part;x_j}=-\frac{\&part;P}{\&part;x_i}+\frac{\&part;\left({\mathrm{\&tau;}}_{ij}\right)}{x_j}\\ \frac{\&part;\left(\mathrm{\&rho;}T\right)}{\&part;t}+\frac{\&part;\left({\mathrm{\&rho;u}}_{j}T\right)}{\&part;x_j}=\frac{\&part;\left(\frac{k\&part;T}{c_P\&part;x_j}\right)}{\&part;x_j}+S_T\end{array}\right.$

其中， ${\mathrm{\&tau;}}_{ij}=\mathrm{\&mu;}(\frac{\&part;u_i}{\&part;x_j}+\frac{\&part;u_j}{\&part;x_i})-\frac{2{\mathrm{\&mu;\&delta;}}_{ij}\&part;u_i}{3\&part;x_j};$

i，j＝1，2，3；u_i和u_j为流体速度，x_i为流场变量，t为时间，ρ为流体密度，P为流体微元上的压力，τ_ij为粘性应力，c_P为比热容，T为温度，k为流体的传热系数，S_T为预设的粘性耗散项，μ为动力粘度，x_j为位置坐标，为导数符号。

8.根据权利要求7所述的肺泡收缩和扩展过程的流场数值模拟测量系统，其特征在于，所述模拟模块还用于：

根据下列公式在所述肺泡的动网格模型中模拟所述肺泡的收缩和扩展过程：

$\left(\begin{array}{c}\frac{{\left(\mathrm{\&rho;}\mathrm{\&phi;}V\right)}^{n+1}-{\left(\mathrm{\&rho;}\mathrm{\&phi;}V\right)}^n}{\mathrm{\&Delta;}t}+{\&Integral;}_{\&part;V}\mathrm{\&rho;}\mathrm{\&phi;}\left(\stackrel{\&RightArrow;}{u}-{\stackrel{\&RightArrow;}{u}}_g\right)\&CenterDot;d\stackrel{\&RightArrow;}{A}={\&Integral;}_{\&part;V}\mathrm{\&Gamma;}\&dtri;\mathrm{\&phi;}\&CenterDot;d\stackrel{\&RightArrow;}{A}+{\&Integral;}_V{S}_{\mathrm{\&phi;}}dV\\ \frac{\&part;}{\&part;t}\underset{\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;\&Integral;}dV=\underset{\&part;\mathrm{\&Omega;}\left(t\right)}{\&Integral;\&Integral;}D\&CenterDot;ndS\\ \frac{dV}{dt}={\&Integral;}_{\&part;V}{\stackrel{\&RightArrow;}{u}}_g\&CenterDot;dA=\underset{l}{\overset{n_l}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l\end{array}\right.$

其中， $V^{n+1}=V^n+\frac{dV}{dt}\mathrm{\&Delta;}t,$ ${\stackrel{\&RightArrow;}{u}}_{g,l}\&CenterDot;{\stackrel{\&RightArrow;}{A}}_l=\frac{{\mathrm{\&delta;V}}_l}{\mathrm{\&Delta;}t};$

ρ为流体密度，Δt为时间步长，φ为变量符号，▽为偏导数符号，为流体速度，为所述动网格模型中运动网格节点的速度，Γ为扩散系数，V为控制体的边界，为边界的面积矢量；Ω(t)为曲面的体积，dS为曲面的面积；为微元面积，l表示X方向,n_l表示法向量，表示矢量；V为运动网格的单位体积，Vⁿ为n时刻所述运动网格的体积，dV/dt为控制体积对时间的导数，n为控制体积面的个数，S_φ为预设的广义源项，A_l为第l个面的面积矢量，δV_l表示第l个面在时间Δt内扫过的体积，D为肺泡所在管道的直径，S为控制面积。

9.根据权利要求7所述的肺泡收缩和扩展过程的流场数值模拟测量系统，其特征在于，所述测量模块还用于：

根据下式在所述流场模型中测量进出所述肺泡时的气流流量及所述肺泡在扩张和收缩时的气体流量：

$Q_a\left(t\right)=\frac{27}{16}{\mathrm{\&pi;\&beta;fR}}_a^2{\mathrm{\&lambda;}}^{2}cos(ft-\frac{\mathrm{\&pi;}}{2})$ 及 $Q_d\left(t\right)=\frac{{\mathrm{dV}}_t}{dt};$

其中，β＝(C+1)^1/3-1，C＝(V_max-V_min)/V_min；

Q_a为进出所述肺泡时的气流流量，Q_d为所述肺泡在扩张和收缩时的气体流量，β为肺泡的扩张幅度，C为常数，t表示时间，cos表示余弦函数，V_max和V_min为所述几何模型的最大体积和最小体积，f为呼吸频率、t为呼吸时刻、R_a为肺泡的半径、λ为分子的平均自由程，V_t为肺泡总体积。