真实人体呼吸道湍流仿真方法分析

崔 岩， 王玉龙， 董榕波， 李 蕾， 李 海， 马松云， 刘银水

(1.华中科技大学 机械科学与工程学院， 湖北 武汉 430074；2.中国人民解放军95829医院， 湖北 武汉 430022；3.华中科技大学 同济医学院附属协和医院儿科， 湖北 武汉 430030；4.亚琛工业大学 机械工程学院通用力学研究所， 德国 亚琛 52062)

引言

吸入给药是利用肺泡具有较大表面积、毛细血管丰富、物质交换短和药物吸收快等优点，将药物粉末/液滴分散，经由口腔进入呼吸道，最终进入血液循环以达到预防和治疗呼吸系统疾病的一种高效治疗手段[1]。目前吸入给药存在药物输送效率较低的突出问题，会急剧增加患者的医疗成本，若过量用药则会损害患者的身体健康。药物颗粒在呼吸道内的沉积主要受惯性沉积、湍流扩散沉积以及颗粒布朗运动的影响[2]，其中湍流扩散沉积起主导作用。因此，精准的求解人体呼吸道的湍动能场，是研究提升吸入给药效率的根本前提。

人体支气管树结构复杂，呼吸道截面变化巨大，其流动雷诺数存在较大的变化范围(Re=500～8000)[3]。因此，需借助可靠的湍流模型预测包含层流到湍流过渡转变、剪切流和回旋流等复杂流态。目前，何种湍流模型能够高效精准求解呼吸道流场依旧是一个开放性问题[4]。使用直接数值模拟(Direct Numerical Simulation, DNS)的方法求解瞬时Navier-Stokes(N-S)方程，以此对湍流进行完全解析是最为精确的。但其存在仿真耗时长，所需计算资源庞大等问题，难以高效的求解呼吸道流场[5]。若仅考虑对湍动能场起主导作用的大尺度的涡，而对小尺度的涡采用数学模型加以刻画，即大涡模拟仿真(Large Eddy Simulation, LES)，则可有效降低DNS仿真中的网格数量。JAYARAJU等人[6]和LAMBERT等人[7]通过使用LES模型对简化的上呼吸道模型仿真后发现，LES所计算的速度和湍动能场与实验数据保持较高的一致性。LIN等人[8]分别使用DNS和LES模拟分析了同一呼吸道模型的内部流场特性，结果表明LES能够较好的模拟流场中的回旋流，且与DNS得到的仿真结果高度吻合。但是，LES本质上还是一种DNS方法，其流场仍需大量网格来表征[9]，计算耗时较长。然而，从临床角度出发，医生希望能在短时间内得到呼吸道流场的仿真模拟结果，以此对呼吸系统疾病做出精准的诊断。

雷诺平均N-S方程(Reynolds-averaged Navier-Stokes， RANS)则不求解瞬时的N-S方程，转而计算时间平均湍流场。如此，流场表征所需的网格数量可显著降低，进而有效提升了流场的计算效率[10]。常用的RANS模型包括涡黏模型(两方程的k-ε模型和k-ωSST模型，以及四方程的Transition SST模型等)和雷诺应力模型(七方程的RSM-SSG模型)。ISLAM等[11]采用不同涡粘模型模拟呼吸道内部流场后发现，k-ωSST模型计算上呼吸道内的湍流时要优于标准的k-ε模型。XU等[12]分别考虑了吸气和呼气两个过程，采用不同的RANS模型得到的结果与PIV实验结果进行比较后发现Transition SST模型更适用于吸入过程的湍流场仿真，而k-ωSST模型则更适用于呼出过程的仿真。ELCNER[13]通过对简化的呼吸道模型进行仿真模拟，分析比较了k-ε、k-ω和LES三种湍流模型。结果表明LES的结果与实验结果的吻合度最高，其余湍流模型计算得到的喉部下游流动区域的流场与实验结果存在较大区别。虽然k-ε和k-ω湍流模型均能得到流场的全局特征，但难以较好的求解流速较低时的低雷诺数区域和过渡流区域。

结合Spezial-Sarka-Catski(SSG)模型[14]的雷诺应力模型(RSM)摒弃了各项同性涡粘假设，通过求解雷诺应力传输方程和耗散率方程，以封闭雷诺平均N-S方程。RSM模型比两方程模型更严格地计算了流线曲率、涡流、回旋流的变化对流场的影响，因此它在计算复杂湍流流动时具有更高优势。SOMMERFELD等人[15]通过分析对比RSM模型和RANS模型后发现，使用k-ε模型计算得到的呼吸道流场难以精准的表征湍流边界层，使用k-ωSST计算的流场能够很好的捕捉到粘性底层与湍流层之间的混合流动区域的细节，各向异性的雷诺应力(RSM-SSG)湍流模型则能更好的计算近壁面的湍流场。

