把“黏液”和“多孔介质”放在一起再用COMSOL去算我脑子里第一反应就不是一个清爽的层流问题而是一坨黏糊糊的东西被挤进海绵里的过程。这个标题里的“蠕动流”其实是一语双关既指低雷诺数下的蠕变流动也指蠕动波驱动流体输运。做生物黏液输运、微流控芯片、药物递送、土壤修复或者油气渗流的朋友大概率会遇到这两件事同时出现流体黏到不行流动又发生在孔隙或者纤维构造里。这篇文章从物理模型怎么选、参数怎么估到COMSOL里怎么搭模型、怎么调求解器把整条路走一遍顺便把我在实际跑模型时踩过的坑都交代出来。只要你手里有COMSOL基础就能照着这条思路开始做自己的算例。1. 动手建模前先想清楚物理场景1.1 蠕动流到底指的是哪一个“蠕动”这个话题最容易翻车的地方就是“蠕动流”三个字在不同行业里指的是完全不同的东西。在流体力学里creeping flow指的是惯性可以忽略的黏性主导流动也就是Stokes流。雷诺数很小的时候NS方程里的非线性惯性项可以扔掉剩下的是压力、黏性应力和外力之间的平衡。很多微流控、渗流问题都属于这一类。另一层意思是peristaltic flow也就是蠕动波驱动的流动。食物在食道里往下走、输尿管把尿液往下送靠的都是管壁上一圈一圈往前进的收缩波。这类问题必须考虑移动边界管壁在振动流体被“挤”着往前走。标题里“黏液遇见多孔介质”这个画面通常意味着两种含义同时成立流体本身很黏雷诺数低所以宏观上有蠕变特征同时又被蠕动波推着穿过多孔结构所以边界是动的、介质是不均匀的。我下面的做法就是把这两层含义都纳入模型。1.2 黏液为什么不能按水来算很多人第一次建模习惯性把流体物性填成水的数据结果算出来完全不对。水的黏度大概是 (10^{-3}) Pa·s而所谓“黏液”是一个很宽的范围胃黏液、宫颈黏液、藻酸盐溶液黏度从零点几到几百 Pa·s 都有。同样是重力的作用下水的唇形下落和蜂蜜的下落完全不同原因就在黏度差了三个数量级。我通常会给一个物性对比表帮自己做判断流体类型动态黏度Pa·s典型特征水约 0.001惯性主导湍流容易发生轻油约 0.010.1层流常见惯性可忽略时接近蠕动流蜂蜜/糖浆约 110黏性主导停止推动后很快“冻结”生物黏液约 0.11000高度剪切变稀静止时很稠黏度高带来的直接后果是速度场对压力的响应几乎瞬时能量主要耗散在剪切流动上惯性项可以安全忽略。这意味着从数值角度这类问题反而没那么容易发散因为流动方程的扩散性很强不容易产生高频振荡。但这不意味着好算——如果黏液还是剪切变稀的非牛顿流体那就有一层新的麻烦我在第三部分专门讲。1.3 多孔介质区Darcy、Brinkman还是统一方程到了多孔介质方程选择是个关键分支。Darcy定律是最简单的模型速度与压力梯度成正比(u -(\kappa/\mu) \nabla p)。它没有粘性剪切项描述的是平均意义上的渗流适合孔隙很小、流动很慢、颗粒直径远小于宏观特征尺度的情况。Brinkman方程则在Darcy项之外保留了Darcy流动的思想额外增加了一个类似NS方程里黏性扩散的项使得它可以在自由流动区域和多孔介质区域之间平滑过渡。它的形式大概是[ \frac{\rho}{\varepsilon_p}\left(\frac{\partial u}{\partial t} (u \cdot \nabla)\frac{u}{\varepsilon_p}\right) \nabla \cdot \left[-pI \frac{\mu}{\varepsilon_p}(\nabla u (\nabla u)^T)\right]\left(\frac{\mu}{\kappa} \beta_F \rho |u|\right) u ]如果你用的是COMSOL的CFD模块里面有一个专门的“自由和多孔介质流动”接口默认就是Brinkman统一方程把自由流动区和多孔区做成同一个物理场界面不再需要额外指定通量连续性条件。我强烈建议优先用这个接口而不是手动把Darcy和层流拼接起来。后者的界面条件处理起来非常容易出问题尤其是在动边界存在的情况下。2. 关键参数估算决定成败的几张“底牌”2.1 特征尺度和无量纲数先算一遍我见过太多人几何模型画得漂漂亮亮材料参数一拍脑袋填进去算出来云图好看但不知道自己在算什么。建模之前建议先把无量纲数算清楚这是判断“这个模型该不该这么简化”的底牌。以我最近做的一个算例为例一根半径 (R 1) mm的管中间有一段多孔材料长度 (L 20) mm蠕动波波速 (c 5) mm/s流体密度 (1000) kg/m³黏度 (0.1) Pa·s。先算雷诺数[ Re \frac{\rho c R}{\mu} \frac{1000 \times 0.005 \times 0.001}{0.1} 0.05 ][ Re 0.05 \ll 1 ]这已经说明惯性作用可以忽略除非你特别想研究入口处的局部加速效应否则没必要纠结NS方程的非线性项。再看多孔区的渗透率 (\kappa)。用Kozeny-Carman公式估算[ \kappa \frac{\varepsilon_p^3 d_p^2}{180(1-\varepsilon_p)^2} ]取孔隙率 (\varepsilon_p 0.5)颗粒直径 (d_p 0.1) mm则[ \kappa \frac{0.125 \times (10^{-4})^2}{180 \times 0.25} \approx 2.78 \times 10^{-11} \ \text{m}^2 ]然后算Darcy数[ Da \frac{\kappa}{R^2} \frac{2.78 \times 10^{-11}}{10^{-6}} \approx 2.78 \times 10^{-5} ]这个数很小说明什么问题说明多孔介质本身的渗透阻力远远大于自由流动区的黏性阻力。流体会倾向于绕过或者只在阻力最小的通道里流过压力大部分都消耗在多孔段上。这个判断直接影响你后处理应该关注什么不是管壁附近的剪切而是多孔段前后的压降、以及进入多孔段的有效流量。2.2 蠕动波参数怎么写进方程蠕动波最常见的是行波形式[ h(x,t) R A \sin\left(\frac{2\pi}{\lambda}(x - c t)\right) ]其中 (A) 是壁面径向振幅(\lambda) 是波长(c) 是波速。这里需要区分两个概念如果壁面是刚性波导那么只有这个位移是真实的运动边界如果整个管子在蠕动那还要考虑壁面厚度和材料的弹性。COMSOL里处理行波壁面通常把A设成管径的5%15%太大会导致动网格畸形太小又看不出泵送效果。波数 (k 2\pi/\lambda) 和波速 (c) 共同决定壁面行波的相位速度。在多孔耦合问题里我们关心的是“壁面的波动是否有效地把流体推过多孔区”。如果波速远大于流体实际渗流速度壁面波动的大部分能量会耗散在自由流区多孔区几乎没有反应如果波速和渗流速度匹配合适就会形成明显的泵送效应。我建议把这几个参数全部设置成全局参数后面做参数化扫描会非常方便R 1e-3 [m] 通道半径 L 20e-3 [m] 总长度 A 1e-4 [m] 壁面振幅 lambda_w 10e-3 [m] 蠕动波波长 c_w 5e-3 [m/s] 蠕动波波速 rho 1000 [kg/m^3] mu0 0.1 [Pa*s]2.3 阶段化简先做稳态还是直接上瞬态还有一张底牌是“要不要一上来就做瞬态”。蠕动壁面是运动的本质上必须瞬态求解但你可以分两个阶段走第一阶段关掉动网格把多孔区两端设成固定压力边界算一个纯压力驱动的稳态渗流场作为后续瞬态计算的初始值第二阶段再开启变形几何把壁面位移加上从稳态解继续瞬态推进。不要一上来就从静止状态直接开算蠕动波否则前几个周期都是在冲刷初始暂态白白浪费计算量还容易在最初几个时间步引起数值振荡。我在实际项目里基本都是这个流程实测下来比直接全耦合要稳得多。3. COMSOL里一步步搭起耦合模型3.1 物理场接口选择优先“自由和多孔介质流动”进入COMSOL 6.x新建模型时选择CFD模块下的“自由和多孔介质流动fp”接口。这个接口的好处是把自由流动区用层流/NS方程描述多孔区用Brinkman方程描述两者在一个接口里自动衔接界面上不再需要手动加“应力连续”或“通量连续”的条件省掉很多麻烦。如果你用的是基础版没有CFD模块也可以用“层流”加“Darcy定律”两个接口手动耦合但在界面处必须定义好压力连续、法向速度连续动网格情况下还要处理网格运动与Darcy域变形的协同。我的经验是除非你特别清楚自己在干什么否则别走手动耦合这条路。几何模型的搭建也比较直白二维轴对称是首选计算量小又能反映管道径向分布如果你的通道不是严格的圆管再考虑二维平面或三维。整个模型可以分成三段入口自由流区保持可变形状态壁面施加蠕动波位移中间多孔介质区孔隙率 (\varepsilon_p 0.5)渗透率 (\kappa 2.78 \times 10^{-11}) m²这个区域通常不参与网格变形出口自由流区同样可变形保留蠕动波的延续。中间多孔区和前后自由流区之间是两个交界面交界面处网格要加密但不要让多孔区也跟着壁面一起“蠕动”。我通常的做法是把变形区域限制在自由流区让多孔区完全固定只在交界面上把位移设成0保证不撕裂网格。3.2 材料定义Newton到Carreau的过渡如果黏液是牛顿流体材料设置很简单填密度和黏度就行。但实际生物黏液多半是剪切变稀的COMSOL里可以用“Carreau模型”定义剪切依赖黏度[ \mu_{eff} \mu_\infty (\mu_0 - \mu_\infty)\left(1 (\lambda \dot{\gamma})^2\right)^{(n-1)/2} ]其中 (\mu_0) 是零剪切黏度(\mu_\infty) 是极限剪切黏度(\lambda) 是松弛时间(n) 是幂律指数(n 1) 表示剪切变稀。我常用的一组参数是参数值说明(\mu_0)10 Pa·s静止时很稠(\mu_\infty)0.01 Pa·s高速剪切时接近水(\lambda)1 s松弛时间(n)0.4剪切变稀程度这组参数的好处是能明显观察到剪切变稀对渗流的影响。壁面蠕动振幅增大剪切率升高表观黏度下降流量会有非线性增加这个现象在牛顿流体模型里永远看不到。许多人不理解为什么数值收敛差其实就是因为在某个剪切率区间黏度曲线斜率太陡让雅可比矩阵估计失真。因此如果你只是在试着验证模型建议先把流体设成牛顿流体把整个流程走通再回到Carreau模型。这是很经典的经验——一步一步加复杂度。3.3 蠕动壁面与变形几何的设置这一步是本算例的“灵魂”。在COMSOL的“变形几何”接口下对入口/出口自由流区的壁面施加“指定位移”。假设管道轴向是 (x)径向是 (r)那么径向位移写为[ d_r A \sin\left(\frac{2\pi}{\lambda_w}(x - c_w t)\right) ]轴向位移可以设成0也可以加一个很小的分量以模拟蠕动的轴向拉伸但我不建议一开始就加轴向位移会让网格变形形态变得不确定。先把纯径向行波跑通再考虑轴向耦合。需要注意一点COMSOL里的坐标变量在边界上使用 (x) 和 (r)但如果你设定的是空间位移它会反映成网格运动。同时“壁面”边界条件要确保流体的速度跟随壁面即设置“壁面运动”为“移动网格速度”否则壁面的速度是不会传递给流体的。多孔介质区域建议设置为“固定网格”不参与变形几何计算。在自由流区与多孔区的交界面上指定位移为0。如果要减少交界面的网格畸变可以在交界面两侧预留一小段缓冲自由流区让壁面的位移从完好差值逐渐衰减到0而不是直接生硬切到多孔区边界。3.4 网格剖分别让动网格毁在第一步动网格问题里网格质量是最容易出问题的特别是在壁面振幅较大的情况下。我的网格策略分三块自由流区近壁面边界层至少要23层边界层网格但高黏度流体湍流少见边界层不像空气动力学那么苛刻。网格厚度可以适当放宽因为剪切层比较厚而不是贴在壁面上。沿轴向的网格密度要和波长匹配一个波长内至少保证1020个单元。如果单元太少行波壁面会呈现出多边形折角造成人为的波形失真影响结果可信度。多孔区不需要太密但孔隙尺度肯定没法直接解析这里用的是宏观Brinkman方程所以局部细化到能平滑描述压力梯度即可。另外在变形几何里可以用“自动平滑”的网格位移方式但如果振幅相对半径超过10%建议改成“超弹性平滑”它可以更好地处理大变形代价是每个时间步需要额外求解一次网格更新方程计算量会上升。选择一个合适网格策略后最好先做一次单周期试算看看最大网格偏斜度变化趋势如果偏斜度持续增大后面迟早会负体积不如趁早调整振幅或平滑方式。3.5 求解器设置与瞬态策略物理场和网格都设好后打开“研究”用“瞬态”研究求解器我习惯这么配置时间格式BDF向后差分最大阶数设为2。高阶格式对非线性、移动网格问题并不总更稳2阶足够。初始步长设成周期长度的1/100比如T_wave lambda_w / c_w # 行波通过一个波长所需时间 dt_init T_wave / 100后续最大步长控制在1/20周期以内保证瞬态精度。一致性初始化打开避免初始时刻因边界条件突变产生伪振荡。如果出现非线性迭代发散先把“恒定阻尼牛顿”的阻尼因子降到0.7左右再把时间步长缩小。暂时不要调大网格密度发散首先应该怀疑是时间步长或阻尼问题而不是网格。这是一个非常有效的排查顺序能节省大量等待重算的时间。4. 踩坑记录收敛失败与网格畸变的排查思路4.1 最常见症状非线性求解器不收敛先给一个速查表对照自己的情况去定位症状可能原因对策时间步大量缩短、反复重试初始步长过大初始步长减到周期的1/200流体压力出现高频振荡惯性项打开造成小扰动检查Re是否很小可尝试忽略惯性项多孔区压力异常高渗透率参数过小先用较大渗透率试算再逐步减小网格畸变、负体积振幅/波长比过大减小A或改用超弹性平滑增加缓冲区非牛顿黏度造成振荡Carreau模型斜率过陡先用牛顿流体建模再引入非牛顿项这些坑我几乎都踩过。最典型的是第一次跑初始步长设置成整个周期的1/20结果前三个时间步就疯狂缩短最后卡在 (10^{-8}) 量级的步长上寸步难行。后来才发现是初始场对不上先跑稳态流场作为初值再把初始步长调小问题立刻消失。4.2 动网格负体积别急着改网格做蠕动壁面时最容易出现的就是网格在某一时刻负体积求解器报错退出。新手第一反应通常是加密网格或换网格类型但我建议先检查振幅和波长比。如果 (A/\lambda) 大于某个临界值任何网格细化都救不了因为物理变形本身就太大了网格会在过度挤压的地方发生翻转。解决思路有两条减小振幅。这是最直接的先证明模型逻辑正确再提高振幅。比如把振幅从管径的15%降到5%。在变形区与固定区之间增加缓冲层。多孔介质区不参与变形交界面位移设为0如果一个网格单元一半要动、一半要不动它迟早会翻转。缓冲层让变形逐渐过渡网格形态平滑很多。4.3 多孔区压力虚高先检查渗透率量级有时候模型能跑通但后处理发现多孔区里压力超出了物理预期比如压降达到几百千帕而流量仍然很小。这通常不是求解器问题而是渗透率参数没校准。回到我前面算的例子渗透率 (2.78 \times 10^{-11}) m²对应的是0.1mm粒径、50%孔隙率的多孔介质。如果你随手填一个 (10^{-14}) m²的数据那就相当于在堵死的水泥块里渗流压力当然会异常高。判断渗透率是否合理可以用Kozeny-Carman反推一下等效粒径看看是否符合你的物理场景。4.4 界面处速度不连续如果用的是我自己不推荐的手动耦合“层流Darcy”可能会发现自由流区和多孔区交界面两侧的速度有突变。原因在于Darcy定律本身不求解速度场的剪切结构它给的是达西速度和自由流区NS方程的表面速度在定义上就不同。如果你一定要这样耦合就需要在界面处正确实施法向应力连续和速度连续条件。用“自由和多孔介质流动”接口就不会有这个烦恼因为多孔区的Brinkman方程保留了剪切项界面处理是自动完成的。反正我实测下来这个接口省力不是一点半点。5. 后处理技巧从云图里挖出有效渗透率5.1 速度场和压力场到底该看什么跑完瞬态周期解后后处理不是看看动画就完事。首先要看的是“一个完整蠕动周期内的平均流量”。因为蠕动波是周期性的瞬时流量在正负之间振荡真正有意义的是净流量——也就是流体到底被输运了多少。在COMSOL里可以定义“派生值”里的“线积分”或“表面积分”对截面上的法向速度分量做积分然后随时间做周期平均。这一步如果只取瞬时值会得到非常误导的结果比如某相位的瞬时流量为零但净流量其实很大。压力方面重点提取多孔介质段前后两个端面的平均压力差。这个压差加上净流量就可以通过达西公式反推有效渗透率[ \kappa_{eff} \frac{u_D \mu L}{\Delta p} ]其中 (u_D) 是达西速度也就是净体积流量除以总截面积。这个“有效渗透率”是在蠕动驱动而非恒定压差条件下测到的平均值它会在不同波速、不同振幅下变化。把这个参数作为结果输出比单纯放几张云图有说服力得多。5.2 参数化扫描看振幅、波速、黏度怎么影响输运模型能稳定跑完一个周期后最强的工具就是参数化扫描。我建议至少做三组扫描振幅从 (0.05R) 到 (0.15R)看净流量是否随振幅近似线性增加还是有明显的阈值效应波速从 (0.001) m/s 到 (0.02) m/s看是否存在一个最优波速让通过多孔区的净流量最大黏度从牛顿流体的 (0.01) 到 (1) Pa·s看压降和流量之间的关系是否符合Darcy的比例关系。扫描完把净流量画成曲线通常会有一些有趣的转折点。比如波速过高时壁面动得太快多孔区来不及响应净流量反而下降波速过低时输运太慢又不现实。这个最优波速如果存在就是你的设计参数。5.3 关于“效率”如何量化蠕动泵送的有效性吐一个我自己比较关注的概念泵送效率。蠕动波壁面对流体做功能量一部分变成流体动能大部分被黏性耗散还有一部分是克服多孔介质阻力消耗掉的。在COMSOL里可以通过后处理计算壁面的“法向应力”和“移动速度”点积的积分得到壁面对流体输入的总功再结合净流量和出口压力计算流体获得的压力能与流量乘积两者相除得到一个无量纲效率。这个值在优化蠕动泵、生物输运过程时非常有用也是很多论文的核心输出。坦白说我刚做这个模型时也没有一上来就算效率后来发现审稿人都爱问“你的蠕动参数为什么取这么大依据是什么”有了效率曲线所有参数选择都有据可依。最后分享一点实际操作中的体会如果你要复现这个算例我的建议很简单先别急着把非牛顿、多孔流动、动网格全部堆上去。第一步用牛顿流体渗透率给大一点壁面振幅给到5%管径把整个周期跑通确认网格没有负体积、流量曲线有周期规律第二步再加多孔介质把渗透率降到目标值观察压降变化第三步加剪切变稀的非牛顿模型最后再把振幅和波速提到目标设计范围。每一步都验证一个物理现象出了任何问题都能快速定位。COMSOL这类多物理场软件里九成的不收敛其实都是“复杂度上太快”导致的而不是软件本身做不到。慢慢加、逐层验证是我做过这么多耦合仿真后最大的心得。