基于VOF-DEM 方法的纖維懸浮液沉積特性研究＊

孫珊珊，朱立平，陳曉燕，吳仁民

（南京玻璃纖維研究設計院有限公司，南京210012）

0 前言

纖維懸浮液廣泛存在于造紙、電池隔板、保溫氈材等工業領域，如圖1 所示，它是一種由纖維、水、空氣等組成的氣—液—固三相共存的分散體系，其流動特性會影響纖維集合體的成形狀態，進而影響終端產品性能［1，2］。因此深入開展纖維懸浮液的流動特性的研究，對于相關生產過程的工藝參數優化具有重要意義。目前，纖維懸浮液流動特性的研究方法主要有理論分析、實驗分析和數值模擬3 種。由于纖維成形過程中懸浮液流動大多屬于較為復雜的湍流狀態，很難以單純的數學模型獲得纖維運動的解析解；實驗分析一直是研究流體流動狀況的最直接有效的方法，但實驗過程存在成本高、周期長及對測試儀器裝備要求高等問題，在一定程度上限制了實驗分析工作的深入。近年來，隨著計算機水平和數值模擬技術的快速進步，尤其是穩定可靠的商業軟件得到推廣應用，采用數值模擬技術解決纖維懸浮液流動問題越來越成為研究復雜流動過程的重要手段［3］。

圖1 纖維懸浮液氣液固三相分散體系

對纖維懸浮液流動特性的數值模擬研究最早出現在上世紀80年代初，Jeffery使用橢球懸浮液來模擬計算纖維懸浮液的運動方程，研究中做了許多假設，比如假設纖維對流場的干擾忽略不計，懸浮液被認為是牛頓流體等［4］；后來在非牛頓流體力學發展的推動下，使用數值仿真計算纖維懸浮液流動特性才有了較好的進展。由Ericksen首先提出的，后經許多學者修正的各向同性流體模型，被認為最適合于非牛頓流體懸浮液的本構方程［5］；Andersson的研究指出，當紙漿濃度在10%以上時，較低剪切速率的纖維類似于固體。當剪切速率達到一定值時，紙漿纖維懸浮液就和水的流動性質比較相似［6］。Baloch對漿料在縮小和擴大流道內的流動進行了模擬，使用Taylor Galerkin修正法進行了數值仿真［7］。范西?。?］等對纖維懸浮液攪拌狀態進行了數值仿真，但使用的是統計力學模型，導出的本構方程更適合于較高濃度的漿液。連璉運用液體-固體粒子兩相流理論及計算流體動力學方法，對液-固兩相流動特性進行了分析和研究，其推導的計算方法為液-固兩相流體的研究提供了有效的手段［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——顆?！w粒及顆?！诿娴南嗷ソ佑|作用力，計算方法詳見式（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 倍以上。