基于VOF-DEM 方法的纤维悬浮液沉积特性研究＊

孙珊珊，朱立平，陈晓燕，吴仁民

（南京玻璃纤维研究设计院有限公司，南京210012）

0 前言

纤维悬浮液广泛存在于造纸、电池隔板、保温毡材等工业领域，如图1 所示，它是一种由纤维、水、空气等组成的气—液—固三相共存的分散体系，其流动特性会影响纤维集合体的成形状态，进而影响终端产品性能［1，2］。因此深入开展纤维悬浮液的流动特性的研究，对于相关生产过程的工艺参数优化具有重要意义。目前，纤维悬浮液流动特性的研究方法主要有理论分析、实验分析和数值模拟3 种。由于纤维成形过程中悬浮液流动大多属于较为复杂的湍流状态，很难以单纯的数学模型获得纤维运动的解析解；实验分析一直是研究流体流动状况的最直接有效的方法，但实验过程存在成本高、周期长及对测试仪器装备要求高等问题，在一定程度上限制了实验分析工作的深入。近年来，随着计算机水平和数值模拟技术的快速进步，尤其是稳定可靠的商业软件得到推广应用，采用数值模拟技术解决纤维悬浮液流动问题越来越成为研究复杂流动过程的重要手段［3］。

图1 纤维悬浮液气液固三相分散体系

对纤维悬浮液流动特性的数值模拟研究最早出现在上世纪80年代初，Jeffery使用椭球悬浮液来模拟计算纤维悬浮液的运动方程，研究中做了许多假设，比如假设纤维对流场的干扰忽略不计，悬浮液被认为是牛顿流体等［4］；后来在非牛顿流体力学发展的推动下，使用数值仿真计算纤维悬浮液流动特性才有了较好的进展。由Ericksen首先提出的，后经许多学者修正的各向同性流体模型，被认为最适合于非牛顿流体悬浮液的本构方程［5］；Andersson的研究指出，当纸浆浓度在10%以上时，较低剪切速率的纤维类似于固体。当剪切速率达到一定值时，纸浆纤维悬浮液就和水的流动性质比较相似［6］。Baloch对浆料在缩小和扩大流道内的流动进行了模拟，使用Taylor Galerkin修正法进行了数值仿真［7］。范西俊［8］等对纤维悬浮液搅拌状态进行了数值仿真，但使用的是统计力学模型，导出的本构方程更适合于较高浓度的浆液。连琏运用液体-固体粒子两相流理论及计算流体动力学方法，对液-固两相流动特性进行了分析和研究，其推导的计算方法为液-固两相流体的研究提供了有效的手段［9］。Olson研究了悬浮液湍流流动过程中的颗粒物运动特性，得到的方程获得了较广泛的应用［10］。Perazzo A等使用微流体法来生产高纵横比、柔性微纤维的均匀浓缩悬浮液，并展示了这种微纤维悬浮液的剪切增厚和胶凝行为［11］。Moosaie A提出一种双向耦合拉格朗日力矩近似方法，用于模拟湍流中的布朗纤维悬浮液。流动方程采用欧拉方法求解，采用非牛顿应力张量来考虑纤维对流体流动的影响，该方法可适用于无惯性微纤维的稀相悬浮流［12］。

总体来看，目前国内外对纤维与悬浮液之间的作用机理研究虽取得一系列成果，但仍有许多问题未得到合理解释，尤其是在理论建模方法层面尚不成熟，这在一定程度上限制了悬浮液流动过程中纤维取向分布及其流动特性的有效控制。本文将构建一种基于VOF-DEM的数值模拟方法，并通过实验验证其可行性，这对于研究和理解悬浮液中纤维与流场间的作用机理、纤维集合体成形过程的工艺优化将具有一定的理论意义和工程价值。

1 数学模型

针对纤维悬浮液的气液固三相流动过程，由于气泡及固相纤维体积分数相对较少（低浓），忽略悬浮液内气泡对纤维运动的影响，主要考虑气-液、液-固之间的相互作用［13，14］。

