基于非牛頓流體本構方程的熔噴纖維直徑預測

孫光武， 李杰聰， 辛三法， 王新厚

(1. 上海工程技術大學 服裝學院， 上海 201620； 2. 東華大學 紡織學院， 上海 201620)

熔噴技術是由高溫高速氣流牽伸熔融聚合物進而一步成纖的紡絲技術，制備出的纖維(一般為1～5 μm)在成網簾上逐漸聚集而形成熔噴非織造布。由于纖維較細，具有較大的比表面積，熔噴非織造布往往被用于制造空氣和液體過濾材料、吸油材料等。熔噴非織造布在應用中的表現直接受到其纖維直徑的影響，具有較細纖維的熔噴非織造布的過濾性能往往較好。熔噴非織造生產線的工藝參數對熔噴非織造纖維直徑具有很大的影響，需要反復調整才可獲得最佳的纖維細度，因此，熔噴拉伸理論模型的研究特別是涉及到纖維直徑的預測則顯得愈發重要。因為其能在理論上預測直徑，從而指導實際生產，并能從機制上闡述各工藝參數之間的關系，為商業化生產提供理論依據。

早期的熔噴拉伸理論由Uyttendaele等[1]通過改進熔融紡絲模型得到，這種模型被稱為“自由端”纖維-氣流模型。因在熔噴技術中氣流牽伸纖維，纖維的尾端并無握持力，這與熔融紡絲“非自由端”不同。文獻[2-3]將熔噴拉伸模型分別推廣至二維和三維空間。各類方法如蒙特卡洛法、龍格-庫塔法、有限元法和冪律模型[4]也被廣泛用于求解熔噴拉伸模型方程。Sinha等[5]提出一種簡單的準一維模型，該模型計算量較小，同時也能精確預測纖維的運動和直徑變化。Shambaugh等[6]以靜態結晶動力學方程為基礎，模擬了熔噴纖維在成纖過程中的結晶變化，該方法成功將熔噴纖維的宏觀結構預測與微觀結構預測相結合。不同于上述模型，Sun等[7-8]提出了一種基于拉格朗日方法的拉伸模型，通過將纖維假設成由珠子、彈簧、黏壺串聯的模型，預測每個珠子的位置和物理性能的變化，進而預測纖維的實時位置和各項性能參數。

熔噴纖維拉伸模型經過多年的發展，從簡單到復雜，從一維到三維，目前文獻中已報道的模型可預測任意時刻纖維在整個氣流場中的位置、溫度、速度、直徑，甚至結晶度。然而為了簡化計算，在模型中采用牛頓流體本構方程或者經典的Maxwell本構方程，這些方程由于難以準確描述聚合物熔體復雜的非牛頓黏彈行為，導致模擬的纖維直徑與實際結果差異較大。為此，本文引入4種非牛頓流體本構方程：PTT(Phan-Thien-Tanner)方程、UCM(Upper-Convected Maxwell)方程、Giesekus方程和Rouse-Zimm方程，結合已建立的拉格朗日熔噴拉伸模型[9-11]預測纖維直徑，以定量化分析非牛頓流體本構方程對最終纖維直徑預測結果的影響[12]，從機制上揭示非牛頓流體的黏彈行為對熔噴纖維的成纖作用。

1 理論模型

拉格朗日方法是一種描述流體運動的基本方法，即對流體質點進行標記，流體質點的空間坐標與物理量就表示為質點及時間的函數。建立熔噴纖維拉伸模型的目的是為了研究熔噴纖維在拉伸過程中的細化規律、纖維內外應力變化、纖維在氣流場中的運動情況以及纖維段之間的相互作用。模型的研究對象是纖維片段，其目的是求解出纖維片段的各物理量與空間坐標在流場中隨時間變化的規律，而這與拉格朗日方法一致。

圖1示出采用珠鏈模型描述的熔噴纖維。纖維可假設為由珠子和串聯珠子的熔體單元(類似鏈條)組成。

圖1 纖維珠鏈模型Fig.1 Bead-chain model of fiber

在不考慮纖維相互糾纏的前提下，熔體從噴絲孔逐漸細化落至成網簾上，纖維所受到的作用力主要有4種：氣流拉伸力、黏彈力、表面張力和重力。其中：氣流拉伸力和重力為外應力，主要起推動纖維運動、拉伸并細化纖維的作用；而表面張力和黏彈力屬內應力，主要起抵抗纖維拉伸和彎曲形變的作用。由于整體的動力學模型較為復雜，且除黏彈力之外的其他作用力并非本文的研究重點，在此僅列出數理方程，詳細推導過程可參見文獻[9-11]。

