陳迪林 張 來 李 捷* 陳海文 錢詩智

(武漢理工大學(xué)能源與動力工程學(xué)院1) 武漢 430061) (老道明大學(xué)微納米技術(shù)學(xué)院2) 諾?？?VA23529)

0 引 言

黏彈性流體屬于時變性非牛頓流體，具有黏性和彈性的雙重特性，且廣泛的存在于生活中的多個工業(yè)部門，如血液、油漆、高分子聚合物溶液、膠體溶液等，深入研究非牛頓流體的流動行為具有重要的工程應(yīng)用價值[1].因此，黏彈性流體流動的分析與建模在工業(yè)領(lǐng)域中具有根本的重要性.黏彈性流體的流變響應(yīng)非常復(fù)雜，包括非線性黏性和彈性效應(yīng)的組合，剪切流中法向應(yīng)力差的存在，弛豫現(xiàn)象和記憶效應(yīng)，沖模膨脹特性等.關(guān)于黏彈性流動的數(shù)值分析，使用一個或多個非線性微分模型的分析處理方法在文獻(xiàn)中都是常見的.無論離散化方法(有限元，有限差分或有限體積)，迭代求解方法或所用的本構(gòu)方程，求解工作中的常見困難是所謂的“高Weissenberg數(shù)問題(HWNP)”. HWNP包含在高Weissenberg(Wi)或高Deborah(De)數(shù)下實(shí)現(xiàn)收斂的困難，其中量綱一的量化的參數(shù)Wi和De用以表征黏彈性流動的流動特征.

當(dāng)下，商業(yè)的軟件分析實(shí)際上僅限于純黏性非牛頓現(xiàn)象描述，具有黏彈性流動分析能力的軟件的開發(fā)和使用仍然幾乎完全在學(xué)術(shù)環(huán)境中進(jìn)行，應(yīng)用于特定領(lǐng)域[2].OpenFOAM(open source field operation and manipulation)中新的求解器，對多弛豫時間、高Wi數(shù)問題、混合網(wǎng)格有著更大的適用性，并且支持并行運(yùn)算與數(shù)據(jù)處理，在解決大網(wǎng)格尺寸的問題上有明顯優(yōu)勢.作為一個經(jīng)過充分測試和廣泛使用的免費(fèi)開源CFD軟件包，更具有C++面向?qū)ο缶幊陶Z言的內(nèi)在優(yōu)勢，OpenFOAM軟件包很明顯是針對具體問題，開發(fā)自定義求解器的最有前途的工具.同時，處理非結(jié)構(gòu)化和動網(wǎng)格以及實(shí)現(xiàn)復(fù)雜數(shù)學(xué)模型的靈活性，是其最吸引人的特征.

因此，于本文中，利用OpenFOAM軟件包中的通用黏彈性流體求解器，采用一種基于分離應(yīng)力張量方法[3]和平衡應(yīng)力張量的概念去處理HWNP問題，且適用于任何微分本構(gòu)方程.通過將其預(yù)測值與來自文獻(xiàn)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較來評估所謂的黏彈性流體求解器ViscoelasticFluidFoam，并用于分析平面4∶1平面收縮高Wi數(shù)流動與腰圓形形功能塊流動的仿真分析.

1 DCPP本構(gòu)模型與基本方程

1.1 數(shù)值模型

黏彈性流體等溫不可壓縮流動的控制方程是質(zhì)量守恒方程(連續(xù)性)和動量守恒方程.

·(u)=0

(1)

(2)

式中：ρ為流體的密度；u為速度矢量；p為壓力；τ為應(yīng)力張量.應(yīng)力張量可分為(溶劑)牛頓應(yīng)力張量τs和彈性聚合物應(yīng)力張量τp(或額外彈性應(yīng)力張量)

τ=τs+τp

(3)

τs=2ηsD

(4)

(5)

式中：ηs為溶劑黏度；D為變形率張量；τpk取決于黏彈性流體的本構(gòu)方程.

(6)

1.2 DCPP本構(gòu)模型