1.1 流体-颗粒系统

流体-颗粒系统采用欧拉—拉格朗日模型来表述，网格内固相颗粒的体积分数计算公式为：

式中，

εp——固相颗粒体积分数；

Vp——固相颗粒体积，m3；

Vf——流体相体积，m3。

1.2 颗粒相受力运动模型

本文以玻璃纤维为研究对象，该纤维是一类柔性丝状非球形颗粒：颗粒细长，形状在轴向上占主导；柔软，易绕曲，变形；它们的受力、转动和取向分布复杂，建模难度较大。为了简化计算，将其假设为均质的刚性细长体的组合，其质心就是几何中心，刚性颗粒间使用铰链进行连接，可用于表征纤维的柔性状态，如图2 所示。

图2 纤维颗粒的链式模型

计算域内所有的颗粒运动仿真基于拉格朗日框架，满足牛顿运动定律［15，16］，其运动表达式为式（2）和式（3）所示：

式中：

mp——颗粒质量，kg；

vp——颗粒速度，m/s；

g ——重力加速度矢量，m/s2；

FD——流体对颗粒的作用力，计算方法详见式（7）；

FC——颗粒—颗粒及颗粒—壁面的相互接触作用力，计算方法详见式（4）。

式中：

Ip——惯性动量，kg·m2；

ω——角速度，rad/s；

Tc——接触力扭矩，N·m。

式中：

FCij——颗粒i和颗粒j的碰撞作用力，N；

FCn，ij——颗粒i和颗粒j的法向作用力，N，计算方法详见式（5）；

FCt，ij——颗粒i和颗粒j的切向作用力，N，计算方法详见式（6）。

式中：

kn——颗粒的法向刚度，N/m；

dn——颗粒法向相撞所产生的弹性形变，m；

ηn——法向阻尼系数，N·s/m；

vn——法向碰撞相对速度，m/s。

式中：

kt——颗粒的切向刚度，N/m；

dt——颗粒切向相撞所产生的弹性形变，m；

ηt——切向阻尼系数，N·s/m；

vt——切向碰撞相对速度，m/s。

式中：

CD——颗粒曳力系数，由当量球形颗粒的阻力系数修正得到；

ug——气相场速度，m/s；

A ——杆状颗粒段的截面积，m2；

ρ——杆状颗粒段的密度，kg/m3；

us——杆状颗粒段的速度，m/s。

1.3 气—液界面与流动模型

基于气-液两相不发生互相穿插现象，采用VOF描述气液界面的多相流动。对增加到模型里的每一附加相，就引进一个变量：即计算单元里相的容积比率［17］。

（1）容积比率方程（连续相方程）

跟踪相之间的界面是通过求解一相或多相的容积比率的连续方程来完成的。对第 q 相，方程如下：

（2）动量方程

通过求解整个区域内的单一的动量方程，作为结果的速度场是由各相共享的。动量方程取决于通过属性 ρ 和 µ 的所有相的容积比率。

2 模拟对象描述

以某纤维毡材的湿法成形工艺过程为研究对象，浆液由装置上方进入形成均匀分散的纤维悬浮液，装置的底部布置透水滤网，浆液通过滤网的漓水作用使纤维在网的表面沉积，滤网下方布置真空抽吸装置，以增强滤网处的液、固分离作用，最终得到纤维沉积层。实验设备结构如图3 所示，该设备整体为圆柱状，直径200 mm，总高度325 mm，上部为纤维悬浮液空间，高度245 mm，中部为滤网，下部为抽真空区域，高度80 mm，底部为脱水出口，直径50 mm。

图3 湿法成形实验示意图

