无心跳供体肺支气管内气体二次流动特性分析

胥 义 邓如意 刘 晶, 刘道平 赵晓刚

1(上海理工大学生物系统热科学研究所，上海 200093)2(上海理工大学制冷与低温工程研究所，上海 200093)3(上海市肺科医院胸外科，上海 200433)

无心跳供体肺支气管内气体二次流动特性分析

胥 义1*邓如意1刘 晶1,2刘道平2赵晓刚3

1(上海理工大学生物系统热科学研究所，上海 200093)2(上海理工大学制冷与低温工程研究所，上海 200093)3(上海市肺科医院胸外科，上海 200433)

无心跳供体(NHBD)肺有望解决临床供肺严重不足的问题，肺内低温通气被公认为是NHBD肺在体低温保护的有效方法之一。通过建立3D非对称四级支气管模型，借助计算流体动力学(CFD)方法，分别研究在3种被动呼吸频率(0.125、0.25、0.5 Hz)条件下肺支气管内的二次流结构特性。研究结果表明，被动吸气时主支气管截面上不存在二次流结构，且截面上最大无量纲速度差分别为0.67(@0.125 Hz)、0.5(@0.25 Hz)和0.3(@0.5 Hz)，但在被动呼气时存在明显的二次流结构，且其截面最大无量纲速度差增大到1.25、1.1和1.06；随着支气管级数的增加，无论是吸气还是呼气状态，其对应截面上都存在明显的二次流动结构，但呼吸频率对其截面上最大无量纲速度差变化的影响不大(维持在1.0左右)；总体来讲，左肺各支气管内的流动结构比右肺要复杂一些，这与左右支气管的不对称结构有关。该研究明确了低温保护性气体在NHBD肺支气管内的流动总是伴随着二次流的特性，这对于进一步探索其临床保护工艺研究具有重要的指导意义。

无心跳供体肺(NHBD)；支气管；二次流；低温通气

引言

无心跳供体(non-heart beating donor，NHBD)肺有望成为肺移植供体的重要来源[1]。尽管表面冷却(topical cooling，TC)[2-4]和肺血管低温灌注[5]均可以实现供肺的快速降温保护，但这两种方法均需要对供体进行一定的创伤处理。在未征得家属同意和完成必要的手续之前，这种操作是违反伦理的。近年来，诸多研究团队对NHBD肺进行热缺血期间的“在体”低温通气保护报以厚望，即保持肺器官在体内的原来位置状态下进行机械通气降温保护[6-10]。这需要对低温气体在肺内部的传热传质机理有充分的认识，才可以制订出合适的低温保护方案。但由于肺支气管结构的复杂性，有关肺内气体流动的实验研究难度非常大，目前主要是借助计算流体(computational fluid dynamics，CFD)技术来进行理论分析[11-12]。

二次流通常指黏性流体在弯管处形成的双旋流现象，通常存在于具有弯曲结构的流道内[13]。从解剖学结构来看，肺支气管树是由一系列“二分叉”结构组成，且分叉气管存在一定程度的弯曲。这种结构必然使流体在通过分叉处的流动过程中，由于离心力的作用，导致支气管内的气体流动产生二次流现象。从目前的文献调研来看，有关肺内气体二次流动的研究，主要集中在探讨固体颗粒物在人体呼吸系统内的传输和沉积机理，用于评估呼吸系统的药物喷雾治疗，以及污染物和有毒气体的吸入对呼吸系统损害等方面[11]。从呼吸机理来看，这些研究均是基于人体正常“主动呼吸”过程，即：呼吸肌收缩，胸腔内形成负压，使得肺泡内也形成负压，外界气体通过支气管吸入，肺泡吸气而膨胀；呼气过程则是呼吸肌舒张，胸腔形成正压，将肺泡的气体又通过支气管排出，基本能确保完全换气。但对于NHBD肺而言，由于呼吸肌不再具有收缩和舒张功能，只能通过在气管出口处形成正压(充气)和负压(抽气)来完成肺内气体交换，是一个典型的“被动呼吸”过程。这两者的初始条件和边界条件有着较大的区别，目前针对这个“被动呼吸”过程中二次流动结构的研究还鲜有文献报道。