氣流拉伸力是纖維在高速流場中受到的主要作用力，氣流力的大小直接影響纖維在氣流中的運動和最終成纖時的纖維直徑。纖維段(i-1,i)上所受到的氣流力Fd,i-1,i可以分解成壓差阻力Fp,i-1,i和摩擦阻力Ff,i-1,i，有

Fd,i-1,i=Ff,i-1,i+Fp,i-1,i

(1)

(2)

(3)

式中：Cp、Cf分別為壓差阻力系數、摩擦阻力系數；ρa為空氣密度，其值為1.29 kg/m3；vrt、vrn為纖維相對速度(纖維實際速度與所在位置氣流速度的差)沿其軸向和徑向的分量，m/s；di-1,i為纖維段(i-1,i)的直徑，m；li-1，i為纖維段(i-1,i)的長度，也是珠子i-1與i的間距，m；ft、fn分別為纖維軸向單位向量、徑向單位向量。

表面張力與纖維的氣流力方向相反，纖維段(i-1,i)的表面張力Fb,i-1,i由下式確定：

(4)

式中：(xi,yi)、(xi-1,yi-1)分別為珠子i、i-1的坐標；b、j為x和y方向的單位向量；θ為纖維表面張力系數，本文取0.7 kg/s2 [1-3]；ki-1,i為纖維段(i-1,i)的曲率。

黏彈力是表征聚合物非牛頓行為的主要作用力，也是本文的重點研究對象。纖維段(i-1,i)的黏彈力Fve,i-1,i[12]可表示為

(5)

式中：ri與ri-1分別為珠子i與i-1的空間位置；N2,i-1,i為纖維段(i-1,i)所受第二法向應力。

根據張量分析原理，在理想拉伸狀態下，對于某一纖維段，有N2=τt-τn。其中τt和τn分別為軸向偏應力張量和徑向的偏應力張量。在本文研究中，引入4種非牛頓流體本構方程分別計算出4種不同的偏應力張量，進而計算出不同的黏彈力，再結合熔噴纖維拉伸模型預測纖維直徑。這4種非牛頓流體本構方程[13]分別為：

UCM方程

(6)

PTT方程

(7)

Giesekus方程

(8)

Rouse-Zimm方程

(9)

張量分解公式

(10)

(11)

式中：Tf,i-1,i為纖維段(i-1,i)的溫度，K。在實際計算中，利用有限差分方法計算上述非牛頓流體本構方程，得到4個不同的第二法向應力值，并分別代入式(5)計算黏彈力。

為使整個模型封閉，還需纖維的運動控制方程和能量守恒方程，分別為：

(12)

(13)

式中：g為重力加速度，取值為9.8 m/s2；Cb為黑體輻射系數，取值為5.67 W/(m2·K4)[10]；hi-1,i為纖維段(i-1,i)與周圍氣流之間的傳熱系數，取值為0.026 W/(m2·K·s)[10]；ψ為纖維黑度，取值為0.4[10]；Ta,i-1,i為纖維段(i-1,i)周圍的氣流溫度，K；mi-1,i為纖維段(i-1,i)的質量，kg；Si-1,i為參與熱對流的纖維面積，m2；Cg為纖維段(i-1,i)的比熱容，J/(kg·K)，可依據下式[14]進行計算：

Si-1,i=πdi-1,ili-1,i

(14)

(15)

Cg,i-1,i=366.9+2.42Tf,i-1,i

(16)

式中，ρf,i-1,i為纖維段(i-1,i)的密度，kg/m3。以聚丙烯為例，其密度與溫度的經驗方程[14]為

(17)

2 算法描述

熔噴纖維拉伸模型在第1節已經建立，在求解時仍缺少計算氣流力所必需的流場速度和溫度，為此需采用Fluent有限元分析軟件模擬熔噴氣流場(設置流場入口速度為60 m/s，溫度為640 K，模頭溫度為640 K，出口為標準大氣環境)，并導出氣流場三維坐標、速度和溫度的數值解[7-9]。模型的求解方法可概括為：首先依據纖維段的初始物理參數值，并結合式(13)與各作用力的數理方程分別計算出氣流力、表面張力和黏彈力，再依據式(12)求解纖維段的速度與位移，最后依據式(15)求解纖維段的直徑。由于采用了4種不同的非牛頓流體本構方程，黏彈力需分別求解4次，以觀察直徑結果的差異。理論模型的求解過程如下。