考慮到歐湘萍等[4]提出的XPP模型預(yù)測穩(wěn)態(tài)拉伸流和剪切流時存在解不唯一的問題,盡管有一些改進(jìn)，但當(dāng)包括非消失的第二正應(yīng)力時，他們的新模型在簡單的剪切流中遭受奇異行為，以上都是與實(shí)際高分子熔體的流變特性是不相符的.為了消除模型存在數(shù)學(xué)和流變學(xué)上的這些缺陷，Clemeur等[5]發(fā)展了DCPP模型，其本構(gòu)模型描述為下述的形式.

(7)

式中：T為黏彈性額外應(yīng)力張量；S為方向張量；Λ為拉伸量，用于表征Pom-Pom大分子的行為；I和G為單位張量、剪切模量；ξ為非線性參數(shù)，用于對第二法向應(yīng)力差的控制.張量S和拉伸量滿足下式.

(8)

(9)

式中：D為變形率張量，Δ和分別為上對流導(dǎo)數(shù)和下對流導(dǎo)數(shù)；λ和G為線性參數(shù)的范疇，λ為取向弛豫時間，λs為拉伸機(jī)制的弛豫時間；q依附于Pom-Pom大分子的支鏈數(shù).其他的參數(shù)λs，q，ξ為非線性參數(shù).

1.3 壓力-速度、動量-應(yīng)力耦合及張力黏度

速度-壓力耦合是通過分離方法完成的，其中連續(xù)性方程(2)采用半離散形式，得到的方程組通過解耦方法求解，使用具有欠松弛的迭代算法，如SIMPLE算法[6].

為了保證動量方程解中的動量-應(yīng)力耦合和數(shù)值穩(wěn)定，所采用的策略包括將黏彈性應(yīng)力分解為與D對應(yīng)，基于張力黏度ηT定義的隱式分量，并進(jìn)行顯式校正.

·τ≈·(ηT·u)+·τcorr

(10)

此外，應(yīng)力傳遞模型寫成下面的通用形式：

(11)

該等式描述了應(yīng)力隨速度場在空間中的傳遞，同時在較小的弛豫時間后趨向平衡值τ*，其中τ*是在沒有運(yùn)輸效應(yīng)的情況下實(shí)現(xiàn)的應(yīng)力狀態(tài)，λ是松弛時間.無論使用任何特定的黏彈性本構(gòu)模型，都可以實(shí)現(xiàn)這種形式.

根據(jù)定義導(dǎo)出，與所選模型無關(guān)的ηT表達(dá)式.

ηT=τ*·D-1

(12)

當(dāng)detD=0時，意味著速度梯度是一個奇異張量.在以上的公式中，動量方程為

-p+·τcorr

(13)

1.4 數(shù)值算法

OpenFOAM中，用于求解黏彈性流體流動問題的方法可以歸納為以下四個步驟.

步驟1通過給定的初始速度u、壓力p、應(yīng)力τ，對壓力梯度與應(yīng)力差進(jìn)行顯式計(jì)算.隨后隱式求解動量方程，得到一個新的速度預(yù)估值u*.

步驟2隨著u*的計(jì)算，得到新的壓力預(yù)估值p*，隨后進(jìn)行速度修正，得到新的滿足連續(xù)性方程的新速度場u**.選擇采用較好選擇的SIMPLE或者PISO算法進(jìn)行計(jì)算.其中，瞬態(tài)流動PISO算法更具有適用性.

步驟3隨著修正后的速度場u**的計(jì)算，代入指定的本構(gòu)方程中，求解得到應(yīng)力張量場τ*.

步驟4瞬態(tài)流動中，前三個步驟每個時間步不斷重復(fù)，能夠生成得到更加精確的解.迭代過程中，u**，p*和τ*不斷地獲得更新.

動量方程(2)在OpenFOAM中的形式如下，具體的代碼含義不做贅述，可以參考文獻(xiàn)[7].

2 結(jié)果與討論

2.1 收縮流幾何和流動條件

通過與Clemeur等[8]的數(shù)據(jù)進(jìn)行比較來測試DCPP模型的實(shí)施.圖1所示的平面突然收縮， 2H=4.10 mm和2h=1 mm，收縮比為4.1∶1. 網(wǎng)格單元數(shù)為23 800，上游和下游流道長度分別為40，80 mm，以使其充分發(fā)展.使下游壁面剪切速率12.4 s-1，得到Weissenberg數(shù)為50.0.基于第一法向應(yīng)力差N1(N1=τXX-τYY)定義主應(yīng)力差PSD，用于參數(shù)比較.

