《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Efficient fluid structure interaction simulation of vocal fold oscillations using a homogenized Lattice Boltzmann Method
编辑推荐：
本文针对人体发声过程中声带振动的流固相互作用(FSI)模拟所面临的高计算成本和传统方法难以处理声带接触的难题，研究人员创新性地提出了一种将六质点模型(6MM)结构动力学映射为时变孔隙率场，并与基于Brinkman-Navier-Stokes方程过滤模型的均匀化格子玻尔兹曼方法(HLBM)进行双向耦合的新策略。该研究在三维喉道通道中实现了稳定的自持振荡，获得了符合生理学特征的结果（基频f0≈248 Hz，平均声门下压≈1.01 kPa），并成功处理了每周期完全的声门闭合。这项研究极大地提高了计算效率（在NVIDIA RTX 2000 Ada笔记本电脑GPU上可达60个振荡周期/小时），为开展以往传统FSI方法难以实现的系统性发声参数研究提供了强大工具。
人类的声音是如何产生的？这背后是一场发生在喉部、涉及气流、软组织与声学相互作用的精妙“舞蹈”。发声的核心在于呼气气流诱导声带的周期性振动，即肌弹性-空气动力学原理。这个过程本质上是一个复杂的流固-声相互作用问题。然而，对这一过程进行深入研究却面临巨大挑战。在体研究极为困难，因为清醒的人体难以在喉部放置传感器，即便是可视化观察也容易引发不适。因此，计算建模成为理解这一生理过程的关键工具。但计算模拟人类发声同样荆棘密布：它需要处理从层流到湍流的过渡流态、声带组织的大变形、以及每个振动周期中必然发生的声带间自接触（碰撞）问题。传统的模拟方法，如任意拉格朗日-欧拉法，在处理接触时面临网格畸变和坍塌的难题；而浸没边界法则需要复杂的界面追踪，计算开销巨大。这些限制使得系统性的参数研究（如比较不同性别、年龄或发声模式的差异）在计算成本上几乎不可行。为了在保证生理真实性的同时实现高计算效率，一项发表在《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》上的研究提出了一种全新的流固耦合模拟框架。
该研究主要应用了几个关键技术方法：1. 采用基于过滤的Brinkman–Navier-Stokes方程的均匀化格子玻尔兹曼方法进行流体动力学求解，该方法能同时处理湍流亚网格尺度模型和作为移动固体边界惩罚项的孔隙介质模型。2. 使用经典的六质点模型来描述声带组织的结构动力学，该模型将每侧声带简化为三个质点，通过弹簧和阻尼器连接，能够捕捉其弹性运动和碰撞。3. 提出了一种创新的耦合策略，即将6MM计算得到的声带位移和速度，映射为HLBM求解域中的一个时变孔隙率场，从而隐式地处理复杂的几何变化和接触，无需重新划分网格。4. 整个计算流程在开源LBM框架OpenLB中实现，并利用现代通用图形处理器进行加速，显著提升了计算效率。
2. 方法
2.1. 格子玻尔兹曼方法
研究人员使用均匀化格子玻尔兹曼方法对描述流体运动的过滤Brinkman-Navier-Stokes方程进行离散求解。该方法通过一个时变的孔隙率场参数d(x,t)来表征流固界面。在完全流体区域d=1，在固体内部d=0，在界面过渡层d∈(0,1)。固体边界的速度通过一个惩罚项（与渗透系数K相关）被施加到流体上，从而实现无滑移边界条件。同时，采用Smagorinsky大涡模拟模型来捕捉过渡湍流。
2.1.3. 变形声带的表示
这是本方法的核心创新。研究人员预先对一个参数化的“主”声带几何进行体素化，生成一个包含从最大会聚到最大发散所有形态的参考网格库。在模拟的每个时间步，根据6MM计算出的声带各段位移，从库中选择对应的几何切片，并将其转换为当前流体网格上的孔隙率场。这种方法避免了每个时间步都进行复杂的几何重构和距离场计算，效率极高。
2.2. 结构模型
结构侧采用经典的六质点模型。每侧声带被建模为三个质点（分别代表下、中、上部），质点间通过线性弹簧和阻尼器连接，并附着在固定的喉部结构上。模型包含了横向和纵向的刚度，能够模拟声带的整体运动和黏膜波。当两侧声带的质点发生重叠时，模型会引入非线性的接触力。
2.3. 模拟设置
研究在一个简化的三维喉道通道中进行模拟。通道包含声门下区域、声门（声带所在处）和声门上区域。在入口施加恒定压力，出口为开放边界。流体网格分辨率为△x=50微米。初始时，两侧声带处于闭合接触状态，然后施加压力差以启动振荡。
2.4. 数据分析
通过监测声门处的横截面积、体积流量、上下游压力差等参数，来分析振荡的稳定性、基频、开放商等生理学指标。同时，也分析了声门上射流的特征和声带接触力。
3. 结果与讨论
3.1. 耦合时间步长的影响
研究首先评估了流体求解器与6MM结构求解器之间的耦合时间步长△tc对结果的影响。发现当△tc小到足以解析碰撞过程时（本研究中≤2微秒），模拟结果稳定且对步长不敏感。最终选择△tc=2微秒以确保精度和效率的平衡。
3.2. 一般特性
模拟成功激发了稳定、自持的声带振动。振荡呈现出典型的相位差：下部先于中部打开，中部先于上部打开，关闭顺序则相反，这与生理观察到的黏膜波现象一致。每个周期都实现了完全的声门闭合。
3.3. 声门参数
从稳定的振荡周期中提取的关键生理参数均落在合理范围内：基频f0≈248 Hz（接近成人语音范围），平均声门下压力≈1.01 kPa，开放商（声门开放时间占周期比例）≈0.64。接触力峰值约为27 mN。
3.4. 声门上射流特征
在声门上方观察到了预期中的脉动射流。射流在脱离声门边缘后产生剪切层不稳定性并发展成涡结构，这些涡的生成和耗散是产生声音源的主要机制。
3.5. 质点模型参数研究
通过系统改变6MM中的质量、刚度、阻尼等参数，研究展示了该耦合模型进行参数化研究的强大能力。例如，增加质点质量会导致基频下降，而增加刚度则使基频升高，这与物理直觉和简化理论模型一致。这证明了该框架非常适合用于探索不同生理或病理条件下的发声机制。
结论与意义
本研究成功地开发并验证了一种用于模拟声带振动流固相互作用的高效计算框架。其核心贡献在于将均匀化格子玻尔兹曼方法与六质点模型通过一个时变孔隙率场进行耦合。这种方法巧妙地绕过了传统方法在处理移动边界和接触时的网格难题，实现了对声带完全闭合的稳健模拟。
该研究的首要重要意义在于其卓越的计算效率。在单个消费级笔记本电脑GPU上，模拟速度可达每小时60个振荡周期，这比许多依赖大型高性能计算集群的传统CFD-FSI方法快了几个数量级。这种效率的飞跃使得大规模、系统性的参数研究变得可行。研究人员可以轻松地探索声带组织属性（质量、刚度）、几何形状、气流条件等大量参数的变化如何影响发声，这对于理解不同性别、年龄、发声模式（如歌唱）的差异，乃至声带病变（如息肉、水肿）的病理机制具有巨大潜力。
其次，该方法在保证生理合理性的同时提供了足够的工程简化。6MM虽然是对复杂多层声带组织的极大简化，但其与高保真流体求解器的耦合仍能再现关键的生理观测指标（如基频、压力、黏膜波、声门闭合）。这种“足够好”的模型与高效率的结合，为开发可用于临床评估或嗓音训练的“数字声带”工具铺平了道路。
总之，这项研究为发声生物力学领域提供了一个强大的新型计算工具。它将模拟的重点从耗时数周才能完成的孤立案例研究，转向了可以快速进行的参数扫描和假设检验，有望加速我们对人类嗓音这一复杂而美妙的生理现象的理解。
