孙 栋 李福生 徐新喜 赵秀国 谭树林

（军事医学科学院卫生装备研究所 国家生物防护装备工程技术研究中心，天津 300161）

引言

上呼吸道是呼吸系统的重要组成部分，是人体与外界环境进行气体交换的主要通道。随着科学技术的迅猛发展以及大气环境的不断恶化，这条与外界物质交换最为频繁的通道越来越多地引起人们的关注。近年来，国内外许多学者从生物流体力学的角度对其进行了大量研究。徐新喜等对循环呼吸模式下人体口喉模型内气流组织特性进行了数值仿真和实验研究［1］；Wang等对上呼吸道内的呼吸气流组织形式进行了数值分析［2］，Nithiarasu等利用无结构网格对上呼吸道中的层流和湍流计算进行了数值模拟［3］，赵秀国等对包括前3级支气管在内的人体上呼吸道模型内的气流运动特性及气溶胶沉积进行了PIV实验研究［4－5］。目前的研究大多假定呼吸道壁面为刚性壁面，而实际上呼吸道内的真实呼吸流是典型的流体-结构交互作用下的流动，这一流固耦合过程对呼吸道内的气流组织形式、气溶胶扩散沉积以及壁面剪切应力的分布有着重要的影响。

流固耦合力学模型的发展和成熟，为研究人体上呼吸道内的气流运动提供了一种崭新的思路。利用流固耦合力学模型，能较好地分析气流和呼吸道壁面的相互耦合作用，避免以往单纯假设呼吸道壁面为刚性壁面，从而使理论数值模拟分析与实际呼吸道内的流动情况更为接近，拓展了人们对呼吸道内气流运动的研究思路和研究方法。本研究通过构建流固耦合力学模型，对流固耦合作用下包括前三级支气管在内的上呼吸道内稳态气流运动特性进行了数值模拟，系统分析了流固耦合作用下上呼吸道壁面的形变特点、壁面剪切应力分布以及呼吸道内的气流运动特点。

1 人体上呼吸道模型和数值仿真方法

1.1 人体上呼吸道模型

研究所采用的上呼吸道模型主要由两部分组成：口喉模型，包括口腔、咽、喉和气管；前3级支气管模型，包括G0～G3级支气管模型。口喉模型的尺寸参考了ARLA（Aerosol Research Laboratory of Alberta）模型［6］和 Stapleton K W 模型［7］的尺寸，前 3级支气管模型采用了Weibel的支气管模型［8］，如图1所示。

1.2 控制方程

图1 人体上呼吸道模型Fig.1 Human upper respiratory tract model

1.2.1 流体控制方程上呼吸道中流体湍流模型采用耦合了Langtry-Menter转捩模型的 Menter's SST两方程模型。Langtry-Menter转捩模型是Menter与Langtry提出的严格建立在局部变量基础上的基于经验关联式的转捩模型，因而能够与CFD代码耦合。该模型基于两个输运方程，即间歇因子的输运方程和转捩动量厚度雷诺数输运方程。前者与SST湍流模型耦合，用以开启转捩点下游的湍动能产生项，从而触发当地的转捩；后者用来捕捉非当地湍流强度的影响，且能够将经验关联式与间歇因子输运方程中的转捩起始准则联系起来，得到张量表达式［9－10］。

间歇因子γ的输运方程为

式中，生成项Pγ和耗散项Eγ的表达式分别为

当地转捩动量厚度雷诺数～Reθt输运方程为

式中，源项 Pθt定义为

式中，U为当地速度，混合函数Fθt的表达式为式中，θ?BL为边界层动量损失厚度，δ为边界层厚度，混合函数Fθt用来关掉边界层内的生成项Pθt从而使得从边界层以外扩散得到，Fwake的作用是保证混合函数Fθt在尾迹区不被激活。

1.2.2 壁面（固体）结构控制方程

将上呼吸道壁材料假设为线弹性材料，而非传统的刚性壁面假设，这样考虑可以使其与实际气流情况更为接近。根据文献［11］，设定气管为均匀各项同性的材料，即服从胡克定律的薄壁弹性模型，其杨氏弹性模量为9 MPa，泊松比为0.4。

固体模型控制方程可表达为式中，σs为是气管壁应力张量，ρs为气管壁密度，as为气管壁加速度。

在流固耦合交界面上，应当满足

赏识教育作为幼儿教育中的首要教学方法，是在幼儿教育过程中通过鼓励性的话语和肢体语言启发幼儿，让幼儿在受教育中得到更多的爱与关心。不仅能够帮助幼儿建立自信心，而且有助于幼儿心智全面健康发展。以下分三个方面对赏识教育在幼儿教育中的应用实践进行探索。

式中，d为位移，n为边界法向，下标s和 f分别表示固体和流体。

1.3 数值仿真方法

