不同因素对软腭被动运动的影响

王澍，姜成惠，姚佳烽，朱桂平

(1. 南京航空航天大学 a. 航天学院； b. 机电学院，江苏 南京 210016；2. 江苏省口腔医院，江苏 南京 210029)

0 引言

腭咽闭合不全和阻塞性睡眠呼吸暂停低通气综合征(OSAHS)都是由于腭咽闭合障碍导致的呼吸道疾病，腭咽闭合障碍会严重影响患者的生活质量和身心健康[1]。腭咽闭合功能受到很多的因素影响，一般认为上气道组织结构异常是导致患者的腭咽闭合障碍的重要原因[2]。上气道结构异常导致的气道狭窄是引发OSAHS的病理学基础，而软腭结构缺陷则会导致腭咽闭合不全。除此之外，神经障碍和学习型障碍也会导致腭咽闭合障碍。

在早期，人们通过从尸体铸造而成的模型或者动物实验的方式来研究人体上呼吸道的流动特性[3]。但是，尸源性研究受组织离体的特性变化影响较大，而动物的气道结构和人类的呼吸道结构也有较大差异，故而以上研究方法仍然有较大的局限性。此后，影像技术得到了巨大的发展，研究人员开始使用CT和MRI技术扫描得到人体在正常生理状态下的上气道结构数据，并构建对应的实体模型用于相应的实验研究[4]。这种实验的针对性较强，一般只能用于测试一种参数。

近些年来，随着高性能计算机的出现以及数值分析软件的不断成熟，计算流体力学被广泛应用于上呼吸道的流场模拟，一般的研究流程为先通过CT或MRI等图像扫描技术获取研究对象的气道模型，再通过仿真软件计算所需要的数据。于驰[5]分别建立了健康人与OSAHS患者的咽腔模型，对比了两者呼吸道压力以及软腭的位移状况。CHOULY F等[6]建立了简化的上气道模型，并使用流固耦合的方法模拟了在呼气过程中舌头和咽腔壁与气流的相互关系，研究OSAHS患者打鼾现象的产生原因。使用数值模拟分析的方法可以有效地克服呼吸道内表面复杂，空间狭小的限制。

1 模型的构建与仿真

1.1 上呼吸道三维模型的构建

对健康无呼吸道疾病的人体在静息状态下的上呼吸道进行CT扫描，得到了人体从声门到嘴唇和前鼻孔的三维模型，如图1(a)所示。为了更好地进行相应的流固耦合计算，通过测量上呼吸道关键部位的尺寸，重建了简化的上呼吸道模型，如图1(b)所示。模型的流体域由口腔、鼻腔和咽腔组成，模型的固体域为软腭的模型。

图1 上呼吸道三维模型重建

1.2 仿真模型

1)流体仿真模型

计算所使用的流体介质为常温下的空气，密度ρ=1.225 kg/m3，动力黏性系数μ=1.8×10-5kg/ms。在流体仿真过程中，希望通过计算得到在某一特定状态下的上气道压力分布，所以忽略上气道和软腭的变形以及相应肌肉的调节变化，将气道壁和软腭视作不可变形的刚性体，将空气视作不可压缩的流体且忽略呼吸过程中的温度变化。所以通过上气道的空气要满足连续性方程和Navier-Stokes方程，可表示为

(1)

(2)

式中：ρ为密度；t为时间；ux、uy、uz分别为速度在3个方向上的分量；V为速度矢量；p为动水压强；μ为动力黏性系数；f为单位质量的质量力。

雷诺数的计算公式为

(3)

式中：ρ、v、μ分别为流体的密度、流速和黏性系数；d为特征长度。根据公式可知在低速的空气流动下，腭咽处为层流流动。基于以上判断，本研究选择的计算模型为不可压缩的黏性湍流Realizablek-ε瞬态模型。

2)流固耦合仿真模型

通过流固耦合仿真可以直观地得到在不同条件下的软腭状态，由于软腭的运动会对流场产生不可忽视的影响，所以使用双向流固耦合对软腭的变形做分析。在流固耦合计算中，流体域依然按照之前的流体仿真模型设置，将软腭设置为线弹性可变形的材料，软腭与流体域的接触面设置为流固耦合面。流固耦合区域需要满足以下方程：

(4)

以上方程表示流体域与固体域的应力、位移的大小相等或者守恒。

3)边界条件设置

使用ANSYS workbench中的双向流固耦合模块计算。在仿真计算中，忽略在实际的腭咽闭合过程中气体成分与空气之间的差异，并且完全不考虑肌肉作用对软腭变形的影响，只计算软腭在气流流场作用下的变形情况。

在流体仿真模型中，入口边界为速度入口，出口设置为压力出口，其余壁面设置为固定不可变形壁面。通过调整入口的气流速度大小计算在不同气流作用下的上呼吸道压力分布，分别记录下从0～2 m/s均匀分布的几组数据。为了模拟软腭在呼吸和发音下不同时间的状态，还设置了一组气流的速度变化为按照正弦波形变化。正方向的速度用于模拟呼气过程，负方向的速度则用于模拟吸气过程。

在本研究中，软腭的弹性模量大小被考虑为影响软腭变形的因素，在计算中使用了5 000 Pa、10 000 Pa、15 000 Pa、20 000 Pa和25 000 Pa的5组材料来分析软腭的变形与材料弹性模量之间的关系，泊松比设置为0.45，材料密度设置为1 050 kg/m3。