此外，由于真实人体支气管树结构复杂，过往相关研究均采用混合网格的方法划分呼吸道内部流场，需要消耗巨大的计算时间和资源。然而，从临床角度出发，医生希望在48 h内给出基于CFD仿真的吸入给药效率分析结果，先前的研究并不能满足临床需求。鉴于上述矛盾，若能采用全结构化网格划分其流场，则能极大的降低网格的数量(相同尺寸的结构化网格数量仅为非结构化网格数量的1/3)，并结合更适合真实人体呼吸道内湍流仿真模型，进而可高效的提升仿真模拟计算的精度和速度。

综上所述，本研究提出使用全结构化网格划分真实人体支气管树模型，并分别采用k-ωSST模型、Transition SST模型、雷诺应力模型(RSM-SSG)和大涡模拟(LES)计算呼吸道内部流场，并与实验测量结果相对比，在确保可靠计算精度的前提下，探寻最适用于真实人体呼吸道流场的快速且精确的计算方法，为呼吸道计算流体力学仿真模拟的临床应用提供参考和依据。

1 湍流模型

1.1 Shear-Stress Transport (SST) k-ω模型

SST模型是The Baseline (BSL)k-ω模型的改进，并在湍流黏度的定义中考虑了湍流剪切应力。这使得该算法在近壁区绕流和旋流方面比标准k-ω和BSLk-ω模型更精确。其湍动能k及耗散率ω的方程如下[16]：

(1)

(2)

式中，ui——i方向的平均速度

uj——j方向的平均速度

Γk——k的有效扩散系数

Γω——ω的有效扩散系数

Yk——k因湍流引起的耗散

Yω——ω因湍流引起的耗散

Gk——k引起的湍动能

Gω——ω引起的湍动能

Dω—— 交叉扩散率

有效扩散系数公式为：

(3)

(4)

式中，σk和σω——k和ω的紊流普朗特数

μt—— 修正后的黏度系数

(5)

式中，α*—— 抑制湍流黏度的低雷诺修正系数

S—— 剪切力张量的常数项

F2—— 湍流普朗特数的融合项

α1—— 常数

1.2 Transition SST模型

(6)

转捩源项定义如下：

Pγ1=Flengthca1ρS[γFonset]cγ3

(7)

Eγ1=ce1Pγ1γ

(8)

Pγ2=ca2ρΩγFturb

(9)

Eγ2=ce2Pγ2γ

(10)

式中，S—— 应变率张量的模

Flength—— 转捩区长度

Ω —— 漩涡强度

Fonset—— 涡量雷诺数ReV的函数

Fturb—— 黏性系数系数比的函数

(11)

(12)

T —— 时间尺度

Fθt—— 由源项控制参数的混合参数

1.3 RSM模型

本研究采用的RSM模型为Spezial-Sarka-Catski的SSG模型，该模型已被证明在平面应变、旋转平面剪切和轴对称膨胀/收缩等一系列基本剪切流中具有优越的性能。典型的线性雷诺应力模型如下[20]：

(13)

Dij—— 扩散项

Pij—— 剪切应力产生项

Φij—— 压力应变项

εij—— 和黏性耗散项

1.4 LES模型

对于LES模型，针对不可压缩流动，经过滤波函数处理后的N-S方程和连续方程为[21]：

(14)

(15)

表1 四种湍流模型的优劣

2 仿真及实验设置

2.1 真实人体支气管树模型

本研究所研究的呼吸道模型是根据一名成年男子的高分辨率计算机断层扫描(CT图片)获得的呼吸道医学影像，通过医学影像三维重构软件，生成的立体光刻(STL)三维几何模型[26]，如图1a所示。三维重构后的支气管模型包括进气管、口腔、气管和支气管(7级)。气管处的水力直径约16～18 mm，第七级支气管处的水力直径约为2～3 mm，气管横截面面积变化之大和支气管树结构之复杂是呼吸道内部产生不同流动状态的根本原因。由于实验的需要，根据人体肺部五叶十八段的生理结构，末端支气管被划分为10个区域并与相应的出口连接。

图1 采用的真实人体呼吸道模型[26]

2.2 网格划分

