Star CCM+多孔介质仿真（二）——阻力系数的计算和设置

试验测量

经验公式

软件参数设置方法

本文主要解决的问题是如何确定多孔介质的惯性和粘性阻力系数。

试验测量

通过实验测出多孔介质的压降 $\Delta P$ 和流量 $Q$ 的多组数据；
根据流通面积将流量 $Q$ 换算成表观速度 $v$ ；
根据实验测量的多孔介质在流动方向上的长度L，将压降转换为压力梯度 $\frac{dP}{dx}=\frac{\Delta P}{L}$ ;
对 $\frac{dP}{dx}$ 和 $v$ 做二次拟合，得到速度一次方和二次方的系数；
根据Forchheimer定律 $\frac{dp}{dx}=av+bv^2$ ，其中 $a$ 即粘性项阻力系数， $b$ 即惯性项阻力系数；
对应的拟合系数就是所需要的阻力系数。

以上说明是以x方向为例，其他方向类似。

经验公式

对于多孔材料中的流动问题，目前采用较多的仍是以Forchheimer为代表的流动阻力模型：

$\frac{\Delta p}{L}=av+bv^2$

式中： $\Delta p /L$ 为压力梯度，v为材料的表观流速，a为粘性项阻力系数，b为惯性项阻力系数。

Ergun基于颗粒填充床的大量实验研究，对Forchheimer模型进行了修正，提出了著名的Ergun方程，此方程由于引入了颗粒直径、孔隙率以及流体物性参数而被广泛接受和验证，并被发展为所谓的Ergun型方程：

$\frac{\Delta p}{L}=A \frac{(1 - \epsilon )^2 \mu }{\epsilon^3 d_p^2} v+B \frac{(1-\epsilon)}{\epsilon^3d_p}\rho v^2$

式中：A为粘性项无量纲系数，B为惯性项无量纲系数， $\epsilon$ 为材料孔隙率， $d_p$ 为颗粒直径。

关于经验系数A、B，不同的研究者给出了各自的经验值，如Ergun给出A=175、B=1.75；Macdonald给出A=180、B=1.8；Irmay给出A=180、B=0.6等。

软件参数设置方法

实验测量可以直接得到粘性项和惯性项阻力系数，其设置方法可参考Star CCM+多孔介质仿真（一）。

对于通过经验公式得到的阻力系数，一种方法是自己先根据材料孔隙率和颗粒直径等参数直接把各阻力系数算出来，然后再添加到软件中，方法同上。

下面是用Star CCM+的Field Function功能来进行参数设置的过程。

Field Function是Star CCM+中存取求解器内单元和边界数据的机制，在Star CCM+中应用非常广泛。虽然使用Field Function没有直接写程序那么灵活，但无须安装特别的编译器，其功能强大、简单易用的特点给工程师们提供了很大的操作空间。

已知参数：

材料孔隙率 $\epsilon=0.3$

颗粒直径 $d_p=1mm$

参数设置步骤：

（1）定义材料孔隙率

在Tools节点下展开Field Functions节点，右键Field Functions节点，选择new -> Scalar左键确定，创建一个新的Field Function。

在Field Functions节点下生成一个名为“User Field Function 1”的新函数，重命名为“Bed Porosity”，把Properties对话框中的Function Name也重命名为“Bed Porosity”，Definition改为0.3，孔隙率是无量纲，Dimension默认就是无量纲的（Dimensionless），不用改。

（2）定义颗粒直径

操作方法与上面相同，名字为“Bed Sphere Diameter”，量纲为“Length”，数值定义为0.001（m）。

（3）定义粘性项阻力系数

同样的方法创建一个新的scalar函数，重命名为“Ergun Viscous Term”，量纲改为 $Mass \cdot Length^{-1} \cdot Time^{-1}$ ，公式定义为

175.0*$DynamicViscosity*pow((1-$BedPorosity),2)/(pow($BedPorosity,3)*pow($BedSphereDiameter,2))

（注：这里经验系数A取175）

输入公式时，已存在的变量或函数会再你输入前几个字面后在下拉列表中显示，如图。然后从列表中选择相应的变量或函数

最终输入的公式如下图，

如果变量写错了，不存在的变量会显示成红色，方便检查错误，比如下面的DynamicViscosity就写错了。一般从列表选是不会出错的。

公式输入完成后，点击OK，这样粘性项就定义好了。

（4）定义惯性项阻力系数

方法同（3），重命名为“Ergun Inertial Term”，量纲为 $Mass \cdot Length^{-4}$ ，定义公式为

(1.75*$Density*(1-$BedPorosity))/ (pow($BedPorosity,3)*$BedSphereDiameter)

（注：这里经验系数B取1.75）

（5）将定义好的阻力系数添加到Region的多孔介质中

选择Region -> Porous -> Physics Values -> Porous Inertial Resistance -> Isotropic Tensor -> Isotropic Component，在Method中选择Field Function

单击Scalar Function右边的<Select Function>，然后弹出一个对话框，在对话框中搜索“ergun”，相关的函数就显示出来了，选择Ergun Inertial Term函数，然后单击OK。

惯性项阻力系数就设置好了

粘性项阻力系数的设置方法相同

材料孔隙率也可以用Field Function的方法设置，如下

另外，多孔材料还有一个参数——迂曲率（tortuosity），软件在计算扩散系数时会用到，和孔隙率一样是输入数据。迂曲率也可以像孔隙率一样用Field Function进行设置，其计算可以根据经验公式用孔隙率来计算详见文献【4】，截图如下。

参考文献：

张新铭, 凌娅, 谷沁洋. 多孔石墨泡沫材料内流动阻力的扩展Ergun方程[J]. 材料导报, 2012, 26(4).
王补宣. 多孔介质中单相对流换热分析的流体渗流模式[J]. 上海交通大学学报, 1999(08):52-55.
Zeng Z , Grigg R . A Criterion for Non-Darcy Flow in Porous Media[J]. Transport in Porous Media, 2006, 63(1):p.57-69.
李滔, 李闽, 张烈辉, et al. 微多孔介质迂曲度与孔隙结构关系[J]. 天然气地球科学, 2018(8):1181-1189.
Star CCM+官方帮助文档。