2 仿真结果与讨论

2.1 上呼吸道内压力分布仿真结果

通过计算可以得到在不同速度入口条件下的上呼吸道内压力分布，在呼气状态下基于该模型的压力分布状态如图2(a)所示，可知在软腭最靠近咽壁处有最大负压，在口腔的入口处有最大的正压。在吸气状态下的压力分布状态如图2(b)所示。可以从气道壁的压力分布中看出，在呼气状态下，口腔处的压力高于咽腔处，口腔的最大压力分布在口腔入口处。咽腔处有较大范围的负压分布，而且最大的负压位置处于软腭与咽壁的最狭窄处。在吸气状态下，咽腔的压力高于口腔。

图2 流体仿真结果

2.2 软腭变形随入口速度和弹性模量的变化

分别计算入口速度为0.5、1.0、1.5和2.0 m/s时的上呼吸道的压力分布，记录其最大正压与最大负压随入口速度变化如图3所示。可以看出上气道中的正压和负压以及压差与入口处的速度呈现某种指数函数的相关关系。根据之前得到的呼气时的压力分布状况，可以得知在软腭的两侧的压差变化随着入口速度的变化在逐渐增大。

图3 上呼吸道压差随速度变化

使用双向流固耦合分别计算从0～2 m/s的几组稳定入口速度下的变形情况，根据之前得到的软腭各个位置处的变形特点，使用软腭处的最大变形来代表软腭的变形程度，得到的结果如图4所示。从图中可以看出在其余条件一致的情况下，软腭处的最大变形与速度也呈指数相关关系。同样可以在图中看出软腭的弹性模量对软腭的变形有较大的影响，在相同情况下，软腭的弹性模量越小，软腭的运动幅度越大。

图4 不同入口速度下软腭的最大位移

为了更好地模拟人真实的呼吸和发音状况，使用正弦速度函数代替之前计算使用的稳定速度入口，入口处的边界速度满足函数v=sin(2πt)，每个周期为1 s，最大峰值速度为1 m/s。对5组不同的弹性模量的软腭材料进行计算，得到在正弦速度作用下软腭最大变形随时间的变化如图5所示。呼气与吸气状况软腭的变形方向相反，在呼气状况下，软腭向咽壁变形；在吸气状况下，软腭向舌根处变形，此处将向咽壁方向的变形定义为正方向，向舌根处的变形定义为负方向。软腭的变形与入口速度相关，软腭在一个周期内的正向最大变形出现在速度的波峰时，软腭在一个周期内反向最大变形出现在速度的波谷时，且在同样的速度下，呼气状态的变形量明显大于吸气状况下的变形。

图5 正弦速度条件下的软腭最大变形

2.3 不同呼吸方式对软腭变形的影响

正常情况下，人主要使用鼻呼吸，口腔气流占总呼吸气流比例大约为4%～8%。但是在某些情况下，人体会适应性地增大口腔气流的占比，口腔气流甚至可达到总气流的70%[7]。不同的呼吸方式对上气道的压降、气流速度等都有很大的影响，也会影响到软腭的变形。在以下研究中，将呼吸方式简化为鼻腔呼气、鼻腔吸气、口鼻腔共同呼气、口鼻腔共同吸气4个过程。

图6(a)为上呼吸道在使用鼻呼气情况下的压力分布，与图2(a)对比可发现压力分布大致相同，最大正压位于口腔处，且最大负压位于咽腔最狭窄处。对比使用口鼻同时呼气和只使用鼻呼气的最大正负压和压差如图7所示。很明显，在只使用鼻呼气的情况下，软腭两侧的压差略大于同时使用口鼻呼气。图6(b)为吸气阶段的气道压力分布，对比图2(b)可以发现其压力分布也大致相同。

图6 鼻呼吸情况下的压力分布

图7 鼻呼气与口鼻同时呼气压力对比

对以上两种呼吸状况分别做流固耦合分析，使用正弦边界条件模拟呼吸过程，得到结果如图8所示。由图可知在使用鼻呼吸的情况下软腭的变形略大于同时使用口鼻。

图8 鼻呼吸和口鼻呼吸模式最大变形对比

3 结语

本文尝试通过流体分析和流固耦合计算找到软腭的受迫运动的影响因素。首先基于CT扫描得到的人体上呼吸道结构图，重建了人体上呼吸道和软腭的模型，并通过流体仿真和流固耦合方法计算了不同的入口速度、软腭弹性强度以及呼吸方式下上呼吸道的压力分布和软腭的变形程度，得到以下结论。

1)不同的呼吸方式对流场的压力分布有较大的影响。在同等速度下，只使用鼻呼气软腭两侧最大压差要比同时使用口鼻呼气大20%～30%。

2)软腭两侧的压力差大小随着气流速度的增大不断增大；且在呼气状态下，口腔侧压力高于咽腔侧，在吸气时，咽腔侧压力高于口腔侧。

3)气流流经口腔和咽腔处产生的压力差是导致软腭变形的一个重要原因。因此，在呼气状态下，软腭向咽后壁移动，而在吸气时向舌后根移动。

4)软腭的变形随其弹性模量的增加而不断减小。

通过仿真分析的方式，分析了多种不同的因素对软腭变形的影响，通过对软腭变形影响因素的不断研究，有望能够更深入地了解腭咽闭合机制并辅助治疗与腭咽闭合相关的疾病。