网格划分是呼吸道仿真的重要环节之一，若网格计算单元质量低(网格非正交性和歪斜度偏高)， 则会引发数值计算发散，进而影响流场的求解精度[27]。由于口腔的模型较为复杂，此处采用非结构化网格划分。其余部位则充分利用ICEM CFD在结构化网格划分方面的优势，将支气管树及10个收集器全部采用结构化网格划分。KOULLAPIS等[9]通过对比不同网格数量(700万～5000万)对结果的影响，发现在保持较好精度的条件下，采用700万非结构化网格就能够较好的解析呼吸道内的流动。对于同一个模型，相比于非结构化网格，使用结构化网格能够大幅降低网格数量，且能更好的控制网格生成质量，使计算结果更易收敛。近壁面流动是影响颗粒沉积的主要原因[28]。为避免近壁面边界层网格厚度过低从而影响壁面压力梯度的求解精度，在网格划分时近壁面贴合的网格层数保持在15层以上(如图2所示)。该支气管树的网格总数约为1200万。

图2 采用结构化网格划分支气管树模型：结构化网格总数为900万，呼吸道壁面边界层网格层数大于15层

2.3 实验设置与边界条件

人体在不同的运动状态或健康状态下，呼吸速率和呼吸量会有较大的差别。根据吸入气体的速度是否恒定(吸入速度是否随时间变化，吸入速度不随时间变化为稳态呼吸，反之)，可以将呼吸状态分为稳态呼吸和非稳态呼吸。本研究主要研究在稳定呼吸状态下持续吸入0.5 s的过程，即Q(t) = 60 L/min，(0

表2 吸入流量Q=31.1 L/min时，10个出口的流速

表3 仿真边界条件和数值方法设置

图3 实验装置：模型的进气管和口腔完全浸没于装有水-甘油混合液体的玻璃水箱中，模型的10个出口置于水箱下部并分别与10个活塞隔膜泵连接[29]。

3 计算结果与分析

图4显示了当t=0.5 s时，采用LES算法得到的进气口到其中一个支气管出口在呼吸道冠状面和矢状面上的时间平均流速分布(仿真中每10个时间步长求一次平均作为该时间区间内的时间平均值)。其中Ⅰ，Ⅱ，Ⅲ和Ⅳ分别表示口腔和气管的矢状面、气管分叉处的冠状面、左主支气管的冠状面和右主支气管的冠状面的具体位置。线段A，B，C，D，G，H和J分别表示支气管树内矢状面或冠状面上的截线。由图4可知，随着支气管横截面面积逐渐减小，气流速度逐渐增大。每个支气管分叉处的速度变化最为剧烈，会出现较为明显的剪切层和回旋流区。为了更详细的对比不同湍流模型计算结果的差异，取t=0.5 s这一时刻点的相关云图和数据(不再赘述)，对支气管树的不同区域进行详细分析。

图4 通过LES模型计算得到的呼吸道内从进口到其中一个出口的冠状面和矢状面时间平均流速分布；A,B,C,D,G,H和J分别表示支气管树内冠状面或矢状面与相应位置处横断面的交线；Ⅰ,Ⅱ,Ⅲ和Ⅳ分别展示了口腔和气管的矢状面、气管分叉处的冠状面、左主支气管的冠状面和右主支气管的冠状面的具体位置。

图5a显示了不同RANS湍流模型和LES模型在口腔和气管矢状面上的计算得到的结果。该处支气管横截面面积较大，气体流动较慢，雷诺数常在3000～5000。从图5a中可看出，气体经过口腔(交线B处)和咽部(交线C处)时速度较低，在口腔上壁出现了流动分离现象。由于咽喉收缩，气体流经此处时速度突然增大，并伴随着剪切流和回旋流。随着气管结构产生弯曲，呼吸道内部高速流从气管后壁转移到气管前壁，并且在前壁产生了一个较小的速度分离区域，该现象在LES模型的仿真结果中较为明显。气管中段区域，LES模型仿真计算的速度从气管前壁到气管后壁过渡的较为均匀。k-ωSST模型的计算结果和RSM-SSG模型的计算结果较为一致。在Transition SST模型的仿真中，并没有在靠近气管前壁面区域产生较大的流速，流速最大处主要集中在气管中部。

对于上呼吸道内部的流场，LES模型能够捕捉到气管内更详细的剪切流和回旋流。不同RANS湍流模型得到的总体速度分布具有较好的一致性，但是在预测喉部的剪切层和回旋流时有较大的差异。k-ωSST模型与RSM-SSG模型在近壁面区域计算的结果与LES模型计算的结果吻合性较好。Transition SST的结果与LES的结果差异较大。这是由于口腔和气管中常存在回旋流和自由剪切流，Transition SST转捩流模型计算这些流动时效果不佳[32]。口腔和气管的时间平均湍动能如图5b所示，湍动能在入口处较小，随着流动的发展，在口腔的上部开始升高。最大的湍动能水平位于喉部气管前壁和喉部下游区域的气管后壁。k-ωSST模型和RSM-SSG模型明显低估了该区域的湍动能水平。

图5 使用不同RANS湍流模型和LES模型在口腔和气管矢状面上的计算结果(t=0.5 s)

为了进一步对比结果，取图4中上呼吸道中的A,B,C和D四条交线处的实验测量值、不同RANS湍流模型计算结果和LES模型计算结果进行比较。为了处理数据方便和更直观地对比各湍流模型计算结果与实验的差异，在每条交线上平均取200个离散点并编号，进而提取各个点上的时间平均速度和时间平均湍动能并进行无量纲化，最后将编号进行归一化处理，均映射到0～1范围之内的区间上，下同，如图6所示。图中D*表示交线上点的归一化坐标。umag和k分别表示仿真得到的交线上的时间平均速度和时间平均湍动能。uT=2.28 m/s表示水-甘油混合液体在Q=31.1 L/min 时，进气管横截面的时间平均速度。4种湍流模型在交线A处得到的速度分布曲线较为一致，但与实验有明显差异。原因是实验在装满水-甘油混合液体的水箱中进行的，水箱体积有限、壁面离进气口较近等因素对实验结果影响较大[15]，从而使得入口处速度不对称。气体进入口腔后，所有湍流模型在B处均计算出了剪切流。随着流动进一步发展，剪切层逐渐转移到上咽(交线C处)内侧，Transition SST模型的结果和RSM-SSG模型的结果与实验结果更为吻合。交线D处的呼吸道横截面突然增大，速度变化较为剧烈，四种模型得到的速度曲线与实验结果具有较高的一致性，但是均不能精确的捕捉到剪切层的位置，其中Transition SST的表现最为不佳。

图6 使用不同RANS湍流模型和LES模型在A,B,C和D四条交线处计算得到的无量纲化时间平均速度曲线(左图)和无量纲化时间平均湍动能曲线(右图)(t=0.5 s)；a),b),c)和d)分别表示A,B,C和D四条交线；D*表示交线上点的归一化坐标；umag和k分别表示仿真得到的交线上的时间平均速度和时间平均湍动能；uT=2.28 m/s表示水-甘油混合液体在Q=31.1 L/min时，进气管横截面的时间平均速度