圖1 收縮流網(wǎng)格

(14)

表1中列出了DCPP模型參數(shù)的對應(yīng)值.

表1 DCPP模型參數(shù)

2.2 DCPP的結(jié)果驗(yàn)證

運(yùn)用Paraview后處理工具，展示了圖2的等色條紋圖.定值描述其變化趨勢，速度值U與主應(yīng)力差PSD的預(yù)測值見圖3.結(jié)果表明沿著中心線y=0并且沿著于y=0.2,y=0.4,y=0.49獲得了良好的收斂，除了收縮的入口有較小的偏差，但與Favero等[9]研究數(shù)據(jù)的誤差不超過5%.主應(yīng)力差與速度曲線都在一個較狹窄的范圍內(nèi)，沿著y=0.49的預(yù)測曲線不是特別理想，幾何拐點(diǎn)的附近是否平滑，可能是導(dǎo)致預(yù)測值與實(shí)際差值的原因.PSD預(yù)測值與實(shí)驗(yàn)值趨勢近似，但y=0線的最大值偏小，誤差可能是由于收縮口前后網(wǎng)格不夠細(xì)化所導(dǎo)致.文獻(xiàn)[10]中網(wǎng)格在收縮處有半徑為3 μm的圓角，而本文中忽略了圓角帶來的影響.總體而言，仿真所得結(jié)果數(shù)據(jù)是令人滿意并良好吻合的.

圖2 速度和主應(yīng)力差等色條紋圖

圖3 收縮前后的速度值和特征主應(yīng)力差

2.3 腰圓形功能塊收縮流動的數(shù)值模擬

求解模擬聚合物在復(fù)雜流動幾何中的流動，關(guān)鍵在于選取合適的黏彈性流體本構(gòu)模型.DCPP是基于Pom-Pom的優(yōu)化模型，能夠很好的規(guī)避剪切流動中的奇異行為，預(yù)測聚合物溶液流動的拉伸和剪切作用.本節(jié)選用DCPP本構(gòu)模型，針對注塑成型過程中腰圓形功能塊收縮流動，對其流函數(shù)、速度、壓力、主應(yīng)力差進(jìn)行了仿真，并與不存在嵌件的流場數(shù)據(jù)進(jìn)行對比分析.

模型收縮比4.1∶1.腰圓形處L=1 mm，總長2 mm，圓弧半徑r=0.5 mm，見圖4.在腰圓形壁面及前后流場進(jìn)行局部網(wǎng)格加密，網(wǎng)格數(shù)30 533.Wi數(shù)50，模型參數(shù)按照2.2選取.選用流量邊界條件，固壁邊界，定值壓力出口，其他設(shè)為零梯度.

圖4 腰圓形功能塊幾何

由圖4可知，為收縮流道不存在功能塊與存在功能塊流場的數(shù)據(jù)對比圖.圖5為收縮口前后的流線圖.很明顯，在功能塊的作用下，流體流線在中心線兩側(cè)發(fā)生分離；在功能塊的上平面與下平面區(qū)域，流線近似平行于軸(X)方向流動，進(jìn)入收縮口之前都緊緊的貼合于功能塊壁面流動，兩邊繞行之后在功能塊后側(cè)進(jìn)行匯合；圖中可見，功能塊的圓弧過渡段，由于不存在幾何尖角，流線在流過功能塊之前方向偏轉(zhuǎn).近二次渦流橫截面上流線更加緊密，是由于功能塊導(dǎo)致的收縮口前的又一次收縮，分子鏈在這一區(qū)域也會有大的強(qiáng)拉伸剪切行為.流體流過功能塊后側(cè)時發(fā)生急劇的方向偏轉(zhuǎn)，是由于橫截面積的急劇收縮及二次渦流的影響，將使流線更加集中，流速逐漸升高，這樣的流動變化與實(shí)際工程中聚合物成型的規(guī)律是極其吻合的.

圖5 流函數(shù)

圖6為y=0與y=0.4上速度U的對比圖.

圖6 速度趨勢