本研究应用CFD方法，研究“被动”呼吸频率(0.125、0.25、0.5 Hz)对NHBD肺支气管三维模型内二次流动结构的影响。对于探索NHBD供肺在低温通气时降温保护效果的评价具有重要意义，有望为临床NHBD供肺低温气体保护提供一定的理论基础。

1 计算方法

1.1 几何模型

图1 支气管树三维模型(①气管；②右主支气管；③左主支气管；④右肺上叶；⑤右肺中叶；⑥右肺下叶；⑦左肺上叶；⑧左肺下叶；T-气管入口段截面；R-右肺主支气管截面；L-左肺主支气管截面；R-s-右肺上叶支气管截面；R-x-右肺下叶支气管截面；L-s-左肺上叶支气管截面；L-x-左肺下叶支气管截面；1-14为各细支气管出口编号)Fig.1 The 3D model of the bronchial tree(①Trachea；②Right primary bronchi；③Left primary bronchi；④Right uper lobe；⑤Right middle lobe；⑥Right lower lobe；⑦Left uper lobe；⑧Left lower lobe；T-cross section of trachea inet；R- cross section of right primary bronchi；L- cross section of left primary bronchi；R-s- cross section of bronchia inside right uper lobe；R-x- cross section of bronchia inside right lower lobe；L-s- cross section of bronchia inside left uper lobe；L-x- cross section of bronchia inside left lower lobe；1 to 14 indicate outlet numbers of bronchia)

依据Horsfield解剖模型[14]，建立如图1所示的前4级支气管三维实体模型。其中，T表示主支气管截面，R表示右支气管截面，L表示左支气管截面；R-s、R-x分别代表右侧上、下两个肺叶区，L-s，L-x分别代表左侧上、下肺叶区。各级管长及管径参数以及三维模型网格划分请参见参考文献[12]。笔者所研究的支气管模型均为刚性光滑圆管模型，而有关支气管黏弹性对支气管内流动结构的影响，将在后续的研究中继续深入探讨。

1.2 控制方程

由于所研究的气体在肺支气管内的流速较慢，进出口压差较小，可以近似为不可压缩流体。并且，支气管树分叉处有曲率的存在，因而其内部气体流动实际上属于弯管内的流动，必然会产生“迪恩涡”。

对于不可压缩流体，弯管内迪恩涡运动的控制方程由连续性方程(即质量守恒方程)和Navier-Stocks方程(即动量守恒方程)组成。

质量守恒方程：

(1)

动量守恒方程：

(2)

u—速度矢量；

ν—空气运动粘度；

带旋流修正的Realizablek-ε方程，其湍动能传输方程为

Gb-ρε-YM+Sk

(3)

(4)

式中，Gk—由层流速度梯度引起的湍流动能；

Gb—由浮力引起的湍流动能；

μ、μt—分别为空气动力学粘度和湍流粘度；

YM—在可压缩湍流中，过度扩散而引起的波动；

C1、C2—常量，分别取C1=1.44，C2=1.9；

σk、σε—k方程和ε方程的湍流普朗特(Prandtl)数，取σk=1.0，σε=1.2；

Sk、Sε—为源项，由用户自行定义。

一般来讲，式(4)左侧最末两项在计算过程中均省略[15]。

边界条件为：理论情况下，肺活量(vital capacity, VC)为人在完成全力吸气后用全力呼出的气体量，如图2中斜线部分的面积。正常成人男子肺活量约为4.0 L，女约为3.0 L[5]。考虑到“被动呼吸”时呼吸肌松弛，不参与呼吸过程，其名义上的肺活量会大大减小[10]，因此本研究选用的肺活量为2 L，入口段横截面流速uin随时间的变化关系为

