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

（南京玻璃纖維研究設(shè)計(jì)院有限公司，南京210012）

0 前言

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

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

對(duì)纖維懸浮液流動(dòng)特性的數(shù)值模擬研究最早出現(xiàn)在上世紀(jì)80年代初，Jeffery使用橢球懸浮液來模擬計(jì)算纖維懸浮液的運(yùn)動(dòng)方程，研究中做了許多假設(shè)，比如假設(shè)纖維對(duì)流場(chǎng)的干擾忽略不計(jì)，懸浮液被認(rèn)為是牛頓流體等［4］；后來在非牛頓流體力學(xué)發(fā)展的推動(dòng)下，使用數(shù)值仿真計(jì)算纖維懸浮液流動(dòng)特性才有了較好的進(jìn)展。由Ericksen首先提出的，后經(jīng)許多學(xué)者修正的各向同性流體模型，被認(rèn)為最適合于非牛頓流體懸浮液的本構(gòu)方程［5］；Andersson的研究指出，當(dāng)紙漿濃度在10%以上時(shí)，較低剪切速率的纖維類似于固體。當(dāng)剪切速率達(dá)到一定值時(shí)，紙漿纖維懸浮液就和水的流動(dòng)性質(zhì)比較相似［6］。Baloch對(duì)漿料在縮小和擴(kuò)大流道內(nèi)的流動(dòng)進(jìn)行了模擬，使用Taylor Galerkin修正法進(jìn)行了數(shù)值仿真［7］。范西?。?］等對(duì)纖維懸浮液攪拌狀態(tài)進(jìn)行了數(shù)值仿真，但使用的是統(tǒng)計(jì)力學(xué)模型，導(dǎo)出的本構(gòu)方程更適合于較高濃度的漿液。連璉運(yùn)用液體-固體粒子兩相流理論及計(jì)算流體動(dòng)力學(xué)方法，對(duì)液-固兩相流動(dòng)特性進(jìn)行了分析和研究，其推導(dǎo)的計(jì)算方法為液-固兩相流體的研究提供了有效的手段［9］。Olson研究了懸浮液湍流流動(dòng)過程中的顆粒物運(yùn)動(dòng)特性，得到的方程獲得了較廣泛的應(yīng)用［10］。Perazzo A等使用微流體法來生產(chǎn)高縱橫比、柔性微纖維的均勻濃縮懸浮液，并展示了這種微纖維懸浮液的剪切增厚和膠凝行為［11］。Moosaie A提出一種雙向耦合拉格朗日力矩近似方法，用于模擬湍流中的布朗纖維懸浮液。流動(dòng)方程采用歐拉方法求解，采用非牛頓應(yīng)力張量來考慮纖維對(duì)流體流動(dòng)的影響，該方法可適用于無慣性微纖維的稀相懸浮流［12］。

總體來看，目前國(guó)內(nèi)外對(duì)纖維與懸浮液之間的作用機(jī)理研究雖取得一系列成果，但仍有許多問題未得到合理解釋，尤其是在理論建模方法層面尚不成熟，這在一定程度上限制了懸浮液流動(dòng)過程中纖維取向分布及其流動(dòng)特性的有效控制。本文將構(gòu)建一種基于VOF-DEM的數(shù)值模擬方法，并通過實(shí)驗(yàn)驗(yàn)證其可行性，這對(duì)于研究和理解懸浮液中纖維與流場(chǎng)間的作用機(jī)理、纖維集合體成形過程的工藝優(yōu)化將具有一定的理論意義和工程價(jià)值。

1 數(shù)學(xué)模型

針對(duì)纖維懸浮液的氣液固三相流動(dòng)過程，由于氣泡及固相纖維體積分?jǐn)?shù)相對(duì)較少（低濃），忽略懸浮液內(nèi)氣泡對(duì)纖維運(yùn)動(dòng)的影響，主要考慮氣-液、液-固之間的相互作用［13，14］。

1.1 流體-顆粒系統(tǒng)

流體-顆粒系統(tǒng)采用歐拉—拉格朗日模型來表述，網(wǎng)格內(nèi)固相顆粒的體積分?jǐn)?shù)計(jì)算公式為：

式中，

εp——固相顆粒體積分?jǐn)?shù)；

Vp——固相顆粒體積，m3；

Vf——流體相體積，m3。

1.2 顆粒相受力運(yùn)動(dòng)模型

本文以玻璃纖維為研究對(duì)象，該纖維是一類柔性絲狀非球形顆粒：顆粒細(xì)長(zhǎng)，形狀在軸向上占主導(dǎo)；柔軟，易繞曲，變形；它們的受力、轉(zhuǎn)動(dòng)和取向分布復(fù)雜，建模難度較大。為了簡(jiǎn)化計(jì)算，將其假設(shè)為均質(zhì)的剛性細(xì)長(zhǎng)體的組合，其質(zhì)心就是幾何中心，剛性顆粒間使用鉸鏈進(jìn)行連接，可用于表征纖維的柔性狀態(tài)，如圖2 所示。