步驟1：設置i=2，并賦予纖維段(i-1,i)初始長度為10 mm，初始直徑為400 μm，初始溫度為640 K。依據式(11)與式(15) ～ (17)計算纖維段的初始黏度、質量、比熱容、密度。導入Fluent計算的流場坐標、速度與溫度數值解。

步驟2：計算開始。讀取珠子i-1和i的空間坐標及其周圍的流場溫度和速度。根據式(1) ～ (3)計算纖維段(i-1,i)受到的氣流力；根據式(4)計算其表面張力；從式(6) ～ (9)中選擇一種非牛頓流體本構方程，并根據式(5)計算纖維段(i-1,i)的黏彈力；依據式(13)計算纖維段(i-1,i)的溫度；依據式(15)計算纖維段的直徑；依據式(12)中纖維段受到的合力的大小和方向，計算珠子i-1和i在Δt時間內(可自行設置Δt值，數值越小則輸出數據越多，結果越精確，程序運行越慢；反之則結果精確性較低，程序運行較快)的位移。輸出各作用力的值、珠子空間坐標和纖維段直徑的計算結果。

步驟3：根據步驟2計算的位移，更新珠子i-1和i到新的位置。

步驟4：判斷珠子i與噴絲孔(原點)的距離是否達到給定值(該值可自由設定，越大則纖維越長，計算耗時越久)。如果是，則執行步驟5；反之，返回至步驟2，以計算珠子i-1和i在下一個Δt后的位置。

步驟5：Uyttendaele等[1]提出內外應力一致的“凝固點”用以描述纖維直徑停止變化的邊界條件。在模型中，本文也采用同樣的方法，即判斷纖維段受到的氣流力與重力之和是否等于黏彈力與表面張力之和。如果是，停止計算纖維段的直徑變化；反之，在噴絲孔處加入珠子i+1。賦予系統中新纖維段(i,i+1)相應的初始物理屬性。返回步驟2。

由于需要研究4種不同的非牛頓流體本構方程對計算結果的影響，所以模型需進行多次計算，每次計算采用一種不同的本構方程。在綜合考慮計算時間與計算機性能的基礎上，本文采用400個珠子模擬纖維在流場中的運動和細化程度，程序采用MatLab編程求解。當第400個珠子出現在噴絲孔處時，所有運算停止。

本文研究采用聚丙烯的物理屬性作為原料參數，但該模型可模擬幾乎所有的熔噴原料及工藝過程，僅需進行相應的修改。如更換其他材料，需獲得該材料的黏度、比熱容和密度方程，即修改式(11)、(16)和(17)；還需修改物理參數值至合適的值，包括應力松弛時間、傳熱系數和表面張力系數。若更改其他熔噴工藝條件，可在Fluent模擬氣流場時，設置不同的壓縮氣流速度、氣流溫度、模頭溫度和外部氣流環境等；同時需依據噴絲孔直徑設置纖維的初始直徑，依據模頭溫度設置纖維的初始溫度。

3 結果與討論

3.1 模擬熔噴過程中纖維應力的變化

圖2示出纖維與噴絲孔距離z對模擬獲得的纖維內外應力的影響。

圖2 不同非牛頓流體本構方程得到的纖維應力隨z坐標的變化Fig.2 Fiber stress calculated by different non-Newtonian fluid constitutive equations and its variation with z coordinate

本文采用4種非牛頓流體本構方程計算得到4種不同的內應力。由圖可知，所有曲線的應力值都是在出現1個峰值后，隨z值的增大而逐漸下降。即隨著纖維遠離噴絲孔，應力值逐漸下降至穩定。在纖維應力值達到最大峰值之前，纖維距離噴絲孔較近，由于噴射的2股高速射流在該區域相互融合形成湍流，導致該區域內纖維受到的各項應力不穩定，此時的曲線波動性較大。同時由于此處的高速射流速度最大，所以產生了最大的氣流拉伸應力，纖維受到這種強烈外應力作用需要抵抗拉伸力，故而內應力也達到峰值。隨著纖維遠離噴絲孔，纖維周圍氣流速度逐漸衰減，其受到的外應力逐漸下降，內應力也隨之下降，然而由于熔融纖維的非牛頓流體黏彈特性，黏彈力會有一定的遲滯，這產生了內外應力的應力差，在圖2中體現為外應力曲線與4種內應力曲線在相同z值處的間距。在纖維距離噴絲孔的間距小于約3 cm的區域內，內外應力差較大，這說明該區域內纖維細化程度較大。在間距超過3 cm后，內外應力差逐漸縮小，這說明纖維細化程度逐漸減小。當間距超過5 cm后，內外應力逐漸趨于穩定，這意味著纖維直徑也逐漸停止變化。