(5)

取通气频率f分别为：高频0.5 Hz，中频0.25 Hz，低频0.125 Hz。各种通气频率下的入口段横截面最大流速uin-max分别为15.6、7.8、3.9 m/s。在图2中，进口处吸气时速度为正值，呼气时为负值，吸气和呼气时均采用0 Pa的定静压边界条件。

壁面边界条件：支气管壁面采用无滑移的刚性壁面边界条件，不考虑支气管内软骨环和黏液的作用以及壁面的弹性变形。

为了便于对比分析，在进行流速分布分析时，支气管内采用横截面的无量纲速度形式表示为

(6)

图2 3种被动呼吸频率下入口流速随时间的变化关系Fig.2 The relationship between the flow velocity and time for passive respiratory frequencies

2 计算结果

2.1 支气管树内总体流速分布

由于3种频率下的总体流速分布基本相似，这里只展示低频(即0.125 Hz)条件下的计算结果，如图3所示。可以看出，吸气时支气管内的最大速度出现在左肺下叶区域内，达到6.94 m/s；在左右主支气管内上边缘的流速较低(约1.13 m/s)，但靠近支气管下边缘的流速较高(4～5 m/s)；左右侧上叶区域内，一旦经过分叉后，管内流速要比上一级减小一些。

对于呼气过程，气体通过末端支气管流入，经过分叉处不断地合流，最后汇集到主支气管流出，其流速分布与吸气过程有着较大差别。例如，呼气时左肺上叶支气管内的流速要比吸气时高，而左肺下叶支气管内与吸气时类似，但仍然具有较高的流速。与吸气时不同的是，呼气时主支气管下端的低速区要比吸气时宽些，中心高速区相对较窄，随着气体在主支气管内不断流动，分层现象逐渐消失，流体混合越来越充分，最终从主支气管流出。

图3 0.125 Hz下吸气和呼气的速度分布。(a)吸气；(b)呼气Fig.3 Velocity distribution at 0.125Hz. (a) Inspiration; (b) Expiration

图4 主气管入口段T截面二次流及流速分布Fig.4 Secondary flow structure and velocity distribution at the T section

2.2 气管入口段(T截面)二次流及流速分布

从图4中发现，3种频率下，吸气时T截面上的速度等值线分布呈现同心圆形状；随着呼吸频率的增加，截面上的阶梯逐渐减少，其对应的最大无量纲速度差分别为0.67(@0.125 Hz)、0.5(@0.25 Hz)和0.3(@0.5 Hz)。呼气时，T截面上的速度分布形状呈现一个中心狭长的速度分布，且各种频率下的分布形状不尽相同，其对应的最大无量纲速度差分别为1.25(@0.125 Hz)、1.1(@0.25 Hz)和1.06(@0.5 Hz)。显然，呼气时的气流分布均匀性比吸气时的要差一些，但无论是“呼”还是“吸”，其气流分布的均匀性都会随着频率的增加而增强。

图5 L、R截面的二次流及流速分布Fig.5 Secondary flow structure and velocity distribution at the L and R sections

从图4中的流线分布可以看出，吸气时T截面上不存在二次流结构，但呼气时却存在明显的二次流结构。在0.125 Hz条件下，旋涡中心处于左上方，旋涡是按顺时针方向，且在左下方和右上方存在两个很小的旋涡；0.25 Hz条件下，旋涡中心也处于左上方，但是旋涡变得更大些，左下方的小旋涡也消失了；在0.5 Hz条件下，大旋涡的中心明显向左下方移动，且在右下方又形成了一个小的旋涡。

2.3 左右主支气管(L、R截面)二次流及流速分布

