eXtended Pom-Pom 黏彈性流體的改進光滑粒子動力學模擬*

許曉陽 周亞麗 余鵬

1)(西安科技大學計算機科學與技術(shù)學院,西安 710054)

2)(南方科技大學力學與航空航天工程系,深圳 518055)

黏彈性流體廣泛存在于自然界和工業(yè)生產(chǎn)中,對其復雜流變特性的研究具有重要的學術(shù)價值和應(yīng)用意義.本文提出一種改進的光滑粒子動力學方法,對基于eXtended Pom-Pom 模型的黏彈性流動進行了數(shù)值模擬.為了提高計算精度,采用一種不含核導數(shù)計算的核梯度修正離散格式.為了防止粒子穿透固壁,提出一種增強型的邊界處理技術(shù).為了消除張力不穩(wěn)定性,將人工應(yīng)力耦合到動量守恒方程中.運用改進光滑粒子動力學方法數(shù)值模擬了基于eXtended Pom-Pom 模型的黏彈性Poiseuille 流和黏彈性液滴撞擊固壁問題,通過與解析解或有限差分方法解的比較以及對數(shù)值收斂性的評價,驗證了改進光滑粒子動力學方法的有效性和優(yōu)勢,并在此基礎(chǔ)上,深入分析了Reyonlds 數(shù)、Weissenberg 數(shù)、溶劑黏度比、各向異性參數(shù)、松弛時間比和分子鏈臂數(shù)等流變參數(shù)對流動過程的影響.

1 引言

黏彈性流體作為非牛頓流體力學領(lǐng)域里的一類重要研究對象,廣泛存在于自然界和工業(yè)生產(chǎn)過程中,如石油輸送和聚合物加工成型等.為了提高產(chǎn)品質(zhì)量和生產(chǎn)效率,對黏彈性流體復雜流變特性的相關(guān)研究已成為工業(yè)生產(chǎn)過程中一個重要的環(huán)節(jié).目前,基于網(wǎng)格的數(shù)值方法,如有限差分法(finite difference method,FDM)、有限體積法和有限元等,已被用于該類問題的數(shù)值模擬.Viezel等[1]使用FDM 模擬了Oldroyd-B 液滴碰撞問題.Li等[2]通過流體體積法模擬了黏彈性液滴在低Weber數(shù)下的流動.盡管基于網(wǎng)格的數(shù)值方法在模擬黏彈性流體方面取得了成功,但其自身依賴于網(wǎng)格生成的質(zhì)量.此外,對于含自由面的黏彈性流動模擬,需借助額外的界面追蹤技術(shù),實施過程相對繁瑣.

近年來,無網(wǎng)格粒子法[3]引起了廣大學者的關(guān)注以及重視.光滑粒子動力學(smoothed particle hydrodynamics,SPH)方法是一種典型的無網(wǎng)格粒子法,1977 年由Gingold,Monaghan[4]和Lucy[5]提出.與基于網(wǎng)格的數(shù)值方法相比,SPH 方法具有無網(wǎng)格特性、自適應(yīng)特性、Lagrangian 特性和易于編程等優(yōu)勢,因此該方法已被成功用于黏性流[6]、動畫[7]、傳熱[8]、多相流[9]、爆炸[10]和流固耦合[11]等領(lǐng)域中.然而,傳統(tǒng)SPH 方法也存在一些缺點,如計算精度低、穩(wěn)定性差和邊界不易處理等.針對這些問題,近年來很多學者已提出了多種不同的改進方案.例如,為了提高計算精度,Liu等[12]用積分再生核思想發(fā)展了再生核粒子方法(reproducing kernel particle method,RKPM);Liu等[13]推導了具有較高計算精度的有限粒子法(finite particle method,FPM);Fang等[14]提出了模擬黏性不可壓縮流動的高精度SPH 方法.為了改善SPH 方法的數(shù)值穩(wěn)定性,Yang等[15]提出了一種新的核函數(shù),克服了黏性液滴形成過程中的張力不穩(wěn)定;Antuono等[16]在密度方程中引入密度耗散項,以解決流動過程中出現(xiàn)的壓力不穩(wěn)定;Lyu等[17]發(fā)展了粒子遷移技術(shù),保證流動過程中流體粒子的規(guī)則分布.為了施加邊界條件,Monaghan等[18]提出了只施加一層粒子在邊界的邊界排斥力法;Morris等[19]發(fā)展了布置多層粒子在邊界的虛擬粒子法;Liu等[20]結(jié)合排斥力法和虛擬粒子法的優(yōu)勢,建立了兩者相耦合的動力學邊界處理技術(shù).雖然這些改進方法已被成功提出并得到一定的應(yīng)用,但它們在黏彈性流體力學領(lǐng)域中的應(yīng)用并不多見.

對于黏彈性流體的SPH 模擬,Fang等[21]將SPH 方法推廣到含自由面的Oldroyd-B 黏彈性液滴撞擊固壁問題的模擬中.Hashemi等[22]研究了兩個相互作用的固體顆粒在Oldroyd-B 剪切流場下的遷移運動,分析了Deborah 數(shù)對固體顆粒運動軌跡的影響.Xu等[23]提出了改進SPH 方法模擬基于Oldroyd-B 模型的黏彈性液滴撞擊固壁和擠出脹大問題.Ozgen等[24]利用SPH 方法捕捉了非連續(xù)剪切增稠流體特有的流動現(xiàn)象.Vahabi 和Kamkari[25]基于SPH 方法數(shù)值模擬了氣泡在Giesekus 黏彈性溶液中的上升和變形.King 和Lind[26]將對數(shù)構(gòu)象公式耦合到SPH 方法模擬了高Weissenberg 數(shù)下的黏彈性圓柱繞流問題.

