机械通气是全身麻醉期间维持呼吸以及危重症患者进行生命支持不可或缺的治疗手段。机械通气肺损伤分为机械性损伤（气压伤、容量伤、不张伤）和生物伤。早期以机械性损伤为主，而后随着炎症反应的发展以生物伤为主^[1]。目前国内外的研究集中在生物伤的病理生理机制^[2]。由于肺内支气管气道结构复杂，难以准确监测肺内气体流动，增加了机械性损伤的研究难度。通过影像学DICOM格式图像，使用计算机三维重建三维气道模型，对模型进行网格划分后使用计算流体力学(computational fluid dynamics，CFD)技术可对气道内气流流动的各类问题进行数值实验、模拟和分析研究，在辅助诊断疾病、治疗方式选择、疗效评估方面有良好的临床应用前景^[3-5]。然而, CFD技术用于模拟仿真呼吸道气体流动的研究还存在很多争论。目前呼吸道的CFD模拟通常采用光滑的圆柱管道形成单代分叉^[6]或多代分叉^[7]，研究表明几何模型形态结构的精确度对于有效预测气流流动和颗粒分布非常重要^[8]。另一关键点在于呼吸道模型边界条件的假设，通常CFD气道模型入口边界被假设为一个恒定或者循环性的速度或压力^[9]，而个性化呼吸道模型的建立和边界条件的设置是准确获得计算结果的关键^[10]。基于影像数据提取的气道结构模型，利用CFD技术施加呼吸曲线，可以定性和定量观察整个呼吸周期内气道树中的流速、壁面压力、壁面剪切应力，也可以计算压降、研究各关键截面的垂直速度和二次流^[11-12]。MRI能良好识别软组织及细微结构，精确还原人体气道，实现气道可视化^[13]。我们通过MRI图像建立人体气道三维模型，在既往研究的基础上，利用CFD技术模拟机械通气下气流流场特性，为研究机械性肺损伤打下基础。

资料和方法

研究对象  本研究已获上海交通大学附属第一人民医院伦理委员会批准（2018KY605），选取1名健康的中年男性志愿者，签署知情同意书。该志愿者为本院麻醉科医师，年龄45岁，身高173 cm, 体重65 kg, ASAⅡ级，无气道结构异常，无MRI检查禁忌证。

设备和方法

气道模型三维重建及网格划分  取仰卧位，使用Systems/Ingenia 3.0T磁共振成像系统（Philips Medical）进行MRI扫描。扫描范围：头颈部加胸部连续扫描。扫描参数设置为3D压脂横断面T1成像；重复时间2.2 ms, 回波时间4.7 ms, 层厚4 mm。将MRI图像以DICOM格式导入Mimics 17.0医学图像软件(比利时Materialise公司)，利用气道分割功能进行图像分割，标记出气道区域，运用区域增长的方式生成相应的蒙板图像，此时得到的蒙板图像较为粗糙，需使用编辑蒙板功能，对每一帧MRI图像中所生成的蒙板进行精细修改。三维模型重建时，需根据不同组织结构的特点来选取不同的图像分割方法，图像分割后进行去噪、平滑、补洞、去除伪影等进一步修饰，最后经过三维计算工具生成高质量的三维可视化模型。将Mimics软件导出的气道三维模型STL文件导入Geomagic Studio 12软件（美国Geomagic公司）进行修整。再将修改后的模型导回Mimics软件，利用3-matic工具进行体网格划分并导出。然后将气道网格模型导入Ansys Workbench软件（美国Ansys公司）, 利用ICEM CFD模块进行网格质量检查、重新划分网格并标记出入口。综合考虑计算精度和时间，采用非结构化四面体划分网格，单元数为50万，节点数为9万。

数学模型及控制方程  由于气道结构的复杂情况以及流体流动的特性，对气道内气流进行如下假设：假设气道内空气为连续的理想气体，温度保持37 ℃恒定不变，因此可忽略能量方程；假设气道内部气流为稳态流动，因此无需求解动量方程中的时间相关项。根据公式计算出本研究模型马赫数 < 0.8，雷诺数 > 2 300，符合亚音速流体的湍流流动。CFD模拟中常用的湍流模型包括单方程模型、大涡流模型、双方程模型（标准k-ε模型、重整化群k-ε模型、可实现的k-ε模型）和雷诺数应力模型等，其中标准k-ε湍流模型应用最广泛，且计算时间较短，结果相对稳定、精确，能够满足本研究的需求。以下为本研究相关的计算公式及控制方程：

马赫数($ Ma$)计算公式：$ Ma = \frac{v}{c}$。$ v$表示流体流速，$ c$表示当地音速。

雷诺数($ Re$)计算公式：$ Re = \frac{{\rho vd}}{\mu }$。$ \rho$表示空气密度，$ v$表示空气最大流速，$ d$表示模型入口直径，$\mu $表示空气动力黏度，温度为37 ℃。