3.2 纖維直徑

為與Uyttendaele等的在線測量實驗進行對比以展示模擬結果的準確程度，在模型中采用了與該在線測量實驗相同的工藝條件(具體工藝條件已在第2節說明)，實驗與模擬的纖維直徑變化示于圖3中。為進一步闡明圖3中纖維直徑的變化機制，表1示出采用4種非牛頓流體本構方程計算獲得的指標結果。纖維最大內外應力差是指纖維受到的內外應力差的最大值，是圖2中外應力曲線分別與4種內應力曲線之間的最大差值；纖維平均細化速度是圖3中纖維直徑變化的平均斜率；凝固點位置是指纖維固化時與噴絲孔的間距，在圖3中表示為曲線停止變化時對應的z值；最終纖維直徑是指纖維經過牽伸后的直徑，在圖3中表示為曲線尾端z為10 cm處對應的纖維直徑。

圖3 不同非牛頓流體本構方程得到的纖維直徑隨z坐標的變化Fig.3 Fiber diameter calculated by different non-Newtonian fluid constitutive equations and its variation with z coordinate

非牛頓流體方程名稱最大內外應力差/kPa平均細化速度/(μm·cm-1)凝固點位置/cm最終纖維直徑/μmGiesekus18.72128.006.4574.12PTT12.54126.865.2779.18UCM14.36127.025.3683.42Rouse-Zimm9.21125.807.6274.76

由圖3可知，Uyttendaele實驗結果[1]與模擬結果變化趨勢較為相符。纖維自噴絲孔擠出后，其直徑變化并不穩定，這是由于噴絲孔附近的氣流場形成的湍流所致。隨后纖維直徑快速下降，這是由于纖維較大的內外應力差所致。從表1可明顯發現，纖維細化速度(即曲線在z小于3 cm區域的平均斜率)與內外應力差高度相關，最大內外應力差從大到小排序為：Giesekus > UCM > PTT >Rouse-Zimm，而纖維平均細化速度也有同樣的順序。當纖維與噴絲孔間距z大于3～4 cm時，由于內外應力差縮小，纖維細化速度減緩。間距z超過5～6 cm后，內外應力差穩定導致纖維直徑也逐漸趨于穩定。這期間，當纖維達到凝固點時，纖維直徑停止變化。采用Rouse-Zimm方程計算得到纖維達到凝固點時距離噴絲孔最遠為7.62 cm，而采用PTT方程計算得到纖維達到凝固點時距離噴絲孔最近為5.27 cm。凝固點出現的位置較遠可使纖維有充分的空間拉伸，更容易產生較細的纖維。

最終預測的纖維直徑從大到小排序為：UCM > PTT > Rouse-Zimm > Giesekus。采用Giesekus方程計算得到的纖維直徑更符合Uyttendaele實驗結果。所有預測的最終纖維直徑均大于實驗測量的結果(約60 μm)。這是由于模型中僅考慮了單個纖維段的受力與形變而沒有考慮纖維段之間的相互纏結影響。在Bresee等[17-18]的實驗觀察中，纖維之間明顯的接觸和纏繞現象會導致纖維直徑減小。

4 結 論

本文以熔噴纖維三維拉伸模型為基礎，采用拉格朗日方法對纖維在空氣中的動力學進行分析，同時結合不同的非牛頓流體本構方程，預測纖維直徑的變化情況。雖然模擬結果與實驗結果基本相符，但是不同的非牛頓流體本構方程的應用使計算出的黏彈力有差異，從而導致纖維在拉伸形變時內外應力差異的變化。纖維直徑的模擬結果顯示，纖維的最終直徑不僅受到纖維在氣流場中的內外應力差的影響，還受到凝固點位置的影響，即：纖維的內外應力差越大，纖維細化速度越快；纖維的凝固點距離噴絲孔越遠，則纖維有更加充分的空間進行拉伸，其最終直徑也越細。本文的數學模型研究揭示了熔噴纖維的細化機制，同時非牛頓流體本構方程的引入，也成功定量對比分析了熔噴纖維黏彈力學行為的非牛頓性，為熔噴纖維拉伸模型非牛頓方程的選擇提供了計算依據。

FZXB