eXtended Pom-Pom(XPP)模型由Verbeeten等[27]在Pom-Pom 模型的基礎(chǔ)上基于支鏈聚合物的分子理論推導得到.與傳統(tǒng)的簡單模型(例如Oldroyd-B,上隨體Maxwell 模型)相比,XPP 模型能夠更合理地描述聚合物溶液的剪切和拉伸行為,且在一定程度上克服了應(yīng)力奇點問題.此外大量實驗[27]表明,XPP 模型能夠合理地表征黏彈性流體的剪切稀化.根據(jù)文獻調(diào)研結(jié)果,目前運用SPH 方法模擬XPP 黏彈性流體的文獻報道并不多見,因此本文重點圍繞基于XPP 模型的黏彈性流動問題展開研究.

本文首先對傳統(tǒng)SPH 方法進行改進.特別地,采用一種不含核導數(shù)計算的核梯度修正離散格式以提高計算精度.為了防止粒子穿透固壁,提出一種增強型的邊界處理技術(shù).為了消除張力不穩(wěn)定性,將人工應(yīng)力耦合到動量守恒方程中.隨后,運用改進SPH 方法數(shù)值模擬了基于XPP 模型的黏彈性Poiseuille 流和黏彈性液滴撞擊固壁問題,通過與解析解或其他方法解的比較和對數(shù)值收斂性的評價驗證了改進SPH 方法的有效性和優(yōu)勢,并在此基礎(chǔ)上,深入分析了Reyonlds 數(shù)、Weissenberg數(shù)、溶劑黏度比、各向異性參數(shù)、松弛時間比和分子鏈臂數(shù)等流變參數(shù)對流動過程的影響.

2 控制方程

在Lagranigian 坐標系下,二維等溫、黏彈性流體的控制方程可寫為

其中,ηs是牛頓溶劑黏度;dαβ是形變率張量,

2.1 XPP 本構(gòu)方程

為了封閉控制方程,需附加一個與彈性應(yīng)力相關(guān)的本構(gòu)方程.本文考慮XPP 模型,其本構(gòu)方程為

式中,

λ1和λ2分別表示聚合物分子鏈的取向松弛和拉伸松弛時間;γλ2/λ1表示拉伸松弛與取向松弛時間之比;G0表示線性松弛模量;ηpG0λ1表示聚合物黏度,ηηs+ηp表示流體總黏度;α為控制各向異性拖曳力的流變參數(shù);Q0反比于分子鏈臂數(shù)Q,即Q0=2/Q;分子鏈拉伸量λ的表達式為

由(6)式—(9)式,可得

顯然,當α0和f(λ,τ)1時,XPP 本構(gòu)方程將退化為Oldroyd-B 本構(gòu)方程.

2.2 狀態(tài)方程

在SPH 框架內(nèi),有兩種常用的求解控制方程的方法:一種是不可壓SPH 方法[28],其中流體壓力通過Possion 方程隱式計算;另一種是弱可壓SPH 方法[23],其中壓力根據(jù)狀態(tài)方程顯式計算.本文采用后者,使用的狀態(tài)方程為

式中,c是聲速,ρ0是初始流體密度.在弱可壓SPH方法中,馬赫數(shù)應(yīng)小于0.1.基于此,聲速c應(yīng)取為不低于流動特征速度值的10 倍,以保證弱可壓流體的流動行為足夠接近不可壓流體的行為.聲速c越大,模擬所用時間步長 Δt越小.

3 SPH 方法

3.1 SPH 離散

對流體控制方程可進行多種不同的SPH 離散.對質(zhì)量、動量守恒方程和XPP 本構(gòu)方程采用的SPH 離散為

其中,

i和j表示粒子編號,m表示粒子質(zhì)量,WijW(|xi-xj|,h)是核函數(shù),h是核函數(shù)影響區(qū)域的光滑長度.

傳統(tǒng)SPH 方法的計算精度低,限制了其向更復雜流動問題的應(yīng)用.為了提高計算精度,許多學者已提出多種修正的SPH 算法,如RKPM[12]和FPM[13]等.對SPH 離散(12)式—(16)式中的核函數(shù)梯度進行修正,即利用如下不含核導數(shù)計算的修正矩陣[29]:

將原核函數(shù)梯度修正為

利用新的核梯度(18)式進行計算,可以提高SPH 方法的計算精度.另外,值得注意的是,本文核梯度修正矩陣不含核函數(shù)一階導數(shù)的計算,因此對核函數(shù)的適用性更高,進行相應(yīng)模擬的穩(wěn)定性也更好.關(guān)于改進SPH 離散的更多內(nèi)容,可參閱文獻[29].

3.2 邊界處理

邊界處理是數(shù)值模擬的一個重要部分,其處理是否得當關(guān)乎SPH 模擬的穩(wěn)定性和精度.目前,較為廣泛使用的邊界方法有排斥力法[18]和虛擬粒子法[19]等.排斥力法只施加一層粒子在邊界,程序設(shè)計簡單,但守恒性較差.虛擬粒子法布置多層粒子在邊界,解決了鄰近邊界處粒子被邊界截斷的問題,但虛擬粒子的速度、壓力等各物理量需通過構(gòu)造偽粒子進行Shepard 插值得到.對于復雜邊界,偽粒子的構(gòu)造較為繁瑣.

