前阵子整理一个超声微流控项目的结题材料我自己重新读到了当初写在报告里的那句话本模型采用声固耦合和两相流耦合多物理场使用的模块包括声流层流、相场、压力声学、固体力学模块。那句话写出来只有几十个字但真正动手做过的人会知道每一种物理场单独拎出来都能折腾一个人好几天而把这几个词放到同一个模型里时最先要回答的是四个物理场到底怎么分工、谁先算、谁给谁提供源项、怎样收敛以及最后算出来的界面运动能不能和实验里液滴移动的现象对得上。这篇文章想聊的就是这类“声固两相流”组合背后的完整建模链路以及我在实际项目里反复试错后沉淀下来的做法。如果你正准备做超声驱动液滴、气泡或微通道内两相界面操控的仿真又恰好被这一长串模块名搞得不知道从哪下手那这套思路应该能帮你省掉不少弯路。1. 谁在振动、谁在流动、界面交给谁模块职责拆解1.1 压力声学加固体力学它们是一对振动搭档很多人看到“声固耦合”四个字第一反应是声波打在固体上然后反弹。这个理解没有错但在这个项目里更准确的说法是固体结构负责把激励源的振动送进液体液体里的声压又会反过来给固体施加动态载荷。也就是说这不是两条独立计算的路线而是一对互为边界条件的搭档。压力声学模块通常用于流体或气体区域求解声压场。它关心的物理量是声压 (p) 和声质点速度支撑的是波传播、驻波形成、辐射与散射这些现象。固体力学模块则用来描述基底、通道壁或其他弹性结构输出的是位移场、应力场和应变场。两者的分界面在哪里声固耦合就发生在哪里一侧是液体声压给固体表面施加压力载荷另一侧是固体表面的法向振动速度把能量回灌进声场。在这类微流控芯片里常见的结构是压电换能器贴在玻璃或硅基底上基底上加工了几十到几百微米深的微通道通道里充满液体和第二相液滴。压电片通电后会产生高频振动振动从固体基底传到液体里在通道中形成驻波场。如果项目里不打算引入压电模块工程上常用一个等效做法在压电片接触位置施加指定频率的边界位移或法向加速度。只要关心的区域离激励源有足够距离这样简化不会明显影响液体里的声场形态却能省掉一套压电本构方程收敛难度会下降不少。1.2 层流加相场它们负责液体怎么动、界面长什么样层流模块解决的是微尺度下液体的流动问题。微通道里的雷诺数通常很低流动处于典型的层流或蠕流状态不太涉及湍流模型。但这套模块名字里有个“声流”前缀它并不是指一个独立的特殊物理场而是指在层流求解中加入由声场产生的时均驱动力从而使液体产生一种缓慢而确定的流动也就是常说的声流。相场模块则是用来追踪第二相界面的工具。它不直接告诉你“这个液滴的边界在哪个坐标点”而是通过一个连续变化的相场变量 (\phi) 描述两种液体在空间中的分布。(\phi) 在主体相里接近一个值在另一相里接近另一个值中间通过一条有厚度的光滑过渡层衔接。这样处理的好处是液滴的合并、分裂、大变形不需要人为维护移动网格拓扑变化会自动发生。缺点是引入了额外的数值参数后面我会专门讲这些参数怎么调。1.3 四套物理场握手实际上只有四个耦合通道拆开单个物理场并不难难点在“握手”的方式。以微通道中一个油滴被声波驱动时的情形为例整套模型的耦合关系可以归纳成下面几条耦合通道从哪个物理场出发作用到哪个物理场传递的物理量声固边界耦合压力声学固体力学声压载荷声固边界耦合固体力学压力声学结构振动法向加速度声致流场源项压力声学层流时均雷诺应力/体积力两相界面输运层流相场速度场输运界面界面反作用相场层流表面张力、局部密度与黏度前两条属于声固耦合范畴处理的是声波如何进入液体第三条是“声流层流”的核心把声场能量转化为宏观流动第四、五条则是标准的层流两相流耦合。标题里说的“两相流耦合多物理场”最容易被忽视的其实就是第五行相场变量影响局部密度和黏度而密度和黏度又直接进入层流方程这两套方程必须联立求解不是先用相场算个界面再把界面作为固定边界去算流场就行的。2. 频率差六个数量级求解顺序必须提前设计2.1 先判断时间尺度再决定要不要“全耦合瞬态”建模之前最不该省的一步是估算不同过程的时间尺度。工程超声频率常常是数百千赫兹到数兆赫兹取一个典型值1 MHz那么声波的一个周期就是1微秒。而微通道里的流动速度通常在毫米每秒量级一个500微米长的通道液体流过需要大约0.5秒。如果想从声波发射一直模拟到液滴被推动几倍直径的距离仿真时长至少要覆盖几十毫秒到几百毫秒这是声周期数量的几万到几十万倍。如果直接在瞬态求解器里同时捕获声波传播和液滴位置变化时间步长必须满足声学周期性否则连声场都会失真。而界面移动需要的时间远远大于这个尺度这就等于要求计算机在每一个微秒级别的时间步上都去求解一套包含四种物理场的大型方程组然后还要重复几十万步。我可以负责任地说常规工作站上基本跑不动而且即使跑得动大量计算量都消耗在重复计算高度振荡的声波上对界面运动几乎没有额外信息贡献。所以这类模型的标准解法是“频率分解”把声场当作一个快速周期过程在频域内计算把声流和两相界面运动当作一个慢速过程在时域内求解。声波对上百万次的振荡最后给液体产生的净效果是一个时间平均值这个平均值就是频域声场可以直接提供的。于是求解链可以从一个看似复杂的瞬态耦合问题降级成“频域声固耦合 带源项的层流两相流”两段任务。2.2 界面要不要跟随声场更新先学会“冻结”有了两步法还不够还有个问题绕不过去两相界面位置不是固定的而声场本身非常依赖边界。液滴从位置A移动到位置B通道里声压分布也会变。严格来说声场和界面运动存在双向耦合。但如果每一步都让声场跟随着界面的微小移动重新求解计算开销仍然很大。我实际采用的办法是分阶段处理。模型搭建初期先假设界面在当前时刻是固定的也就是“冻结界面”假设用这个固定界面算出一个声场再把声场源项放进两相流模型让液滴移动一段距离。等这一段运动结束后再根据新的界面位置重新计算声场如此循环。这个流程听上去是粗糙近似但在很多场景下精度足够。原因是油和水之间的声学特性差异通常不如气液之间那么悬殊界面对声场的影响相对有限界面移动对整体声场结构的改变往往滞后于液滴位置变化。如果研究对象换成气泡这个假设就要小心了气泡在水里的声阻抗差异太大声场遇到气泡会发生强烈的散射和局部驻波那时就必须把声场更新插到界面演化进程里。2.3 声流来源不同模型厚薄也不一样“声流”这个笼统的词实际包含两个物理来源建模方法差别很大。一种是体声流也叫Eckart流由声波在液体介质内部被吸收、能量衰减产生驱动力量级取决于液体的吸收系数。另一种是边界声流也叫Rayleigh流由声波与黏性边界层相互作用产生液体靠近固体壁面时通过边界层的非线性效应形成稳态流动。对于兆赫兹级超声水中的声吸收系数通常很小体声流在宏观尺度下往往不是主要贡献但在微通道这个尺度边界声流会变得非常显著因为通道尺寸和声粘性边界层厚度处在一个能发生强烈相互作用的范围。这里的理论计算很简单边界层厚度 (\delta_v\sqrt{2\nu/\omega})。水在1MHz下运动黏度约 (10^{-6}\ \mathrm{m^2/s})算出来边界层厚度只有大约0.56微米。通道壁附近会出现强剪切层而液滴和颗粒恰恰最容易在这一层里受到横向作用。因此第一种声流可以在层流方程里加一个分布式的体积力源项第二种声流则往往需要在壁面处施加等效滑移速度或用特殊边界层网格来模拟。先明确项目里到底是哪一种声流起主导作用再决定网格和边界条件的复杂度不然很可能花了大力气去画0.5微米的网格结果界面动力学的核心却被忽略了。3. 实操三步先算声固耦合再算声流最后接入相场界面3.1 第一步在频域里把声固耦合调稳我在COMSOL这类商业软件里的习惯是先不要一上来就建完整两相流模型而是先做一个“压力声学频域 固体力学 声结构边界”的模型把单纯的声音传播和结构振动验证清楚。几何可以先用二维表示取微通道的纵向截平面包含上方的基底固体、中间液体通道、下方的换能器激励区域。二维模型跑得快参数调起来很顺手等物理机制确认了再考虑是否升级到三维。材料参数是最先要填的。液体区域要给定密度和声速比如水取998 kg/m³、声速约1480 m/s固体区域如果是玻璃基底要给定密度、杨氏模量和泊松比有些材质还应当加一点结构损耗因子防止频域求解时共振峰过于尖锐导致压力虚高。边界条件方面激励位置可以设置成简谐位移边界或法向加速度边界。频率不要随便取最好先用特征频率研究算一下看通道内能否形成驻波模态。如果激励频率落在驻波模态附近声压幅度才会稳定放大。求解后不能只看云图颜色好看重点检查三点。第一声压极大值位置是否和理论驻波节点一致第二结构位移幅值是否处于合理量级比如几纳米到几十纳米第三声场在边界处有没有异常反射造成的压力堆积。如果这三个检查都过了再往模型里加后面的物理场。3.2 第二步把声场能量转成体积力单独跑一次层流声流从频域声场到层流体积力的这一步是整套建模中最容易出错的地方。声流体积力的本质是声波非线性项的时间平均值常用雷诺应力散度来描述。如果用公式写大致是[ F_i-\frac{\partial}{\partial x_j}\left(\rho_0\langle v_i v_j\rangle\right) ]其中 (v_i)、(v_j) 是声质点速度分量尖括号表示对一个或多个声周期取时间平均。这里最容易犯的错是直接用瞬时声压梯度当作力源那样计算出来的流场会随声波周期来回抖动根本不是一个稳态的声流。在实际软件里实现这个源项有两条路。如果你用的是可以自定义方程的软件可以直接定义声质点速度变量然后写出雷诺应力散度表达式加到层流方程的体积力里。如果嫌麻烦也有简化方案对一维驻波场或边界声流可以先做局部区域的平均得到等效的体积力或边界滑移速度再输入到层流模块。有一个我强烈推荐的验证动作先不要加相场只在纯单相液体里把层流声流跑一遍。你会看到声场驻波节点两侧形成一对对旋转方向相反的流涡这就是典型的声流结构。如果这一步的涡心位置和大小与文献观测对不上说明源项或者边界层处理有问题此时排查成本很低一旦直接跳到两相流模型里再发现问题就很难分辨是声流的锅还是界面参数的锅。3.3 第三步相场两相流开始接管液滴运动声流验证通过后把物理场升级为“层流两相流、相场”接口同时让声场的频域结果保持可用状态。这个阶段要设置两类参数一类是材料本身的物理性质比如两种液体的密度、黏度和界面张力另一类是相场模型的数值参数包括界面厚度参数和迁移率。相场方法的基本控制方程属于Cahn-Hilliard型形式上可以写成[ \frac{\partial \phi}{\partial t} \mathbf{u}\cdot\nabla \phi\nabla\cdot(M\nabla \eta) ]其中 (\eta) 是化学势与界面厚度参数和混合能密度有关。这个方程中的 (M) 就是迁移率控制相场变量扩散的快慢。调参时一个常见误区是把 (M) 调得很大希望界面快速达到平衡结果反而导致质量不守恒或出现细小碎滴。初始液滴位置用一种“相场初始化”功能设置即可比如一个半径50微米、球心在特定位置的油滴。但要注意初始化后的界面并非立刻处于力学平衡这时如果直接施加强声流源项液滴会在表面张力尚未稳定时就被冲变形导致非物理解。我通常会在前几个毫秒关闭声流源项只让表面张力把液滴慢慢松弛到圆形再开启声场驱动。这一步等于给两相流一个“热身期”能明显降低后续数值发散的概率。3.4 完整更新循环怎么设计工程上如果要让界面运动过程中声场持续更新就需要设计一个循环逻辑。简单方案是把时间推进分成若干段每一段内固定界面位置和声场计算一段液滴运动然后暂停流动求解把更新后的相场变量映射回声学模型重新计算声场再继续推进下一段。这种块状更新效率远高于每个时间步都全耦合而且便于在每段结束后检查结果是否合理。软件层面可以用研究步骤里的“辅助扫描”或者脚本驱动来实现也可以手动分段跑。只要分段长度不要太大保证每段内液滴位移远小于声波波长的尺度结果就是可接受的。4. 网格、界面厚度、声粘性边界层三张网如何兼顾4.1 先算一算尺度错配有多严重这套模型有个天然麻烦不同物理过程对网格尺寸的要求差距极大。按一个典型模型来看设液体中声速约1500 m/s频率1 MHz则声波波长约1.5毫米。即使按一个波长布置10个单元声学网格最大也能到150微米这并不苛刻。而相场界面厚度参数通常取1到3微米为了解析界面网格不能比这个值大太多否则界面会变成一条参差不齐的过渡带。两者对比声学对网格的“下限要求”和相场对网格的“上限要求”相差了上百倍。更极端的是声粘性边界层0.56微米。如果真的要在壁面处把这个边界