用Comsol仿真纸片吸水-毛细现象

纸条上含水饱和度随时间的变化.

这个案例来自Comsol案例库。中文标题为《纸条芯吸》。

纸条是一种多孔材料。当干燥的纸片与液体接触时，由于毛细力的作用，液体会逐渐被纸条吸附向上，直到重力与毛细力达到平衡。

Comsol模型与物理场

模型是一个很简单的二维矩形。平面问题，纸片厚度定义为1mm。

这个分析涉及的物理场比较复杂。所以相对的，几何模型就要简单一些。一方面便于收敛，另一方面重点展示使用的物理场和分析技术。

模型中包含两个物理场，分别是【多孔介质相传递】和【达西定律】。预置的多物理场耦合节点为【多孔介质多相流】。 

模型树

这里，我们就依次了解一下这两个物理场描述的是什么。

多孔介质相传递 物理场

首先是【多孔介质相传递】。这个物理场的文档目录位于 CFD Module > User's Guide > Multiphase Flow Interfaces > The Phase Transport Interfaces > The Phase Transport in Porous Media Interface. 可见它属于流体CFD模块的多相流模型中。

相传递物理场，用于模拟多个不相溶的相 在自由流动或多孔介质中的输运现象。这个物理场求解的场变量是每个相的平均体积分数（averaged volume fractions），在多孔介质中 又称为饱和度（saturations）。尽管毛细现象的方程中已经考虑了微观界面效应，但这个物理场不会跟踪描述出不同相之间的界面。

在 多孔介质相传递 物理场中，软件会默认添加【相和多孔介质传递属性】节点，以定义每个项的关键属性。由于是多相流，最简单的比如两个相 分别是空气 和 水，那么它们饱和度之和应该始终等于1：

在【多孔介质相传递】对象的设置里，【受体积约束】选项中，可以选择是哪一个项要根据体积约束条件计算得出。默认是s1.

多孔介质中的相传递与自由流动的一个重要不同点 就在于多孔介质中需要考虑相对渗透率（relative permeabilities）和毛细压力（capillary pressures）。这个纸芯吸水的案例本身就是用来展示毛细现象仿真的。

根据帮助文档，我们先来学习一下多孔介质相输运的控制方程。

软件界面上展示的瞬态方程为：

只能看出是一个对时间的偏导，加一个对空间的偏导。N_i具体是什么，写在了下附的【相和多孔介质传递属性】对象中。

回到文档里，方程写的更清楚一些（这就是连续性方程嘛）：

第一个 \epsilon_p，是无量纲的孔隙率，属于多孔介质本身的材料参数；\rho_i 当然就是密度了。s_i 是前面提到的待求的场变量，即某个相的体积分数。右边散度符号作用的对象中，u_i 这里应该被解释为第i相的体积通量（volumetric flux）。这一体积通量根据扩展的达西定律（extended Darcy's law）确定：

这么大一坨东西……别慌，前面那仨一看就知道是系数；后面的粗体g表示重力加速度矢量，也是常数（在界面上勾选“包含重力”才会考虑这一项）。真正的场变量只有p_i，某一相的压力。

总结到一起，方程写为（这张截图清晰一点，因为来自案例PDF文档）：

至于这么多系数分别是啥，只需要知道都是材料参数就行了。

在帮助文档里，写法略有不同。如果把某一个相的压力设为独立（就是前面设置了受体积约束那里的相，默认为s1），那么其他相的压力都是和它相比的相对压力。

上面的方程，梯度那里就要拆成两项：

到这，还是有很大问题。上面的控制方程，左边对时间求偏导的部分已经清楚了；右边对空间求偏导的部分，场变量却还没联系到相分数 s_i。我们知道，一个偏微分方程你别管前边跟了多少个系数，到最后总要写成对某一个统一的因变量分别求时间和空间偏导的形式。

这说明，在压力 p_i 和相分数 s_i 之间，还差着一个方程，把它俩联系起来。这种填补空缺的方程，在固体力学或传热学中一般都被称为本构方程，constitutive equations. 方程里的系数都被称为constitutive constant本构常数。