从图5流线分布观察，3种通气频率下，吸气时左主支气管截面L上的二次流速度分布类似，没有明显差异，都存在两个反向对称旋涡，也就是迪恩涡结构。右侧的旋涡逆时针旋转，左侧的旋涡顺时针旋转，旋涡中心均偏向上方。呼气时旋涡结构变得不规则，右上方的旋涡较小，而左下方的旋涡较大，在截面的右下方存在一个小旋涡，在这3种呼气条件下，二次流的结构类似，没有明显差异。吸气时右主支气管截面R与左主支气管截面L的二次流结构类似，存在两个反向对称旋涡，但这个旋涡要比L截面上的大；呼气时截面上也存在两个反向对称旋涡，这个对称结构向右上方倾斜。

在3种通气频率下，吸气时左主支气管截面L上的速度等值线分布都呈现向正下方凹陷的形状，高速流体大多聚集在截面底部，低速流体在截面上部流过。随着频率增大的增加，中心高速区域不断扩大，但3种频率所对应的最大无量纲速度差均为0.9左右。然而，右主支气管截面R上的速度等值线分布与左主支气管的等值线分布类似，但随着频率的增大，其所对应的最大无量纲速度差分别为1.2、1.0和0.8。

呼气时，L截面上的速度分布都呈一向右下方凹陷的形状，左上方凹槽处的流速最低，凹槽底部两边分布着两个高速区，这与吸气时高速区的分布略有不同，其所对应的最大无量纲速度差均为1.0左右。而R截面上呼气时的速度分布形状与L截面上的类似，都是向右下方凹陷的形状，只是在凹陷底部的高速区域并没有像L截面上那样被分成两块，且随着速度的增大，高速区域的面积不断扩大，其3种频率所对应的最大无量纲速度差分别为1.1、0.9和0.9。

2.4 左右侧肺上、下叶支气管(L-s、R-s、L-x、R-x截面)二次流及流速分布

图6和图7分别为左右侧肺上、下支气管截面二次流及流速分布情况。可以明显看出，吸气时L-x截面上的速度要比L-s截面上的速度大得多，最大无量纲速度要高出0.4左右。其等值线倾斜方向相反，即L-s截面向右上方倾斜，而L-x截面向左上方倾斜。呼气时这两个截面上的等值线速度分布差异较小，最大无量纲速度相差0.2左右，但是呼气时L-s截面上的速度明显要比吸气时大，且截面上端的凹槽比较突出。

从二次流分布上看，左肺肺上叶L-s和下叶L-x两个截面上的二次流分布在吸气和呼气时存在较大差异。L-s截面吸气时的左侧存在一个大旋涡结构，呼气时截面上则存在两个旋涡，其左侧的旋涡顺时针旋转，其右侧的旋涡逆时针旋转。L-x截面吸气时不存在旋涡结构，但是流体有向右上方汇集的流动趋势，呼气时截面上的旋涡结构比较复杂，不存在完整的旋涡结构，但是有向左下方汇集流动的趋势。

与左肺叶支气管内的流动相比，右肺上叶支气管截面R-s和下叶支气管截面R-x上的无量纲速度分布显得更加均匀一些，其所对应的最大无量纲速度差均为1.1左右。从二次流分布上看：吸气时，R-s截面存在一个大的顺时针方向旋转的旋涡，流体大部分汇集于右下侧； R-x截面上不存在旋涡结构，流体有向左上方汇聚的流动趋势，但有形成旋涡的趋势。呼气时，R-s截面上的二次流分布与吸气时有很大差别，截面正上方存在一个很小的涡结构，R-x截面的左侧存在两个微弱的反向对称旋涡结构。

3 讨论

图3表明，无论是在吸气还是呼气时，由于支气管结构分叉和旋转角度的存在，气体在低速区和高速区充分混合之前，流动就被下级分叉打乱，这是造成左右主支气管内的流动出现分层和二次流现象的重要原因。