圖2 纖維顆粒的鏈?zhǔn)侥Ｐ?/p>

計(jì)算域內(nèi)所有的顆粒運(yùn)動(dòng)仿真基于拉格朗日框架，滿足牛頓運(yùn)動(dòng)定律［15，16］，其運(yùn)動(dòng)表達(dá)式為式（2）和式（3）所示：

式中：

mp——顆粒質(zhì)量，kg；

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

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

FD——流體對(duì)顆粒的作用力，計(jì)算方法詳見式（7）；

FC——顆?！w粒及顆粒—壁面的相互接觸作用力，計(jì)算方法詳見式（4）。

式中：

Ip——慣性動(dòng)量，kg·m2；

ω——角速度，rad/s；

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

式中：

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

FCn，ij——顆粒i和顆粒j的法向作用力，N，計(jì)算方法詳見式（5）；

FCt，ij——顆粒i和顆粒j的切向作用力，N，計(jì)算方法詳見式（6）。

式中：

kn——顆粒的法向剛度，N/m；

dn——顆粒法向相撞所產(chǎn)生的彈性形變，m；

ηn——法向阻尼系數(shù)，N·s/m；

vn——法向碰撞相對(duì)速度，m/s。

式中：

kt——顆粒的切向剛度，N/m；

dt——顆粒切向相撞所產(chǎn)生的彈性形變，m；

ηt——切向阻尼系數(shù)，N·s/m；

vt——切向碰撞相對(duì)速度，m/s。

式中：

CD——顆粒曳力系數(shù)，由當(dāng)量球形顆粒的阻力系數(shù)修正得到；

ug——?dú)庀鄨?chǎng)速度，m/s；

A ——桿狀顆粒段的截面積，m2；

ρ——桿狀顆粒段的密度，kg/m3；

us——桿狀顆粒段的速度，m/s。

1.3 氣—液界面與流動(dòng)模型

基于氣-液兩相不發(fā)生互相穿插現(xiàn)象，采用VOF描述氣液界面的多相流動(dòng)。對(duì)增加到模型里的每一附加相，就引進(jìn)一個(gè)變量：即計(jì)算單元里相的容積比率［17］。

（1）容積比率方程（連續(xù)相方程）

跟蹤相之間的界面是通過求解一相或多相的容積比率的連續(xù)方程來完成的。對(duì)第 q 相，方程如下：

（2）動(dòng)量方程

通過求解整個(gè)區(qū)域內(nèi)的單一的動(dòng)量方程，作為結(jié)果的速度場(chǎng)是由各相共享的。動(dòng)量方程取決于通過屬性 ρ 和 μ 的所有相的容積比率。

2 模擬對(duì)象描述

以某纖維氈材的濕法成形工藝過程為研究對(duì)象，漿液由裝置上方進(jìn)入形成均勻分散的纖維懸浮液，裝置的底部布置透水濾網(wǎng)，漿液通過濾網(wǎng)的漓水作用使纖維在網(wǎng)的表面沉積，濾網(wǎng)下方布置真空抽吸裝置，以增強(qiáng)濾網(wǎng)處的液、固分離作用，最終得到纖維沉積層。實(shí)驗(yàn)設(shè)備結(jié)構(gòu)如圖3 所示，該設(shè)備整體為圓柱狀，直徑200 mm，總高度325 mm，上部為纖維懸浮液空間，高度245 mm，中部為濾網(wǎng)，下部為抽真空區(qū)域，高度80 mm，底部為脫水出口，直徑50 mm。

圖3 濕法成形實(shí)驗(yàn)示意圖

根據(jù)濕法成形設(shè)備結(jié)構(gòu)參數(shù)建立相應(yīng)的幾何模型，如圖4 所示。對(duì)濾網(wǎng)的處理，采用多孔介質(zhì)模型，其過濾過程涉及滲透動(dòng)力學(xué)和多孔介質(zhì)理論，流動(dòng)過程遵循Darcy定律。初始階段纖維懸浮液填充于濾網(wǎng)上部空間，纖維均勻分布在懸浮液內(nèi)部，上部空間與大氣接觸面設(shè)置為壓力邊界條件，出口設(shè)置為流量出口。計(jì)算域的網(wǎng)格劃分結(jié)果如圖5 所示，其中濾網(wǎng)表面的網(wǎng)格進(jìn)行加密處理，經(jīng)網(wǎng)格無關(guān)性檢測(cè)后確定網(wǎng)格數(shù)量為54 萬個(gè)。初始狀態(tài)下假設(shè)纖維在流體中隨機(jī)均勻分布，如圖6 所示。計(jì)算區(qū)域的離散化采用有限體積法，控制方程的離散化全部采用一階迎風(fēng)格式，速度場(chǎng)和壓力場(chǎng)的耦合解法采用經(jīng)典的SIMPLE算法，采用Ansys-Fluent軟件作求解器，顆粒相的運(yùn)動(dòng)狀態(tài)通過二次開發(fā)程序進(jìn)行求解。

圖4 幾何結(jié)構(gòu)模型

圖5 模型的網(wǎng)格劃分