本文提出一種由邊界粒子和虛擬粒子兩類粒子組成的增強型邊界處理技術(shù).首先,在邊界上以粒子初始間距 Δx布置一層邊界粒子.與排斥力法[18]不同的是,邊界粒子不再對鄰近邊界處的流體粒子施加排斥力.相反,邊界粒子參與控制方程中速度、壓力和彈性應(yīng)力的計算,以有效防止粒子穿透邊界.計算過程中,邊界粒子的密度和位置保持不變.速度設(shè)置為0,以施加無滑移邊界條件.壓力和彈性應(yīng)力通過對相鄰流體粒子進行Shepard 插值得到:

式中,A表示壓力或彈性應(yīng)力,i為邊界粒子,j為與粒子i相鄰的流體粒子.

其次,為解決鄰近邊界處粒子被邊界截斷的問題,在邊界外以粒子初始間距 Δx規(guī)則布置三層虛擬粒子.與邊界粒子類似,虛擬粒子的密度和位置在計算過程中也保持不變.但與虛擬粒子法[19]不同的是,虛擬粒子的各物理量不再通過構(gòu)造偽粒子進行Shepard 插值得到.本文每個虛擬粒子均與邊界法向上的一個邊界粒子相對應(yīng).為滿足壓力和彈性應(yīng)力的紐曼邊界條件,虛擬粒子的壓力和彈性應(yīng)力設(shè)置與邊界法向上對應(yīng)的邊界粒子值相同.另外,所有虛擬粒子的速度設(shè)置為0,以滿足無滑移邊界條件.

相較于傳統(tǒng)的排斥力法[18]和虛擬粒子法[19],本文提出的增強型邊界處理技術(shù)不僅程序設(shè)計簡單,而且避免了構(gòu)造偽粒子進行Shepard 插值的繁瑣操作.數(shù)值試驗表明,本文提出的增強型邊界處理技術(shù)能夠獲得精確穩(wěn)定的計算結(jié)果;即便對于沖擊力較大的黏彈性液滴撞擊固壁問題,也能夠有效地防止流體粒子穿透邊界.

3.3 人工應(yīng)力

當用SPH 方法模擬固體彈性變形問題時,會產(chǎn)生張力不穩(wěn)定性.在進行含自由面的黏彈性流動問題模擬時,由于彈性力的作用,也會產(chǎn)生張力不穩(wěn)定性,即粒子在運動過程中形成團塊并最終導致模擬中斷.對于該問題,目前已有一些解決方法,如采用新核函數(shù)法[15]、粒子位移修正法[30]等.本文將Monaghan[31]和Gray等[32]提出的人工應(yīng)力耦合到動量守恒方程中,以消除張力不穩(wěn)定性.人工應(yīng)力的基本思想是在一對相鄰的粒子之間引入一個小的短程排斥力,以防止兩個粒子靠的太近而形成團塊,其表達式為

其中,nW(0,h)/W(Δx,h)為指數(shù)型因子,fijWij/W(Δx,h)為人工應(yīng)力系數(shù),Rαβ為人工應(yīng)力分量,Δx為粒子初始間距.

二維問題中,人工應(yīng)力張量Rαβ的分量可通過以下方式轉(zhuǎn)換得到:

其中,R′xx和R′yy分別是沿x′軸和y′軸人工應(yīng)力張量的對角分量.坐標系(x′,y′)相對于(x,y)旋轉(zhuǎn)θ角度,即

R′xx的表達式為

其中ε是一個介于0 和1 之間的人工應(yīng)力參數(shù),R′yy也是類似的形式.對角分量σ′xx和σ′yy可通過以下關(guān)系式計算:

于是,通過耦合人工應(yīng)力,動量方程(13)的SPH離散形式修正為

3.4 時間積分

為了求解控制方程,采用蛙跳公式進行時間積分.令Xi表示未知變量(ρi,ui,τi)的向量,Bi表示所對應(yīng)方程的右端向量,那么蛙跳公式可概括如下:

在第一個時間步長t0結(jié)束后,Xi向前推進半個時間步,ri向前推進一個時間步:

為了實現(xiàn)每個時間步的一致性,在每個時間步的開端,各粒子的速度向前推進半個時間步,從而保證和位移的一致性:

在隨后時間步的末端,粒子的密度、速度和位移按照標準的蛙跳公式進行推進:

此外,為了保持數(shù)值穩(wěn)定性,時間步長 Δt必須滿足以下條件:

其中,Fi表示每單位質(zhì)量的力,υ0η/ρ表示流體的運動黏度.

4 數(shù)值算例

4.1 黏彈性Poiseuille流

目前,尚未見到基于XPP 模型對黏彈性Poiseuille 流進行SPH 模擬研究的相關(guān)文獻報道,因此首先選取該問題作為數(shù)值模擬的第一個算例,并分析XPP 模型中各流變參數(shù)對黏彈性流動過程的影響.在模擬過程中,為了更好地描述黏彈性流動,引入無量綱Reynolds數(shù)

無量綱Weissenberg數(shù)

以及無量綱溶劑黏度比

這里V和L分別表示流體流動的特征速度和特征長度.Wi比較了彈性力與黏性力,通常由流體的應(yīng)力松弛時間與具體的加工時間的關(guān)系給出.