与时均平均速度场相比，各个模型计算的湍动能明显不同。在进气管处，LES模型低估了湍动能水平，而k-ωSST、 Transition SST和RSM-SSG模型在进气管中心区域得到地湍动能与实验结果较吻合。其中，RSM-SSG模型在近壁面处捕捉到了较大的湍动能，这是由于仿真设置的进口边界条件为大气压力，与进口周围静态环境相比进器管口产生较大的速度变化，从而表现出更高的湍动能。在B,C和D三条交线处，LES和Transition SST模型得到的湍流强度较为接近，能够捕捉到湍动能的整体分布和大小。k-ωSST模型和RSM-SSG模型得到的湍流强度均低于实验。RSM-SSG模型计算的湍动能沿着气管衰减最快，但其仍能够在近壁面捕捉到与实验更接近的湍流强度。

随着气流进入左右支气管，气管分叉处和左右主支气管(结构特征由人体生物学定义)的时间平均速度分布云图如图7所示。在Lizal[26]的体外测量实验中，左肺和右肺的通气量之比为3∶7，这种不相等的通气量导致流速在气管分叉处的左右区域差异巨大。在气管的右壁处可以看到一个高速流动区域，此处的雷诺数能达到6000以上，而左侧分叉处出现了一个较大的回旋流区域。随着左主支气管下游变窄和变短，流速逐渐增大，该区域产生了较大的剪切流区域。右主支气管与下一级支气管分支角度较大，这种几何形状的不对称导致右主支气管内出现射流现象。

图7 使用不同RANS湍流模型和LES模型计算得到的时间平均速度场(t=0.5 s)

气管分叉处和左右主支气管的湍动能云图如8所示，尽管湍流强度在气管中出现了衰减，但是在剪切流和回旋流区域之间出现了较高的湍流水平，最大的湍流强度出现在右主支气管。LES和Transition SST模型均能够较好计算出左右主支气管处的湍动能大小，k-ωSST模型和RSM-SSG模型依旧低估了此处的湍动能水平。如上所述，RSM-SSG虽然不能准确预测呼吸道中心区域的湍流水平，但是在近壁面处能够捕捉到与实验较为吻合的湍流强度。