圖6 纖維懸浮液初始狀態(tài)

3 模擬結(jié)果分析與驗(yàn)證

圖7為經(jīng)氣—液—固三相流動(dòng)模擬后的纖維沉積效果圖，（a）和（b）分別為纖維集合體的主視和俯視圖。從圖7（a）來看，纖維沉積厚度從中心向邊緣呈現(xiàn)逐漸下降的變化趨勢(shì)，中心與邊緣的厚度差異超過1倍，這是由于脫水出口位于設(shè)備中部，使得懸浮液在負(fù)壓抽吸作用下向中心聚集，進(jìn)而導(dǎo)致了纖維的不均勻分布。通過圖7（a）的局部放大圖可觀察到纖維主流朝向與水平面平行，而7（b）中則呈現(xiàn)了纖維在水平面上的分布狀態(tài)，在各角度上呈現(xiàn)無序分布狀態(tài)，經(jīng)交織后形成了大小不一的微孔。

圖7 纖維沉積模擬效果圖

為驗(yàn)證數(shù)值模擬方法的合理性，本文設(shè)計(jì)了小規(guī)模懸浮液流動(dòng)實(shí)驗(yàn)，首先將充分疏解后的玻璃纖維漿液置入實(shí)驗(yàn)設(shè)備，后將經(jīng)脫水后沉積在濾網(wǎng)上的纖維沉積層樣本進(jìn)行烘干，測(cè)試圓形樣本徑向的厚度分布數(shù)據(jù)，用于模型驗(yàn)證。為了降低實(shí)驗(yàn)過程中各種因素波動(dòng)的影響，通過多次實(shí)驗(yàn)取平均值來減少誤差，圖8 是烘干后的纖維沉積層樣本。圖9 為樣本厚度測(cè)試圖，以同心圓上的多點(diǎn)厚度平均值作為該徑向位置的纖維厚度，圖9a中包括6 種徑向位置，分別稱為一環(huán)、二環(huán)、三環(huán)、四環(huán)、五環(huán)、六環(huán)，每個(gè)徑向通過測(cè)厚儀測(cè)量四個(gè)位點(diǎn)的厚度，取其平均值作為該徑向位置的厚度值，重復(fù)上述步驟獲得10 組實(shí)驗(yàn)樣本厚度分布情況取徑向厚度分布平均值。

圖8 烘干后的纖維樣本

圖9 樣本厚度測(cè)試

圖10 為各樣本徑向位置厚度分布曲線圖，從圖10 可以看出，大部分樣本都呈現(xiàn)環(huán)數(shù)越大厚度越小的情況，某些樣本中存在個(gè)別不符合單調(diào)性的點(diǎn)，這可能和實(shí)驗(yàn)誤差有關(guān)。平均值曲線表現(xiàn)出離抄片圓心越遠(yuǎn)厚度越小的現(xiàn)象，因此得到這樣的結(jié)論：隨機(jī)分布的纖維隨著脫水過程沉積在濾網(wǎng)上，中心處纖維沉積最多，從中心向邊緣擴(kuò)展，纖維沉積厚度呈現(xiàn)逐漸下降的趨勢(shì)，這在定性上與模擬結(jié)果一致。

圖10 各樣本徑向位置厚度分布曲線

本文對(duì)模擬結(jié)果的沉積纖維層徑向厚度進(jìn)行計(jì)算，求取每個(gè)徑向位置處厚度的平均值，并與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，數(shù)據(jù)如圖11 所示。從圖中可以看出，模擬的厚度分布結(jié)果與實(shí)驗(yàn)具有相同變化趨勢(shì)，從中心到邊緣處，沉積層的厚度逐漸降低。從厚度分布數(shù)據(jù)來看，模擬結(jié)果與實(shí)驗(yàn)結(jié)果最大偏差為13%，與實(shí)驗(yàn)結(jié)果較為吻合，表明所建立的數(shù)值模擬方法的可靠性。

圖11 纖維沉積層厚度分布模擬驗(yàn)證

4 結(jié)論

本文構(gòu)建了纖維懸浮液氣—液—固三相耦合流動(dòng)數(shù)學(xué)模型，并通過纖維氈材的成形實(shí)驗(yàn)對(duì)模擬方法進(jìn)行了驗(yàn)證。主要結(jié)論如下：

（1）基于VOF方法描述氣液兩相流動(dòng)，采用DEM方法求解纖維運(yùn)動(dòng)，能夠較為準(zhǔn)確地計(jì)算出纖維懸浮液的三相流動(dòng)過程，模擬結(jié)果與實(shí)驗(yàn)結(jié)果達(dá)到規(guī)律性一致，最大量化偏差為13%。

（2）纖維氈材成形實(shí)驗(yàn)中，初始空間分布均勻的纖維，在文中所采取的實(shí)驗(yàn)設(shè)備條件下，由于脫水出口位于設(shè)備中部，經(jīng)沉降后纖維厚度分布從中心向邊緣呈現(xiàn)逐漸下降的變化趨勢(shì)，兩者的最大差異達(dá)1 倍以上。