连续方程：$ \frac{{\partial {u_j}}}{{\partial {x_j}}} = 0$。

动量方程：$ \frac{{\partial \left( {{u_i}{u_j}} \right)}}{{\partial {x_j}}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial {x_i}}} + \frac{\mu }{\rho }\frac{{{\partial ^2}{u_i}}}{{\partial {x_i}\partial {x_i}}}$。

湍动能k方程：$ \frac{{\partial \left( {\rho k} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho k{u_i}} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] + {G_k} + {G_b} - \rho \varepsilon - {Y_M}$。

湍动能耗散率ε方程：$ \frac{{\partial \left( {\rho \varepsilon } \right)}}{{\partial t}} + \frac{{\partial \left( {\rho k{u_i}} \right)}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _\varepsilon }}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \right] + {C_{1\varepsilon }}\frac{\varepsilon }{k}\left( {{G_k} + {C_{3\varepsilon }}{G_b}} \right) - {C_{2\varepsilon }}\frac{{\rho {\varepsilon ^2}}}{k}$。

湍流黏度系数：$ {\mu _t} = \rho {C_\mu }\frac{{{k^2}}}{\varepsilon }$。

方程中$ {u_i}, {u_j}$(i, j=1, 2, 3)表示x、y、z方向的速度矢量，$ p$表示压力，$ G_k$表示由平均速度梯度而产生的湍流动能，$ G_b$表示由浮力产生的湍动能，$ Y_M$表示可压缩湍流扩散产生的波动，$ {\sigma _k}$和$ {\sigma _\varepsilon }$分别是与k方程和ε方程对应的湍流普朗特数，$ {C_\mu }$、$ {C_{1\varepsilon }}$、$ {C_{2\varepsilon }}$、$ {C_{3\varepsilon }}$均为常数。

CFD模拟计算  进入Fluid flow (CFX)模块，选择标准k-ε模型进行模拟计算，分析不同气流流速下的气流流场特性。CFX求解器设置：离散格式选择高精度，最大迭代次数设置为100，收敛残差标准 < 10E-4，时间步设为默认自动时间步。设置各种边界条件：入口边界条件仿照机械通气吸气相，分别设置入口流速为1、2、5和10 m/s；出口边界条件采用0 Pa的定静压边界条件，即1个大气压；壁面边界条件假设气道壁面为无滑移的光滑壁面，忽略气管环以及气管壁的弹性变形。通过CFX求解器进行模拟计算后，再利用后处理器完成计算结果的统计和图形化处理，创建各时间步的流速图、流线图、管壁压力图、管壁切应力图及动画演示等。

结果

本研究通过气道MRI影像资料对气道模型进行三维重建（图 1），并使用CFD模拟计算机械通气吸气相的气流流场变化（图 2~4）。

气道流线图（图 2A）  入口流速1 m/s气道流线图整体上较为顺畅光滑，没有发现明显的涡流，对应的流量为6 L/min，即我们平时临床工作中最常见的分钟通气量，可将其作为对照与其他流速状态下的气道情况进行对比。其他入口流速情况下的流线分布较为相似。在分叉处附近存在曲率变化较大的区域，甚至出现回旋线，表明出现涡流。与1 m/s相比，其余3种入口流速的气流在气管和两侧支气管分叉附近的流线曲率变化较其他地方更为明显，随着入口流速增大，分叉处回旋线增多，而且从动画上看，回旋线流速明显减慢甚至停滞。

气道涡流体绘制图（图 2B）  分叉处附近存在曲率变化较大的区域，甚至出现回旋线，表明出现涡流。1、2 m/s的入口处存在涡流，分叉处的涡流相对较少，而5、10 m/s的分叉处涡流较多。

气管壁剪切应力云图（图 2C）  随着流速增大，壁剪切应力逐渐增大。狭窄处会出现高壁剪切应力，剪切应力随着狭窄程度的增大而增加。

气道流速云图（图 3）  各种入口流速条件下的气道流速分布大体一致，从入口到出口的流速逐渐降低。我们发现在入口处有一高流速区域。在气道分叉附近，有个别区域比周边区域的流速明显降低，甚至接近于低限。1 m/s的流速降低区域在分叉处，2、5、10 m/s的流速降低区域分布在分叉周边靠近管壁的区域。

气道壁压力云图（图 4）  气道壁压力由气管到支气管逐渐降低。流速越高，壁压力的绝对值越高。在气道分叉处，壁压力高于其他区域。

讨论