图8 使用不同RANS湍流模型和LES模型计算的时间平均湍动能场(t=0.5 s)

图9显示了G,H和J三条交线处的无量纲化的时间平均速度和时间平均湍动能曲线，在G处可以发现，LES模型计算的结果和实验测量结果具有更好的一致性，Transition SST模型得到的结果与LES计算的结果最为接近。k-ωSST和RSM-SSG模型在此处的仿真结果均略高于实验结果。随着气流进入左右主支气管，气管内开始出现了较为强烈的剪切流动，LES模型的计算结果与实验测量结果最为接近。湍动能方面，LES模型与Transition SST模型在G处计算得到的湍动能分布与实验较为一致，二者均能捕捉到此处湍动能的最大值和最小值。k-ωSST模型和RSM-SSG模型依旧低估了呼吸道中心处的湍动能水平，但是RSM-SSG仍然能很好的捕捉到近壁面的湍流强度。综合比较三处交线的湍动能曲线可以发现，LES计算低级支气管处的湍动能具有很大的优势，其计算的湍动能与实验测量值保持较高的一致性。k-ωSST模型在H处计算的湍动能与LES模型的结果最接近，Transition SST模型和RSM-SSG模型计算的湍动能与实验结果的吻合度较低。

图9 使用不同RANS湍流模型和LES模型在交线G,H和J三处计算的无量纲化时间平均速度曲线(左图)和无量纲化时间平均湍动能曲线(右图)(t=0.5 s)；a),b)和c)分别表示G,H和J三条截线；D*表示横坐标归一化区间；umag和k分别表示仿真得到的截线上的时间平均速度和时间平均湍动能；uT=2.28 m/s表示水-甘油混合液体在Q=31.1 L/min时，进气管横截面的时间平均速度。

4 结论

本研究使用真实人体呼吸道模型，采用全结构化网格划分支气管树，并分别使用k-ωSST模型、Transition SST、雷诺应力模型(RSM-SSG)和大涡模拟(LES)计算呼吸道流场，通过与实验测量结果比较发现：

(1) 采用结构化网格划分支气管树能够极大的降低网格数量，进而可有效提升计算效率。

(2) 在上呼吸道区域，由于上呼吸道的横截面面积较大，气流在此处的流动处于发展阶段，雷诺数相对较低。使用不同湍流模型在该区域计算得到的速度场中，Transition SST模型在流动混合区域的计算与实验结果存在差异，k-ωSST模型、雷诺应力模型(RSM-SSG)和大涡模拟(LES)的计算结果与实验结果的吻合度较高。在该区域湍动能场计算结果的比较中发现，LES不仅能较为准确地预测湍动能的分布状态，而且能捕捉到喉部下游的湍流细节；Transition SST模型虽在预测近壁面湍动能时表现欠佳，但在预测呼吸道中心区域的湍动能水平时与LES最接近；k-ωSST和RSM-SSG湍流模型的计算结果对湍动能的预测均低于实验测量结果，但RSM-SSG模型在近壁面处的湍动能计算结果与实验保持一致。由于药物颗粒的湍流扩散沉积主要与近壁面湍流场计算精度相关，因此RSM-SSG的计算结果是可接受的。

(3) 在支气管区域，由于该区域的速度变化剧烈，常伴随较强的剪切流和回旋流，雷诺数相对较高。LES和RSM-SSG在该区域的流动状态和湍动能计算结果与实验结果的差异较小。k-ωSST和Transition SST在此区域的速度场计算结果与实验测量结果较为接近，但近壁湍动能场计算结果与实验结果差异较大。

综上所述，在计算真实人体呼吸道内部流场时，LES的计算结果与实验结果最为接近。但当受限于计算效率时，结合结构化网格的巨大优势，使用RSM-SSG模型代替LES模型，可以有效节约计算时间，能够最大程度的减小计算精度与LES模型的差异对颗粒湍流扩散沉积的影响，从而为计算流体力学仿真模拟大规模应用于呼吸道疾病的临床应用提供参考和奠定基础。

致谢

作者十分感谢科技部国际合作司中国与斯洛文尼亚科技人员交流项目对该研究的支持,以及COST (European Cooperation in Science and Technology; Simulation and pharmaceutical technologies for advanced patient-tailored inhaled medicines; www.cost.eu)对本研究提供的帮助。