毛细压力模型

这个将相分数 s_i 和相压力 p_i 联系起来的关系，就是毛细压力模型。在Comsol软件界面上，如图所示，有四种可选。

用户定义咱就不提了。第二、第三两个模型，分别用提出模型的科学家的名字命名。诶你看看人家的研究工作，一篇论文青史留名。

排在前面的Genuchten论文写于1980年；排在后面的Brooks-Corey俩人论文写于1966年。都是“经典物理学”啊。

使用这些预定义毛细压力模型时，用户需要指定哪个相是润湿相。如果与体积约束相相同，则其他相压力使用加法，否则使用减法。

这里，毛细压力 p_c 具体的计算公式，由【湿相饱和度】（本例为s1）（见鬼，此处官方案例PDF里写错了，官方笔误写成了s2）和【入口毛细压力】p_ec 共同决定。Brooks-Corey公式嘛……实在很简单，就是一个指数形式的拟合。

\lambda_p 是孔隙分布指数，本例中取2. 这东西完全就是个拟合用的系数嘛。

但是，可别觉得Brooks和Corey俩人就做了这点工作。注意到这里 s_w头顶上有一横\bar，它可不直接是某个相分数，而是润湿相的有效相分数。另外，前面瞅着长得很像常数的相对渗透率 \kappa，其实也是这些相分数的函数。

这套公式，软件界面和PDF文档倒还算比较一致。

好。到此为止，虽然公式写的乱花渐欲迷人眼，但总算是把控制方程给封闭起来了。控制方程最终可以化简成只关于相分数 s_i 的瞬态偏微分方程。

达西定律 物理场

第二个物理场就是【达西定律】了。达西定律物理场，属于 多孔介质流动 模块。这个模块中的物理场接口用于描述多孔介质中的输运现象。

什么是多孔介质？

多孔材料是由固体（基质）和充满液体或气体的空腔（孔隙）组成的。从纳米材料到多孔反应器，从电子元件的冷却到大规模的岩土工程应用，多孔材料的尺寸和应用范围都很广泛。很多工业领域都会用到多孔介质，例如燃料电池、纸浆干燥、食品生产、过滤、混凝土、陶瓷、吸湿剂、纺织物等等。土壤是最常见的天然多孔介质，里面可以包含多相流动和输运。

通过多孔介质的流动

这些多孔介质的共同点是材料的宏观尺寸比微观平均孔径大得多，因此必须采用宏观方法建立模型。

在多孔介质中，最重要的两个参数是孔隙率和渗透率。孔隙率 \epsilon 描述了空洞/孔隙占总体积的比例。

而渗透率（permeability），\kappa，单位(m^2) 指定了流体通过多孔介质的能力。

达西定律

而 描述多孔介质中液体流动的基本定律，就是达西定律（Darcy's law）。它描述了速度场 u 与压力梯度之间的关系——就是乘系数的线性关系。

分子\kappa 前面说了是渗透率；分母\mu 当然是流体的动力粘度。p是孔隙压力。

达西定律只适用于非常低的速度，或低雷诺数的流动。当雷诺数 Re>10时，或在相对较高的Knudsen数 Kn>0.1下，达西定律不再有效。这时，过渡区域可以使用 non-Darcian Ergun equation；雷诺数Re>1000时，Ergun's equation 可以用Burke-Plummer方程来近似，描述湍流。好在，由于流体在孔隙中的摩擦阻力中损失了相当大的能量，多孔介质中的流速一般都非常低。

在多孔介质中，流体中的切应力对整体流动的影响通常可以忽略不计。由于多孔介质的跨尺度特征，不可能对每个孔进行单独建模，因此需要将多孔介质中的流体进行均匀化。达西定律与孔隙流体的连续性方程和状态方程一起，提供了一个完整的数学模型用于描述多孔介质流动。其中，主要的驱动力就是压力梯度。

达西定律物理场接口中，控制方程和上面很像，都是连续性方程的形式。