圖1 給出了黏彈性Poiseuille 流的幾何模型和初始狀態(tài).當t0時,流體位于兩塊水平固定且間距為L的無限長平板之間,流體初始速度為0.當t ＞0時,流體受到水平外力F的作用開始流動.模擬中,沿x軸方向采用周期邊界條件,故計算區(qū)域可取為Lx×Ly0.5×1.水平外力F1.0,粒子密度ρ1,各流變參數(shù)取值如下:Re=2,Wi=1,β=0.1,α=0.01,γ=0.9,Q=4,其中特征長度L=1,特征速度V0.25 m·s-1.粒子初始間距為 Δx0.02 m,對應(yīng)于1274 個流體粒子、64 個邊界粒子和192 個虛擬粒子.核函數(shù)采用5 次樣條函數(shù),光滑長度h0.9Δx.聲速c10V,時間步長 Δt1.0×10-5.對于聲速c,還選取較大的 20V和40V進行數(shù)值實驗,所得模擬結(jié)果和 10V的模擬結(jié)果幾乎相同,表明聲速c取值為特征速度的10 倍,已能夠保證弱可壓流體的流動行為充分接近不可壓流體的行為.

圖1 黏彈性Poiseuille 流的幾何區(qū)域Fig.1.Geometric region of viscoelastic Poiseuille flow.

圖2(a)給出了利用本文改進SPH 方法模擬基于XPP 模型的黏彈性Poiseuille 流關(guān)于y軸的速度分布.可以看出,黏彈性Poiseuille 流的速度關(guān)于y=0 對稱;且由于受到彈性力作用的影響,速度有明顯的過沖現(xiàn)象,即速度在t=0.8 時過沖到最大,隨后在t=2.1 時降低至最小,最終在t=8時達到穩(wěn)態(tài).空間點1—3(如圖1 所示)的速度隨時間的變化情況如圖2(b)所示.很明顯,3 個空間點均展現(xiàn)出了速度過沖現(xiàn)象;且空間點離邊界越近,速度值越小.此外,由于空間點1 和3 關(guān)于中心點2 對稱,因此它們的速度值相同.從圖2(c)可知,空間點2 的彈性剪切應(yīng)力一直為0,而空間點1和3 的彈性剪切應(yīng)力大小相等,方向相反.顯然,空間點離邊界越近,彈性剪切應(yīng)力值越大.

圖2 基于XPP 模型的黏彈性Poiseuille 流的SPH 模擬(Re=2,Wi=1,β=0.1,α=0.01,γ=0.9,Q=4)(a)速度分布圖;(b)空間點1—3 速度u 隨時間的變化;(c)空間點1—3 彈性剪切應(yīng)力τxy 隨時間的變化Fig.2.SPH simulation of viscoelastic Poiseuille flow based on XPP model(Re=2,Wi=1,β=0.1,α=0.01,γ=0.9,Q=4):(a)Velocity profile;(b)time change of velocityu at points 1 to 3;(c)time change of elastic shear stressτxyat points 1 to 3.

4.1.1 改進SPH 方法的有效性和優(yōu)勢

目前,基于XPP 模型的黏彈性Poiseulle 流尚沒有瞬態(tài)解析解.但當f(λ,τi)1和α0時,XPP 本構(gòu)方程將退化為Oldroyd-B 本構(gòu)方程,而基于Oldroyd-B 模型的黏彈性Poiseuille 流具有瞬態(tài)解析解[33].因此,為了驗證本文改進SPH 方法模擬XPP 黏彈性流體的準確性和可靠性,增加了f(λ,τi)1和α0的數(shù)值試驗,而其他參數(shù)保持不變.將得到的u數(shù)值解與基于Oldroyd-B 模型的u解析解進行比較,如圖3 所示,SPH 結(jié)果與解析解吻合很好,從而驗證了本文提出的增強型邊界處理技術(shù)是合理的,以及本文改進的SPH 方法模擬XPP 黏彈性流體是有效的.

圖3 本文改進SPH 數(shù)值解與傳統(tǒng)SPH 數(shù)值解和解析解的比較Fig.3.Comparison of improved SPH solutions with original SPH solutions and analytical solutions.

為了展現(xiàn)本文改進SPH 方法的優(yōu)勢,這里利用傳統(tǒng)SPH 方法對f(λ,τi)1和α0的數(shù)值實驗進行了測試,計算得到的u數(shù)值解如圖3 所示.同時,為了更好地展示不同曲線標志的差異性,還增加了特定區(qū)域的局部放大圖.可以看出,相比于傳統(tǒng)SPH 方法,利用本文改進SPH 方法得到的u數(shù)值解更接近解析解.另外,為了定量地體現(xiàn)本文改進SPH 方法的優(yōu)勢,進一步引入L-2 范數(shù)誤差:

其中,R代表SPH 數(shù)值解或解析解.圖4 給出了利用本文改進SPH 方法和傳統(tǒng)SPH 方法得到的u數(shù)值解與解析解的L-2 范數(shù)誤差隨時間變化的比較.很明顯,相比于傳統(tǒng)SPH 方法,利用本文改進SPH 方法得到的u數(shù)值解與解析解的L-2 范數(shù)誤差更小,這表明了本文改進SPH 方法相較于傳統(tǒng)SPH 方法具有更高的計算精度.

圖4 利用傳統(tǒng)SPH 和改進SPH 方法得到的u 數(shù)值解與解析解的L-2 范數(shù)誤差隨時間變化的比較Fig.4.Comparison of the time change of L-2 norm error obtained based on original SPH solutions and improved SPH solutions.

另外,為了評價本文SPH 方法模擬XPP 黏彈性流體的數(shù)值收斂性,以圖2 為例(Re=2,Wi=1,β=0.1,α=0.01,γ=0.9,Q=4),本文特別增加了粒子初始間距分別為 Δx0.025,0.0125和0.01 的三個數(shù)值實驗,而其他參數(shù)均保持不變.所有數(shù)值試驗的時間步長均設(shè)置為 Δt1.0×10-5,以保證數(shù)值穩(wěn)定性.圖5(a)和圖5(b)分別給出了利用不同粒子初始間距 Δx得到的點2 處的速度和點1 處的彈性剪切應(yīng)力隨時間的變化情況.可以看出,利用不同粒子初始間距 Δx得到的點2 處的速度和點1 處的彈性剪切應(yīng)力隨時間的變化均相同,從而表明了本文SPH 方法具有良好的數(shù)值收斂性.