图4中，吸气时T截面速度分布和二次流结构域呼气时差别很大。这是由于笔者所建立的支气管模型为光滑直圆柱型管段，不存在曲率，吸气时也就不会出现由于管段曲率引起的二次流现象，所以说此时主支气管截面上不存在二次流结构。但呼气时，由于主支气管属于出气管，其内部的流动将受到左右主支气管内气体合流的影响，在分叉处又有曲率的存在，使得从左右两侧支气管进入主气管内的流动有二次流的出现。

从图5～7可以看出，左肺L截面各支气管内的流动结构比右肺R截面要复杂一些，这与左右支气管的不对称结构有关。从解剖学上看，左主支气管长而细，右主支气管短而粗。当主支气管内的气体流动到第一个分叉以后，但是由于右主支气管短粗，分层流动明显，流体还没来得及混合，就被下一个分叉将这种分层流动分开，继续分层流动，导致右肺各支气管截面无量纲速度分布均匀。

图6 左右侧肺上叶截面(L-S, R-S)二次流及流速分布Fig.6 Secondary flow structure and velocity distribution at the L-s and R-s sections

图7 左右侧肺下叶截面(L-X,R-X)二次流及流速分布Fig.7 Secondary flow structure and velocity distribution at the L-x and R-x sections

但是，由于左主支气管长而细，进入左主支气管内的流体一开始就出现明显的分层流动现象。随着高速区的流体与低速区的流体渐渐混合，还没有完全混合时就遇到下一级分叉的出现，又被分成两束，流向左肺上叶支气管和左肺下叶支气管，最终导致左肺各支气管内的无量纲速度分布和二次流结构比较复杂。

通过以上分析讨论，可以明确气体在支气管内的流动伴随着二次流是一种必然现象。正是由于这种二次流的存在，才使得气体在支气管内更容易实现均匀混合。这很有利于认识低温气流在支气管内的传质特性，对于探索无心跳供体肺低温通气的临床保护工艺研究具有重要的启发意义。

4 结论

本课题利用数值模拟的方法，通过对NHBD供肺支气管内的气体流动特性进行数值模拟分析，得到主支气管以及有代表性支气管截面上的二次流动和速度分布。通过研究发现：

1)由于支气管的非对称结构存在，使得支气管内的气体流动必然存在较为复杂的二次流动结构。

2)在被动的“吸气”和“呼气”过程中，各级支气管内的流速分布和二次流结构均有较大的差别。

3)由于所建立的理论模型采用的是规则圆柱型支气管树，忽略了支气管实际形状、力学性能等的影响，使得本研究结论的推广有一定的局限性。因此，迫切需要建立基于临床医学影像3D重构实体模型，并结合一定的实验验证，才会使该项研究更有实际指导意义。

[1] Gomez A, David, Varela U,etal. Present state of non-heart-beating lung donation [J]. Current Opinion in Organ Transplantation, 2008, 13(6): 659-663.

[2] Masaru K, Isao M. Effective 6-hour preservation in non-heart-beating donor canine lungs with topical cooling: assessment from histopathological aspects [J]. Surgery Today, 2005, 35: 389-395.

[3] Van WC, Neyrinck AP, Geudens N,etal. Retrograde flush following topical cooling is superior to preserve the non-heart-beating donor lung [J]. European Journal of Cardio-Thoracic Surgery, 2007, 31: 1125-1132.

[4] Kutschka I, Sommer SP, Hohlfeld JM,etal.In-situtopical cooling of lung grafts: early graft function and surfactant analysis in a porcine single lung transplant model [J]. European Journal of Cardio-Thoracic Surgery, 2003, 24: 411-419.

[5] 丁嘉安, 姜格宁. 肺移植 [M]. 上海：上海科学技术出版社, 2008.

[6] Han JQ, Zhang K, Cui J,etal. Acceptable Warm Ischemia Time of Tracheal Grafts From Non-heart-beating Donors in Rats [J]. Transplantation Proceedings, 2011, 43(10): 3638-3642.

[7] Van Raemdonck DE, Jannis NC, Rega FR,etal. External cooling of warm ischemic rabbit lungs after death [J]. Ann Thoracic Surgery, 1996, 62: 331-337.