多孔介质中的多相流

多孔介质多相流，是Comsol中的一个内置多物理场耦合节点。它将【多孔介质相传递】和【达西定律】物理场耦合起来。流体的密度和动力粘度，是多个相的平均值。

其中，在【多孔介质相传递】物理场中，模型输入的压力由【达西定律】给出。

在【达西定律】物理场中，流体相的密度和动力粘度由多物理场耦合节点计算平均值给出。

嗯。其乐融融，很是和谐。

边界条件与网格

在【多孔介质相传递】物理场中，添加了重力条件（这没啥可说的）。由于在相传递中，指定了s1（即流体相）为根据体积约束计算得出，所以后面的边界条件只需针对相s2（即空气相）来定义。

所以，左右两侧和下边界，都定义“无通量”条件。因为假设左右两侧不允许流体流出，而下边界对于气体相来说，也是没有流入的。

在上边界，定义一个质量通量边界条件，描述气体相 s2的流出。

此处 th指的是厚度；而p_lm是达西定律物理场中的拉格朗日乘子（Lagrange multiplier）. 这个东西是啥我们稍后再讲。

在【达西定律】物理场中，两侧边界定义无流动条件；下表面压力为零；上表面压力是这样一个表达式。

这里的含义是 静水压力减去毛细孔压。通过乘上之前在模型中定义的阶跃函数，当纸条吸满水时，毛细孔压将被设置为 0。

这个模型的网格很有意思。因为纸条宽度方向上没有变化，所以宽度方向只划分了一个单元；而因为靠近下表面处吸水速率比较快，所以网格划分时使用了线性分布的尺寸函数。网格画出来长这样：

毕竟不是结构力学分析。不同的物理场对网格的要求还真是差别非常大。

弱约束与拉格朗日乘子

前面在多孔介质相传递物理场里，顶部的通量处定义了一个p_lm 拉格朗日乘子。这是什么东西呢？

相关文档页位于：

首先，在达西定律物理场顶端的压力处，勾选了“使用弱约束”。

对于第一类边界条件，即直接在边界上指定待定场变量数值的条件，如果是强约束，则有限元列式中会在网格节点上逐点应用约束，在每个节点上指定对应自由度的数值，从而在总刚度矩阵中消去相应的自由度；

而如果改用弱约束，则会以单元的形函数为权重，将Reaction terms（这里可以翻译为反力，但其实含义更广泛）直接显式的包含在方程组中。这里的拉格朗日乘子通常有明确的物理意义。例如，对于结构分析来说，如果在位移边界条件处指定了弱约束，那么拉格朗日乘子就相当于施加强制位移时在节点上的力；对于其他分析类型来说，其物理意义可能相当于通量。例如在这里，达西定律物理场里，显式地指定了上边界的压力，那么拉格朗日乘子就代表相对应的质量通量。

彩云小译机翻的稀烂。凑合看吧。

对于很多求解项目来说，反力或通量如果是用场变量求偏导得到，在边界上可能不够精确。而使用弱约束和拉格朗日乘子，能直接得到更精确的边界通量数值。

需要注意的是，在轴对称分析中，第一类（狄利克雷）边界条件对应的拉格朗日乘子并不等于单位面积的通量，而是等于单位长度*旋转一周（not equal to the reaction flux per area but rather per length and full revolution）。另外，在使用弱约束的拉格朗日乘子时，其单位可能与对应的通量单位不同，需要统一单位制后使用。

计算结果

最后来欣赏一波Comsol后处理画出的美丽云图。

达西定律求解出的压力场，t=600s

不同时刻的含水饱和度：

不得不说，Comsol这颜色表真好看，见一次爱一次。

JupiterAuroraBorealis，木星极光。多美的名字。偏个题，看NASA官网2004年发布的木星极光照片，是不是就这个配色。用它的反转色卡来描述吸水结果，真漂亮。

吸水量，解析解与数值解对比：

两条曲线的表达式分别为：

好。这个案例就到此为止。