流固耦合问题一般分为两类，一类为流-固单向耦合，一类为流-固双向耦合。本研究对流固耦合采取双向耦合的方式，其基本思路是：分别计算流场和固体结构，然后通过中间平台来交换耦合量。在每次大迭代中，分别进行1次流体和固体计算，并交换两次数据（每个方向各1次），直到最终收敛。其中，在固体结构计算上，由于人体上呼吸道的生理结构极其复杂，咽喉-气管部位又受到周围组织的非线性、非固定的力学约束，仿真非常困难。为简化起见，模型通过合理设置壁面的弹性系数以及在模型的口部添加固定位移约束，以达到仿真的目的。而在流场的计算上，设置其边界条件为：在低强度呼吸条件（呼吸流量30 L/imn）下，口喉模型的入口边界采用速度入口条件，入口边界假定速度均匀分布，支气管出口采用压力出口边界条件，出口相对压力为0。

1.4 数值仿真模型验证

应用粒子图像速度仪（particleimage velocimetry，PIV），对人体在呼吸流量为30 L/min的低强度呼吸条件下上呼吸道内的稳态气流运动特性进行了试验研究；将试验结果与数值仿真结果进行对比，验证数值仿真方法的准确性。

图2 人体上呼吸道喉部气流组织形式试验结果与数值仿真结果对比。（a）试验结果；（b）仿真结果Fig.2 Comparison of airflow patterns between experimental results and large eddy simulation results of the larynx.（a）Experimental results；（b）Large eddy simulation results

图2为上呼吸道喉部的试验结果与数值仿真结果的对比。图2（a）为PIV试验各点气流速度的测量结果，图2（b）是仿真气流的速度矢量。可以看出，两者均是在气流进入咽部时，截面积变小而气流速度增大，整体气流组织形式较为一致。同时试验速度最大为8.51 m/s，仿真速度最大为7.85 m/s，且测量数值与仿真数值误差最大不超过10％，吻合较好，从而说明数值仿真方法的准确性。

2 结果与讨论

2.1 上呼吸道形变分析

图3为上呼吸道壁面的位移变化云图。图3（a）为上呼吸道壁面的总体形变对比，显示了人体上呼吸道壁面的总体位移特点。其中，总体位移是指上呼吸道壁面上各点受到气流作用后相对初始时刻位置的变化，由3个方向的位移组成：前后位移dx、左右位移dy、纵向位移dz，即d=dx+dy+dz。图3（b）显示了上呼吸道壁面纵向形变特点，图3（c）和图3（d）显示了上呼吸道壁面的前后形变特点。通过图3（a）可以看出，在呼吸流量为30 L/min的稳态气流作用下，人体上呼吸道整体向后运动，咽喉部位向后的位移量较小，气管支气管向后的位移量较大，且气管支气管的位移量由上往下呈逐级递增的状态，气管前壁面受到拉伸，后壁面受到压缩。由图3（c）和图3（d）可以看出，咽喉部位受到气流运动的作用，造成前壁面的中间前后位移小、两边前后位移大，后壁面的中间前后位移大、两边前后位移小，造成前壁凹、后壁凸的现象，在咽部的后壁面中心处形成一个小凸包。通过图3（b）可以看出，在稳态气流的作用下，上呼吸道前壁面向下运动、后壁面向上运动，且咽喉部位后壁面向上运动的位移量最大，而气管前壁面向下的位移量最大。

图3 呼吸流量为30 L/min时上呼吸道壁面形变。（a）形变前后对比；（b）纵向位移；（c）和（d）前后位移Fig.3 The deformation of human upper respiratory tract model with breathing intensity Q=30 L/min.（a）Comparison of pre-and post-deformation；（b）Longitudinal displacement；（c）and（d）Lateral displacement

2.2 上呼吸道壁面剪应力分布分析

图4 呼吸流量为30 L/min时上呼吸道壁面的Von Mises stress分布。（a）口喉模型前壁正视图；（b）口喉模型后壁正视图；（c）气管支气管前壁正视图；（d）气管支气管后壁正视图Fig.4 The Von Mises stress distribution of human upper respiratory tract model with breathing intensity Q=30 L/min.（a）Front view of anterior wall in human mouth-throat model；（b）Front view of posterior wall in human mouth-throat model；（c）Front view of anterior wall in trachea to triple bifurcation region；（d）Front view of posterior wall in trachea to triple bifurcation region

图4为人体上呼吸道壁面剪应力分布云图。其中，（a）和（b）为口喉模型上剪应力的分布状况，（c）和（d）为气管支气管上剪应力的分布状态。可以看出，口喉模型中受到的壁面剪切应力较大，气管支气管受到的壁面剪切应力较小，剪切应力最大处可达30.34 Pa，且气管支气管所受壁面剪应力由上往下呈现出逐级递减的状态。此外，口喉模型前壁面所受剪切应力小于后壁面所剪应力，气管前壁面所受剪切应力大于气管后壁面所受剪切应力。通过图4中的（a）和（b）可以看出，气流进入咽部时，在咽部的后壁面中心处以及咽部两侧形成剪应力集中区，剪应力较大，而咽部的前壁面所受剪应力较小。由于喉部结构形状复杂，后壁面中心处及两侧产生应力集中区；同时，由于声门处的截面减小和气流转向，在声门上方的后壁面上产生两个接近对称的剪应力集中区。通过图4中的（c）和（d）可以看出，气流到达声门后，在声门两侧产生较大的剪应力，在中间部位产生较小的剪应力，由于气流在声门处向气管的前壁面产生喷射，造成气管前壁面上的剪应力明显高于后壁面上的剪切应力。在气管、支气管的分叉处出现剪切应力集中区，该处的剪切应力明显大于其他部位。