根据湿法成形设备结构参数建立相应的几何模型，如图4 所示。对滤网的处理，采用多孔介质模型，其过滤过程涉及渗透动力学和多孔介质理论，流动过程遵循Darcy定律。初始阶段纤维悬浮液填充于滤网上部空间，纤维均匀分布在悬浮液内部，上部空间与大气接触面设置为压力边界条件，出口设置为流量出口。计算域的网格划分结果如图5 所示，其中滤网表面的网格进行加密处理，经网格无关性检测后确定网格数量为54 万个。初始状态下假设纤维在流体中随机均匀分布，如图6 所示。计算区域的离散化采用有限体积法，控制方程的离散化全部采用一阶迎风格式，速度场和压力场的耦合解法采用经典的SIMPLE算法，采用Ansys-Fluent软件作求解器，颗粒相的运动状态通过二次开发程序进行求解。

图4 几何结构模型

图5 模型的网格划分

图6 纤维悬浮液初始状态

3 模拟结果分析与验证

图7为经气—液—固三相流动模拟后的纤维沉积效果图，（a）和（b）分别为纤维集合体的主视和俯视图。从图7（a）来看，纤维沉积厚度从中心向边缘呈现逐渐下降的变化趋势，中心与边缘的厚度差异超过1倍，这是由于脱水出口位于设备中部，使得悬浮液在负压抽吸作用下向中心聚集，进而导致了纤维的不均匀分布。通过图7（a）的局部放大图可观察到纤维主流朝向与水平面平行，而7（b）中则呈现了纤维在水平面上的分布状态，在各角度上呈现无序分布状态，经交织后形成了大小不一的微孔。

图7 纤维沉积模拟效果图

为验证数值模拟方法的合理性，本文设计了小规模悬浮液流动实验，首先将充分疏解后的玻璃纤维浆液置入实验设备，后将经脱水后沉积在滤网上的纤维沉积层样本进行烘干，测试圆形样本径向的厚度分布数据，用于模型验证。为了降低实验过程中各种因素波动的影响，通过多次实验取平均值来减少误差，图8 是烘干后的纤维沉积层样本。图9 为样本厚度测试图，以同心圆上的多点厚度平均值作为该径向位置的纤维厚度，图9a中包括6 种径向位置，分别称为一环、二环、三环、四环、五环、六环，每个径向通过测厚仪测量四个位点的厚度，取其平均值作为该径向位置的厚度值，重复上述步骤获得10 组实验样本厚度分布情况取径向厚度分布平均值。

图8 烘干后的纤维样本

图9 样本厚度测试

图10 为各样本径向位置厚度分布曲线图，从图10 可以看出，大部分样本都呈现环数越大厚度越小的情况，某些样本中存在个别不符合单调性的点，这可能和实验误差有关。平均值曲线表现出离抄片圆心越远厚度越小的现象，因此得到这样的结论：随机分布的纤维随着脱水过程沉积在滤网上，中心处纤维沉积最多，从中心向边缘扩展，纤维沉积厚度呈现逐渐下降的趋势，这在定性上与模拟结果一致。

图10 各样本径向位置厚度分布曲线

本文对模拟结果的沉积纤维层径向厚度进行计算，求取每个径向位置处厚度的平均值，并与实验结果进行对比，数据如图11 所示。从图中可以看出，模拟的厚度分布结果与实验具有相同变化趋势，从中心到边缘处，沉积层的厚度逐渐降低。从厚度分布数据来看，模拟结果与实验结果最大偏差为13%，与实验结果较为吻合，表明所建立的数值模拟方法的可靠性。

图11 纤维沉积层厚度分布模拟验证

4 结论

本文构建了纤维悬浮液气—液—固三相耦合流动数学模型，并通过纤维毡材的成形实验对模拟方法进行了验证。主要结论如下：

（1）基于VOF方法描述气液两相流动，采用DEM方法求解纤维运动，能够较为准确地计算出纤维悬浮液的三相流动过程，模拟结果与实验结果达到规律性一致，最大量化偏差为13%。

（2）纤维毡材成形实验中，初始空间分布均匀的纤维，在文中所采取的实验设备条件下，由于脱水出口位于设备中部，经沉降后纤维厚度分布从中心向边缘呈现逐渐下降的变化趋势，两者的最大差异达1 倍以上。