圖5 利用不同粒子初始間距Δx 得到的SPH 數(shù)值解(a)空間點2 處速度u 隨時間的變化;(b)空間點1 處彈性剪切應(yīng)力τxy 隨時間的變化Fig.5.SPH solutions obtained by different initial particle spacings Δx :(a)Time change of velocityu at point 2;(b)time change of elastic shear stressτxyat point 1.

4.1.2 XPP 流變參數(shù)對流動過程的影響

下面討論XPP 本構(gòu)模型中各流變參數(shù)對流動過程的影響.由于空間點離邊界越近,速度值越小,彈性剪切應(yīng)力值越大,故討論各流變參數(shù)對速度的影響時,只考慮中心點2;當討論各流變參數(shù)對彈性剪切應(yīng)力的影響時,只考慮中心點1.

1)Re的影響

Re對速度u和彈性剪切應(yīng)力τxy的影響如圖6(a)和圖7(a)所示.由圖6(a)可見,Re越大,速度過沖達到的最值越大,達到最值所需的時間越長.對于Re=1,2,3,4 和5,速度過沖達到的最值分別為0.385,0.575,0.710,0.890 和1.000,達到最值所需的時間分別為0.510,0.575,1.105,1.255和1.450.另外,能觀察到速度過沖現(xiàn)象隨著Re的增大而減弱.當Re=2 時,有2 次明顯的速度過沖現(xiàn)象;但當Re增大到4 以上時,卻只觀察到了1 次,這是因為Re越大,流體慣性力越大,彈性力相對變小.Re對速度的穩(wěn)態(tài)值也有較大的影響.Re越大,速度達到的穩(wěn)態(tài)值越大.Re=1 時的速度穩(wěn)態(tài)值為0.115,明顯低于Re=4 和5 時的速度穩(wěn)態(tài)值(分別為0.625 和0.920).對于彈性剪切應(yīng)力:Re越大,彈性剪切應(yīng)力過沖現(xiàn)象越弱,但過沖達到的最值越大,達到最值所需的時間越長.此外,Re對彈性剪切應(yīng)力穩(wěn)態(tài)值的影響不大,表明Re并非影響彈性力的主要因素.

2)Wi的影響

Wi是黏彈性流動的一個重要參數(shù).圖6(b)和圖7(b)分別描述了Wi對速度u和彈性剪切應(yīng)力τxy的影響.由圖6(b)可見,Wi越大,速度過沖達到的最值越大,達到最值所需的時間越長.對于Wi=0.1,0.2,0.5,1.0,2.0,3.0,4.0 和5.0,速度過沖達到的最值分別為0.220,0.310,0.425,0.530,0.715,0.815,0.900 和0.985,達到最值所需的時間分別為0.355,0.410,0.520,0.750,0.955,1.115,1.425和1.655.另外,速度過沖現(xiàn)象隨著Wi的增大先增強后減弱.從圖6(b)可以看出,Wi=1 時速度過沖現(xiàn)象最為顯著,發(fā)生有2 次.低于或高于該值,速度過沖現(xiàn)象均減弱.當Wi降到0.2 以下時,僅觀察到1 次,這是因為Wi越小,彈性力越小.當Wi增到4 以上時,也僅觀察到1 次,這是因為XPP模型可合理地描述流體的剪切變稀行為;當Wi較大時,流體發(fā)生了剪切變稀.Wi對速度穩(wěn)態(tài)值也有顯著的影響,即Wi越大,速度達到的穩(wěn)態(tài)值越大.由圖7(b)可見,由于剪切變稀,隨著Wi的增大,彈性剪切應(yīng)力的穩(wěn)態(tài)值反而變小,但過沖達到的最值先增大后減小,達到最值所需的時間越長.

圖6 不同流變參數(shù)下空間點2 處速度u 隨時間的變化情況(a)Re;(b)Wi;(c)β;(d)α;(e)γ;(f)QFig.6.Time change of velocityu at point 2 under different rheological parameters:(a)Re;(b)Wi;(c)β;(d)α;(e)γ;(f)Q.

圖7 不同流變參數(shù)下空間點1 處彈性剪切應(yīng)力τxy 隨時間的變化情況(a)Re;(b)Wi;(c)β;(d)α;(e)γ;(f)QFig.7.Time change of elastic shear stressτxy at point 1 under different rheological parameters:(a)Re;(b)Wi;(c)β;(d)α;(e)γ;(f)Q.

3)β的影響

β表示黏彈性流動中牛頓溶劑黏度與總黏度的比值.圖6(c)和圖7(c)分別給出了β對速度u和彈性剪切應(yīng)力τxy的影響.由圖6(c)可見,隨著β的增大,速度過沖現(xiàn)象逐漸變?nèi)?這是因為β越大,牛頓溶劑黏度越高,流體越接近于牛頓流體.對于β=0.1,0.3,0.5,0.7 和0.9,速度過沖達到的最值分別為0.565,0.430,0.365,0.310 和0.235.不同β值達到最值所需的時間差異不大,穩(wěn)態(tài)值也近乎相同.由圖7(c)可見,β極大地影響著彈性剪切應(yīng)力的穩(wěn)態(tài)值.β越小,彈性剪切應(yīng)力的穩(wěn)態(tài)值越大,這是因為β越小,聚合物黏度越高,彈性力越大.可見,β是影響彈性力大小的一個主要因素.