既往研究表明，利用CFD技术能够揭示正常生理病理条件下的肺内气流特性。Luo等^[14]采用低雷诺数k-ω湍流模型计算模拟出从喉部到气管支气管的气流情况。Brouns等^[8]通过对气管狭窄的CFD模拟研究发现气管狭窄的流量和压力符合幂指数函数关系。Freitas等^[15]基于包括六级支气管的几何模型，采用Lattice-Boltzmann方法，对吸气和呼气过程的气流进行了模拟，并且应用颗粒影像速度的方法对仿真结果进行了实验验证。Cebral等^[16]利用虚拟支气管镜法以及CFD模拟的方法，对获取的Wegener肉芽肿病患者的气管支气管模型进行模拟计算，发现在气管狭窄处存在低压力和高壁面剪切应力的现象。Mihaescu等^[17]利用CFD技术模拟研究声门下狭窄患儿在吸入空气和氦氧混合气时的气流分布。我们前期应用CFD技术对不同麻醉插管头位下上气道气流流场特性进行了比较研究^[18]，证实了CFD应用于麻醉相关呼吸动力学研究的可行性。本研究拟通过MRI图像建立气道三维模型，应用CFD技术模拟分析机械通气吸气相气流流速变化所导致的气流流场特性改变，为进一步探索气道的病理生理以及各种通气模式下气道结构和气流流场变化打下基础。

我们发现不同流速下气道流线图和气道流速云图分布较为相似。气道模型从入口到出口的横截面积逐渐增大，因而从入口到出口气道流速会逐渐降低。边界层中流速较低的流线会遮挡中心区域流速快的流线，在气道分叉处出现回旋线，即涡流。入口流速越大，分叉处涡流越明显。从气道涡流体绘制图可看出，流速较低时入口处易形成涡流，而流速高时在气道分叉处易形成涡流。气道入口处的高流速区域可能与入口气流垂直方向集中进入气道有关，此处气道结构略有曲度，从而导致入口气流刚进入气道即触壁。气道分叉处有个别区域比周边区域流速明显降低，流速越低，降速区域越靠近分叉处，流速越高，降速区域分布在气道分叉周边靠近管壁。说明低流速时，气流缓慢匀速进入气道分叉，气体流通性能高，而在高流速时，分叉处对气流有部分反射，从而在管壁附近形成涡流。我们发现，气道流体流速分布与入口流体流速成正相关，并且随着气道结构的深入，呈梯度下降。观察气道横断面与纵切面图可见高速气流集中在靠近气道前壁的区域，而靠近气道后壁的区域流速相对较低。

比较不同入口流速下气管壁剪切应力云图可发现，对于固定的气道结构，当流速增大时，相当于气道相对变窄，壁剪切应力也逐渐增大。壁剪切应力是气流正切于气道内壁的力。对于正常形态的气道，壁剪切应力非常小，然而对于存在病变的气道，气流可能会导致很高的壁剪切应力。与涡流不同的是，高壁剪切应力出现在气道入口附近，这可能由于入口垂直气流对管壁的冲击以及入口处相对的狭窄增大了进气气流的流动阻力。狭窄处会出现高壁剪切应力，剪切应力随着狭窄程度的增大而增加，这与以往的研究结果相同^[8]。

比较气道壁压力云图可以发现，气道壁压力与入口流速呈正相关，入口流速快，则模型内流体流速高，对管壁的压力大；反之则流速低，对管壁的压力小。壁压力分布均匀，最大的壁压力出现在受到气流直接冲击的气道分叉处。气道分叉处的壁压力高于其他地方，这可能与气流的直接冲击有关，高壁压力和气流冲击可能会导致炎症产生^[19]。

本研究的局限和不足为：(1)由于MRI建立三维气道模型CFD模拟机械通气的研究还处于探索阶段，目前只对1例健康志愿者进行计算流体力学模拟，未来将进一步建立气道模型数据库，加大样本研究；(2)模拟结果与既往文献观察结果类似，但未建立气道实体模型来校验模拟结果的准确度；(3)提取的气道模型仅包含气管和支气管，未来将建立更多级数的气管树模型；(4)求解过程采用的是稳态计算，今后将进行包含整个吸气和呼气过程的瞬态模拟研究。

已有学者通过气道的三维重建以及CFD模拟，研究机械通气对气道壁剪切应力和应变率的变化，来改善重症患者的机械通气策略^[19]。本研究利用CFD技术对人体气道气流流动情况进行模拟研究发现，气道流速、壁压力、壁剪切应力以及涡流的形成都与入口流速以及气道结构有着密切的关系，有助于我们从形态学和力学的角度，重新认识气道功能与解剖结构之间的相互关系，在以后的研究中，我们拟深入探索机械通气性肺损伤发生发展的生理和病理基础；同时，对各种类型的机械通气（气管导管、气切导管、喉罩、双腔管等不同的通气工具，压力控制通气、容量控制通气等不同的通气模式^[20]，各种手术麻醉状态下的机械通气）进行模拟计算，比较不同的通气方案的优劣，从而提高机械通气的安全性。