2.3 上呼吸道内流场分析

图5为呼吸流量30 L/min时人体上呼吸道模型内不同截面的轴向速度分布。其中，图5（a）为口喉模型内对称垂直截面内的气流流场分布，图5（b）为气管支气管对称截面内的气流速度分布。可以看出，在口喉模型内，气流分别在咽部和喉部形成两个速度增长点，并且在咽部产生分离现象，在咽部外形成分离区，见图5（c）中C-C＇截面。气流在声门位置受到几何结构的限制，产生湍流喷射的现象。声门下方的气管内发生流动分离，靠近气管内壁气流速度较高，见图5（c）中 D-D＇和E-E＇截面。随着与声门距离的增加，气管内、外壁的气流速度差逐渐减少，但是气管内壁的气流速度还是要高于外壁的气流速度，如图5（c）中F-F＇截面所示。气流在支气管分叉处发生分离，靠近支气管内壁速度较高，靠近外壁速度较低，见图5（c）中 H-H＇截面。在 G2和 G3级支气管内，产生了非对称流量的气流流动现象。

图5 呼吸流量为30 L/min时上呼吸道不同截面内的速度分布。（a）口喉垂直截面；（b）气管支气管垂直截面；（c）横截面Fig.5 Axial velocity contours at different cross section in human upper respiratory tract model with breathing intensity Q=30 L/min.（a）The vertical section of human mouth-throat model；（b）The vertical section of trachea to triple bifurcation region；（c）The cross sections

3 结论

1）在呼吸流量为30 L/min的稳态气流作用下，人体的上呼吸道整体向后运动，咽喉部位向后的位移量较小，气管、支气管向后的位移量较大，且气管、支气管的位移量由上往下呈逐级递增的状态，气管前壁面受到拉伸，后壁面受到压缩。

2）口喉模型中受到的壁切面剪应力较大，气管、支气管受到的壁面剪应力较小，且气管、支气管所受壁面剪切应力由上往下呈现出逐级递减的状态；口喉模型前壁面所受剪切应力小于后壁面所受剪切应力，气管前壁面所受剪切应力大于气管后壁面所受剪切应力。

3）在口喉模型内，气流分别在咽部和喉部形成两个速度增长点，并且在咽部发生分离现象；气流在声门位置产生湍流喷射的现象，声门下方的气管内发生流动分离，靠近气管内壁气流速度较高；气流在支气管分叉处发生分离，靠近支气管内壁速度较高。

［1］徐新喜，赵秀国，谭树林，等.循环呼吸模式口喉模型内气流运动特性数值模拟［J］.力学学报，2010，42（2）：182 －190.

［2］Wang Ying， LiuYingxi， SunXiuzhen， etal. Numerical analysis of respiratory flow patterns within human upper airway［J］.Acta Mecheanica Sinica，2009，25（6）：737 －746.

［3］Nithiarasu P，Liu CB，Massarotti N.Laminar and turbulent flow calculations through a modelhuman upperairways using unstructured meshes［J］.Communications in Numerical Methods Engineering，2007，23（12）：1057－1069.

［4］赵秀国，徐新喜，谭树林，等.人体上呼吸道内稳态气流运动特性的PIV初步试验研究［J］.实验流体力学，2009，23（4）：60－64.

［5］徐新喜，赵秀国，谭树林，等.人体上呼吸道内气流运动特性的数值模拟分析［J］.计算力学学报，2010，27（5）：881－886.

［6］Grgic B，Finlay WH，Heenan AF.Regional aerosol deposition and flow measurements in an idealized mouth and throat［J］.Aerosol Science，2004，35：21 －32.

［7］Stapleton KW，Guentsch E，Hoskinson MK，et al.On the suitability of k-εturbulence modeling for aerosol deposition in the mouth and throat：a comparison with experiment［J］.Aerosol Science，2000，31（6）：739 －749.

［8］Weibel ER.Morphometry of the human lung［M］.New York：Academic Press，1963：1 － 3.

［9］Menter FR，Langtry RB，Likki SR，et al.A correlation based transition model using local variables Part 1-Model Formulation［R］.ASME 2004-GT-53452.

［10］Langtry RB，Menter FR，Likki SR，et al.A correlation based transition model using local variables Part 2-Test Cases and Industrial Applications［R］.ASME 2004-GT-53434.

［11］Wlofgang AW，Timon R.Fluid-structure interaction in lower airways of CT-based lung geometries［J］.International Journal for Numerical Methods in Fluids，2008，57（5）：653 －675.