4)α的影響

α是黏彈性溶液中與各向異性拖曳力相關(guān)的一個流變參數(shù),它對速度u和剪切應(yīng)力τxy的影響如圖6(d)和圖7(d)所示.由圖6(d)可見,α并不影響速度過沖達到的最值和達到最值所需的時間,但對速度達到的穩(wěn)態(tài)值產(chǎn)生一定的影響,即α越大,速度穩(wěn)態(tài)值越大.由圖7(d)可見,α對彈性剪切應(yīng)力過沖行為和穩(wěn)態(tài)值的影響均不大.

5)γ的影響

圖6(e)和圖7(e)分別描述了γ對速度u和彈性剪切應(yīng)力τxy的影響.可以看出,γ并不影響速度過沖達到的最值,但在一定程度上影響彈性剪切應(yīng)力過沖達到的最值.γ越大,彈性剪切應(yīng)力過沖達到的最值越大.這是因為γ越大,流體拉伸松弛時間越長.γ對速度達到的穩(wěn)態(tài)值有一定的影響,即γ越大,速度穩(wěn)態(tài)值越小.γ也輕微影響彈性剪切應(yīng)力達到的穩(wěn)態(tài)值.隨著γ的增大,彈性剪切應(yīng)力穩(wěn)態(tài)值輕微地增大.

6)Q的影響

Q表示聚合物分子鏈的臂數(shù).Q對速度u和剪切應(yīng)力τxy的影響如圖6(f)和圖7(f)所示.由圖6(f)可見,Q并不影響速度過沖達到的最值和達到最值所需的時間,但它輕微影響速度達到的穩(wěn)態(tài)值,即Q越大,速度穩(wěn)態(tài)值輕微地變小.另外,由圖7(f)可知,Q也輕微地影響彈性剪切應(yīng)力的過沖行為,但其對應(yīng)力穩(wěn)態(tài)值的影響并不大.

4.2 黏彈性液滴撞擊固壁

接下來,研究基于XPP 模型的黏彈性液滴撞擊固壁問題,該問題含有復雜的自由面.值得注意的是,對于含自由面的黏彈性流體的數(shù)值模擬,若直接使用SPH 方法進行模擬,則會產(chǎn)生張力不穩(wěn)定性(見圖8),并最終導致模擬過程的中斷.為了解決該問題,我們特別采用了第3.3 節(jié)介紹的人工應(yīng)力.

圖8 XPP 黏彈性液滴產(chǎn)生張力不穩(wěn)定性的SPH 結(jié)果Fig.8.SPH results of a XPP viscoelastic droplet with tensile instability.

圖9 給出了液滴撞擊固壁問題的計算模型.初始時刻的液滴為圓形,直徑d00.02 m,圓心位于(0 m,0.04 m).固壁幾何尺寸取為{(x,y):-0.3 ≤x≤0.3,y0}.當t=0 時,液滴具有初始下落速度V1 m·s-1;當t＞ 0 時,液滴在重力加速度g的作用下開始加速下落,然后撞擊固壁,發(fā)生鋪展變形.對于該問題,其鋪展變形主要受初始下降速度和重力作用的驅(qū)動,表面張力對其影響很小[21],因此本文并不考慮表面張力.

圖9 液滴撞擊固壁問題的計算模型Fig.9.Computational model of droplet impacting solid wall problem.

模擬中,液滴密度ρ1000 kg·m-3,取向松弛時間λ10.02 s.液滴總黏度η4 Pa·s-1,其中牛頓溶劑黏度ηs0.4 Pa·s-1,聚合物黏度ηp3.6 Pa·s-1.分別選取液滴直徑d0和下降速度V作為特征長度和特征速度,則無量綱Reynolds數(shù)Re=5,Weissenberg數(shù)Wi=1,溶劑黏度比β=0.1.其他XPP 流變參數(shù)為:α=0.01,γ=0.9,Q=4.為了解決張力不穩(wěn)定性,人工應(yīng)力參數(shù)ε0.5.粒子初始間距設(shè)置為Δx=0.0002,對應(yīng)于7845 個流體粒子、301 個邊界粒子和903 個虛擬粒子.核函數(shù)采用5 次樣條函數(shù),光滑長度h1.5Δx.聲速c12.5 m·s-1,時間步長 Δt1.0×10-6s.對于聲速c,分別取c20 m·s-1和 40 m·s-1進行數(shù)值實驗,所得模擬結(jié)果和 12.5 m·s-1的模擬結(jié)果沒有顯著差異,表明聲速c取值為 12.5 m·s-1,已能夠保證弱可壓流體的流動行為充分接近不可壓流體的行為.對于該算例,本文在聯(lián)想深騰6800 型高性能服務(wù)器上,使用1 個CPU(Intel Xeon Platinum 8273 @ 3.0 GHz)模擬20 萬個時間步,所耗時間約為5.3 h.

圖10 給出了利用本文改進SPH 方法模擬基于XPP 模型的黏彈性液滴在撞擊固壁之后6 個不同時刻的結(jié)果.值得注意的是,本文選用無量綱時間TtV/d0來記錄液滴的鋪展歷程.可以看出,液滴在撞擊固壁之后的流動過程分為三個階段:第一階段,液滴在撞擊固壁后迅速鋪展.該階段的液滴雖具有初始向下的速度,但在撞擊固壁之后,液滴左半部分流速向左,右半部分流速向右(見T=1.55 和1.85).第二階段,液滴鋪展到最大寬度后,受彈性力的作用逐漸開始收縮.在該階段,液滴的速度方向發(fā)生了轉(zhuǎn)變,即液滴左半部分流速向右,右半部分流速向左(見T=2.35,2.85 和3.45).第三階段,液滴收縮到最小寬度后,由于能量的耗散和重力的影響,再次緩慢鋪展(見T=5.00).