y=0為中軸線上的速度，在功能塊的作用下，無論是近壁面(y=0.4)還是中心線(y=0)上，都經(jīng)歷了一個巨大的速度損失；相較而言，在中心線流道上，速度U在功能塊與收縮口之間的管段速度下降較少，更大的速度損失是在收縮口之后的管段；在近壁面處，速度U在收縮口前后的速度U幾乎是同比下降，趨勢近似；中心線趨勢變化的原因：流體繞過功能塊頂部后又進(jìn)入了一個比較狹窄的通道內(nèi)，大分子鏈有一個剪切到拉伸的過程，且塊后局部流體受到上下分開的流道內(nèi)的流體擠壓，速度進(jìn)一步降低，導(dǎo)致了收縮口后中心線上的速度損失更大；近壁面處，由于流體都是通過寬的流道隨著流線而過，流動中受到的阻力更小，因此速度損失較小，幾乎是同比下降.

圖7a)為y=0和y=0.4在收縮口前后的壓力降(圖中y1為收縮流道預(yù)測值，y*為功能塊收縮流道預(yù)測值).由結(jié)果可得在收縮口前后遭遇了大的壓力損失，無功能塊的近壁面壓力曲線與中心流道壓力曲線重合，存在功能塊的情況下的預(yù)測值與收縮流道的趨勢一樣，也遇到了重合的預(yù)測曲線；可見對于壓力p的預(yù)測值皆獲得了良好的收斂，與實(shí)際情況下的壓力變化吻合.

圖7 腰圓形功能塊收縮流動的壓力與主應(yīng)力差對比

圖7b)為y=0與y=0.4在收縮口前后的主應(yīng)力差PSD的預(yù)測值.預(yù)測的數(shù)值在中心線流道y=0上獲得了良好的收斂，靠近壁面的預(yù)測值則不是那么理想，在收縮口下游拐角后流道有大的偏差與波動；功能塊的出現(xiàn)導(dǎo)致了主應(yīng)力差有部分降低.根據(jù)預(yù)測值的表征，具有彎角的功能塊使得收縮口后半段的應(yīng)力差波動越加穩(wěn)定，并沒有y1=0.4后半段波動那么明顯.功能塊后段流域,第一法向應(yīng)力差與拉伸量同比變化，且取值較大.這是由于在功能塊后側(cè),收縮流道中的二次回流，將增大分子鏈拉伸量，表征為聚合物分子彈性的增強(qiáng)，漩渦產(chǎn)生了一個強(qiáng)剪切作用的區(qū)域；在流過功能塊與壁面的小收縮區(qū)域后，分子鏈部分彈性增強(qiáng)，當(dāng)還未完全松弛的情況下，又經(jīng)過回流區(qū)域的強(qiáng)拉伸區(qū)域及收縮拐角的剪切拉伸作用，聚合物彈性得到進(jìn)一步增強(qiáng).因此，主應(yīng)力差的最大值出現(xiàn)在這一區(qū)域；另外一方面，在靠近壁面處的第一法向應(yīng)力差的最大值較純收縮流動的最大值有所滯后，主要是由于通道截面的多次急劇變化以及回流的作用，使得分子鏈向平衡態(tài)舒緩的過程變慢，剪切作用的疊加與拉伸舒緩的雙重效果，共同形成了極值的滯后；在收縮口下游區(qū)域，流道得到充分發(fā)展，流場中剪切拉伸作用逐漸弱化，分子鏈進(jìn)一步松弛，靠近壁面流體，由于受到拐角的擾動更大，向平衡態(tài)舒緩的速度更慢，而中心線附近流體則舒緩較快，二者都表征為第一法向應(yīng)力差的逐漸減小.

3 結(jié) 束 語

本文中網(wǎng)格收斂與流動趨勢、數(shù)值預(yù)測方面表現(xiàn)良好，確認(rèn)所采用的方法是成功的；同時針對聚合物加工過程中常見的腰圓形功能塊收縮流動，對流函數(shù)、U、PSD等聚合物濃厚體系表征量進(jìn)行了預(yù)測分析，并與不存在功能塊的流場數(shù)據(jù)進(jìn)行對比，揭示其內(nèi)部規(guī)律與機(jī)理，能夠更好的服務(wù)于聚合物成型、模具優(yōu)化、缺陷防范的相關(guān)研究；本文選用開源OpenFOAM工具，其所具有的優(yōu)勢與內(nèi)在屬性，能夠高靈活性的求解黏彈性流體流動的相關(guān)問題，為黏彈性流動的開發(fā)應(yīng)用提供巨大潛力.