[8] Egan TM. Non-heart-beating donors in thoracic transplantation [J]. Journal of Heart Lung Transplant, 2004, 23: 3-10.

[9] Takahiro O, Alicia C, Salvatore P,etal. High-flow endobronchial cooled humidified air protects non-heart-beating donor rat lungs against warm ischemia [J]. Journal of Thoracic and Cardiovascular Surgery, 2006, 132: 413-419.

[10] 彭富裕, 胥义, 赵晓刚，等. 无心跳供体肺原位降温保护系统的研究 [J]. 制冷学报, 2012, 33(4): 74-78.

[11] Werner H, Modeling inhaled particle deposition in the human lung—a review [J]. Journal of Aerosol Science, 2011, 42: 693-724.

[12] 刘晶，胥义，刘道平, 等. 无心跳供体肺支气管内气体三维流动的数值模拟研究[J]. 中国生物医学工程学报, 2014,33(3)：320-328.

[13] Nadim N, Chandratilleke TT. Secondary flow structure and thermal behaviour of immiscible two-phase fluid flow in curved channels[J]. International Journal of Thermal Sciences, 2014, 82: 9-22.

[14] Horsfield K, Dart G, Olson DE. Models of the human bronchial tree [J]. Journal of Applied physiology, 1971, 31(2): 207-217.

[15] Lateb M, Masson C, Stathopoulos T,etal. Comparison of various types of k-εmodels for pollutant emissions around a two-building conguration[J]. Wind Eng Ind Aerodyn, 2013, 115:9-21.

Analysis on Secondary Flow Characteristics inside Bronchia of Non-Heart Beating Donor Lung

Xu Yi1*Deng Ruyi1Liu Jing1,2Liu Daoping2Zhao Xiaogang3

1(InstituteofBiothermalScienceandTechnology,UniversityofShanghaiforScienceandTechnology,Shanghai200093,China)2(InstituteofRefrigerationandCryogenicEngineering,UniversityofShanghaiforScienceandTechnology，Shanghai200093,China)3(ShanghaiPulmonaryHospital,Shanghai200433,China)

Non-heart beating donor (NHBD) lung is expected to be the most potential donor source for clinical lung transplantation, and low temperature ventilation was normally thought to effectively protect the donor lung in vivo. A 3D asymmetric bronchial model with four levels was reconstructed, and computational fluid dynamics (CFD) technique was adopted here for investigating the sencond flow structures inside NHBD lung bronchia under three kind respiratory frequencies, i.e. 0.125 Hz,0.25 Hz and 0.5 Hz. The results show that no second flow structures exist inside main bronchia during inspiration process, and the maximum difference of dimensionless velocity across the main section is 0.67@0.125 Hz, 0.5@0.25 Hz and 0.3 @0.5 Hz, respectively. However, there are obvious differences for the expiration process, and the maximum difference of dimensionless velocity across the main section is 1.25, 1.1 and 1.06, respectively. Due to the asymmetric structure of bronchial and bifurcation rotation angles, the second flow structures always present for further daughter bronchias but the maximum difference of dimensionless velocity of them do not change significantly with increasing breathing rates (about 1.0). Generally, the flow structures inside left side are more complicated than the right side. The current work indicates that the second flow structures always exist when low temperature ventilation is utilized to pretect NHBD lung, which may be beneficial to develop the optimization method for preserving NHBD lung in situ.

non-heart beating donor (NHBD); bronchia; secondary flow structure; low temperature ventilation

10.3969/j.issn.0258-8021. 2015. 04.007

2014-11-25， 录用日期:2015-06-10

国家自然科学基金(50906056)；上海市自然科学基金(13ZR1428600)

R318

A

0258-8021(2015) 04-0429-09

# 中国生物医学工程学会会员(Member, Chinese Society of Biomedical Engineering)

*通信作者 (Corresponding author)，E-mail: xu_hongyi@263.net