圖10 基于XPP 模型的液滴撞擊固壁問題的SPH 模擬(Re=5,Wi=1,β=0.1,α=0.01,γ=0.9,Q=4)(a)T=1.55;(b)T=1.85;(c)T=2.35;(d)T=2.85;(e)T=3.45;(f)T=5.00Fig.10.SPH simulation of droplet impacting solid wall problem based on XPP model(Re=5,Wi=1,β=0.1,α=0.01,γ=0.9,Q=4):(a)T=1.55;(b)T=1.85;(c)T=2.35;(d)T=2.85;(e)T=3.45;(f)T=5.00.

圖11 給出了XPP 液滴的鋪展寬度隨時間的變化.為了便于分析,這里也用液滴直徑d0對液滴鋪展寬度d(T)進行無量綱處理.從圖11 可清晰地觀察到液滴撞擊固壁之后的三個流動階段,其中液滴在T約為2.405 時達到最大鋪展寬度(～2.255);在T約為4.015 時收縮到最小寬度(～1.875).對于該算例,由于引入了人工應(yīng)力,因此為了驗證本文改進SPH 方法模擬含自由面的黏彈性流動問題的有效性,這里也展示了利用FDM 方法[34]模擬該問題的數(shù)值結(jié)果(見文獻[34]中的圖9(a)),其中自由界面的追蹤通過標記單元法來實現(xiàn).從圖11可以看出,利用本文改進SPH 方法得到的結(jié)果和FDM 解是基本符合的,從而表明改進SPH 方法在解決張力不穩(wěn)定性的同時,能夠得到可靠的模擬結(jié)果.值得注意的是,兩種數(shù)值結(jié)果在液滴收縮階段有微小的差異,可能歸因于:FDM 方法使用Possion方程求解壓力,而本文考慮弱可壓縮流體,使用狀態(tài)方程求解壓力.為了評價本文改進SPH 方法的數(shù)值收斂性,還特別增加了粒子初始間距分別為Δx=0.0001,0.00025 和0.0004 的數(shù)值試驗,而其他參數(shù)均保持不變.得到的XPP 液滴鋪展寬度隨時間的變化曲線如圖11 所示.可以看出,通過3 組不同粒子初始間距和原始粒子初始間距(Δx=0.0002)得到的液滴鋪展寬度結(jié)果基本一致,從而驗證了本文改進SPH 方法模擬黏彈性自由表面流是數(shù)值收斂的.

圖11 不同粒子間距下液滴鋪展寬度隨時間的變化以及SPH 解和FDM解[34]的比較Fig.11.Time changes of droplet spread width obtained by different initial particle spacings and comparison between SPH and FDM[34] solutions.

對于該問題,由于液滴不僅具有初始下落速度,而且受重力的作用加速下落,因此液滴會對固壁邊界產(chǎn)生較大的沖擊力.然而,由于液滴較大的黏度(η4 Pa·s-1),因此在撞擊壁面過程中并沒有出現(xiàn)明顯的小液滴飛濺.模擬過程中,也并未發(fā)現(xiàn)任何流體粒子流出固壁邊界(見圖10),表明本文提出的增強型邊界處理技術(shù)能夠有效地防止流體粒子穿透邊界.此外,SPH 結(jié)果與FDM解[34]是基本吻合的(見圖11),這進一步表明了本文提出的增強型邊界處理技術(shù)是合理的.總之,相較于傳統(tǒng)的排斥力法和虛擬粒子法,本文為SPH 方法的邊界處理提供了一種合理的、可供替代的邊界處理方案,且該方案具有程序設(shè)計簡單、避免構(gòu)造偽粒子進行Shepard 插值繁瑣操作的優(yōu)勢.

下面討論XPP 本構(gòu)模型中各流變參數(shù)對液滴鋪展寬度隨時間的變化的影響.

7)Re的影響

Re對液滴鋪展寬度的影響如圖12(a)所示.可以看到,隨著Re的增大,液滴的最大鋪展寬度顯著增加,達到最大鋪展寬度所需的時間相對滯后.對于Re=2,5 和10,液滴的最大鋪展寬度分別約為1.735,2.215 和2.705,達到最大鋪展寬度所需的時間分別約為1.850,2.005 和2.495.這是因為Re越大,液滴慣性力越大,導致液滴撞擊固壁后有更快的鋪展行為.另外,Re越大,液滴的收縮行為越弱.從圖12(a)可知,Re=2 的收縮行為明顯強于Re=10 的收縮行為.這是因為Re越大,液滴彈性力相對越小.由于Re越大,液滴的鋪展越快,導致液滴的最小收縮寬度和最后鋪展寬度也相應(yīng)地更大.

8)Wi的影響

圖12(b)展示了Wi對液滴鋪展寬度的影響.可以看出,增加Wi導致液滴達到的最大鋪展寬度顯著增加,達到最大鋪展寬度所需的時間也相對滯后.對于Wi=0.5,1,2,4 和8,液滴達到的最大鋪展寬度分別約 為2.005,2.215,2.455,2.700 和2.875,達到最大鋪展寬度的時間分別約為2.015,2.355,2.500,2.765 和2.955.另外,隨著Wi的增大,液滴的最小收縮寬度相應(yīng)增加.不同Wi得到的液滴最小收縮寬度分別約為1.805,1.845,1.975,2.110 和2.400.這是因為,隨著Wi的增大,黏彈性流體的剪切稀化增強,導致液滴的鋪展更快,因此液滴的鋪展寬度相應(yīng)增加.

圖12 不同流變參數(shù)下液滴鋪展寬度d(T)/d0 隨時間的變化(a)Re;(b)Wi;(c)β;(d)α;(e)γ;(f)QFig.12.Time change of droplet spread widthd(T)/d0 under different rheological parameters:(a)Re;(b)Wi;(c)β;(d)α;(e)γ;(f)Q.

9)β的影響

β對液滴鋪展寬度的影響如圖12(c)所示,β越大,液滴的收縮行為越弱,但最終鋪展寬度相差不大.特別是,對于β=0.9 時的液滴已觀察不到收縮現(xiàn)象.這是因為隨著β的增大,牛頓溶劑黏度增加,聚合物黏度降低,導致液滴的彈性力變小.在這種情形下,液滴表現(xiàn)出類似于牛頓流體的鋪展行為.另外,改變β只改變牛頓溶劑黏度占總黏度的比,而不改變總黏度的大小,因此液滴的最終鋪展寬度相差不大.

10)α的影響

各向異性流變參數(shù)α對液滴鋪展寬度的影響如圖12(d)所示.可以觀察到,改變α并不影響液滴達到的最大鋪展寬度,但在一定程度上影響液滴的收縮行為,即α越大,液滴的收縮行為越弱,液滴的最小收縮寬度因此相應(yīng)地增加.對于α=0,0.01,0.1,0.3 和0.5,液滴的最小收縮寬度分別約為1.835,1.837,1.875,1.925 和1.995.由于變?nèi)醯囊旱问湛s行為,液滴的最后鋪展寬度隨著α的增大而增大.

11)γ的影響

如圖12(e)所示,降低γ導致液滴的最大鋪展寬度顯著增加,達到最大鋪展寬度所需的時間也相對滯后.對于γ=0.9,0.7,0.5,0.3,0.2 和0.1,液滴的最大鋪展寬度分別約為2.215,2.235,2.265,2.335,2.400 和2.495,達到最大鋪展寬度的時間分別約為2.355,2.357,2.475,2.500,2.515 和2.535.另外,γ在一定程度上也影響著液滴的收縮行為,即γ越小,液滴的最小收縮寬度越大,這是因為降低γ導致液滴的拉伸松弛時間變短,聚合物分子鏈得不到有效的拉伸,因此液滴的彈性力變小,導致液滴的收縮行為變?nèi)?由于變?nèi)醯囊旱问湛s行為,液滴的最后鋪展寬度隨著γ的減小而增大.

12)Q的影響

圖12(f)展示了聚合物分子鏈臂數(shù)Q對液滴鋪展寬度的影響.可以看出,臂數(shù)越多,對流動行為的影響越弱.特別是,當Q＞ 4 時,已觀察不到其對液滴鋪展寬度的影響.另一方面,臂數(shù)越少,液滴的收縮行為越弱,達到的最終鋪展寬度越大.對于單臂的情形(Q=1),得到的液滴鋪展寬度明顯高于Q＞ 4 得到的液滴鋪展寬度.

5 結(jié)論

本文提出一種改進SPH 方法數(shù)值模擬了基于XPP 模型的黏彈性流動.為了提高計算精度,采用了一種不含核導數(shù)計算的核梯度修正離散格式.為了防止粒子穿透固壁,提出了一種增強型的邊界處理技術(shù).為了消除張力不穩(wěn)定性,將人工應(yīng)力耦合到動量守恒方程中.運用改進SPH 方法數(shù)值模擬了基于XPP 模型的黏彈性Poiseuille 流和黏彈性液滴撞擊固壁問題,深入分析了不同流變參數(shù)對流動過程的影響,結(jié)論如下.

1)對于黏彈性Poiseuille 流,通過比較SPH數(shù)值解和解析解以及評價數(shù)值解的收斂性,表明了本文改進SPH 方法模擬XPP 黏彈性流動問題是準確、有效的,且相對于傳統(tǒng)SPH 方法具有更高的計算精度.對于黏彈性液滴問題,利用增強型的邊界處理技術(shù)可有效地防止粒子穿透固壁;借助人工應(yīng)力可有效地消除張力不穩(wěn)定性.

2)XPP 流變參數(shù)對黏彈性Poiseuille 流動過程有著重要的影響,即Re,Wi和α越大,速度穩(wěn)態(tài)值越大;γ和Q越大,速度穩(wěn)態(tài)值越小;β越大,速度過沖現(xiàn)象越弱,但并不影響速度穩(wěn)態(tài)值的大小.Wi和β越大,彈性剪切應(yīng)力穩(wěn)態(tài)值越小;γ越大,彈性剪切應(yīng)力穩(wěn)態(tài)值輕微增大;Re,α和Q對彈性剪切應(yīng)力穩(wěn)態(tài)值的影響不大.

3)XPP 流變參數(shù)對黏彈性液滴撞擊固壁流動過程有著重要的影響,即Re和Wi越大,液滴鋪展越快,鋪展寬度越大.β越大,液滴收縮行為越弱,但并不影響液滴的最終鋪展寬度.α越大,液滴收縮行為越弱,鋪展寬度越大.γ越大,液滴收縮行為越強,鋪展寬度越小.Q越大,對液滴鋪展寬度的影響越弱.

4)本文改進SPH 方法能夠有效而準確地描述基于XPP 模型的黏彈性流體的復雜流變特性和自由面變化特征.