基于SPH方法的瞬態非等溫黏彈性Couette流的數值模擬

許曉陽, 趙雨婷

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

0 引 言

黏彈性流動廣泛存在于自然界和工業加工中,如擠壓成型、注塑成型等[1]．這些流動一般呈現非等溫狀態,因此有必要考慮溫度模型[2]．深入研究黏彈性流體的非等溫流動機理具有重要的學術價值和應用意義,但由于其多物理場耦合的復雜性和黏彈流變特性的復雜性,實驗和解析求解均難以完成,因此數值模擬成為一種可替代的研究途徑．

目前,基于網格的數值方法已被應用于此類問題的數值模擬中．Meburger等[3]用非結構網格有限體積法(finite volume method,FVM)模擬了非等溫黏彈性流體在高Weissenberg數下的流動．Moreno等[4]用變分多尺度穩定有限元法分析了溫度對黏彈性流動的影響．Gao[5]用相場模型對非等溫聚合物填充過程進行了有限元計算．

SPH[6-7]方法是一種Lagrange型的無網格數值方法．與基于網格的數值方法相比,SPH方法不僅免去了生成網格的繁瑣,而且避免了傳統Lagrange方法中網格扭曲、重構等問題[8]．近年來,隨著SPH方法計算精度的提高和穩定性的改善,該方法已被較多地應用于潰壩[9]、爆炸[10]、液體晃蕩[11]和流固耦合[12]等問題的數值模擬中．然而,對于黏彈性流動問題的數值模擬,由于其非線性黏彈本構方程的引入,極大地增加了問題求解的復雜度,因此相關文獻研究報道較少．2005年,Ellero等[13]首次進行了等溫條件下黏彈性流動的SPH模擬,其中采用Oldroyd-B模型來描述流體的黏彈性行為．Fang等[14]將SPH方法推廣到了含自由面的瞬態等溫黏彈性液滴下落問題的數值模擬中．Murashima等[15]提出了多尺度SPH方法研究支化聚合物熔體繞圓柱形障礙物的流動．楊波等[16]運用SPH方法數值模擬了基于Oldroyd-B模型的平面突然起動Couette流．Xu等[17]提出了一種改進SPH方法,對等溫黏彈性液滴和擠出脹大問題進行了研究．King等[18]將對數構象公式與應力分裂技術相結合,提出了模擬黏彈性流動的不可壓SPH方法．Duque-Daza等[19]將SPH與流體流動模擬的等效彈性勢能結合在一起,建立了黏彈性流體建模與計算的簡化方法．Vahabi等[20]用弱可壓SPH方法模擬了Oldroyd-B黏彈性溶液中兩個初始圓形的氣泡在上升過程中的相互作用．根據文獻調研結果,目前應用SPH方法對黏彈性流動問題開展的模擬研究大多是在等溫情況下進行的,而工業生產中的黏彈性流動通常是非等溫流動,因此有必要開展非等溫黏彈性流動的數值模擬．

本文在文獻[17]的基礎上,運用SPH方法對瞬態非等溫黏彈性流動問題進行了數值模擬研究．首先,模擬了等溫情況下基于Oldroyd-B模型的黏彈性Couette流動;隨后,將其擴展到非等溫情況下進行模擬,其中選用Reynolds指數模型來評估黏度和松弛時間的溫度依賴．通過與有限體積方法解的比較和對數值收斂性的評價,驗證了SPH方法模擬非等溫黏彈性流動問題的準確性和有效性．討論了非等溫流動相較于等溫流動的不同流動特征,分析了溫度依賴系數、Péclet數等對黏彈性流動過程的影響．

1 控 制 方 程

在Lagrange坐標系下,非等溫黏彈性流體的控制方程組可表示為[21]

(1)

(2)

(3)

式中,d/dt為物質導數,ρ為流體密度,t為時間,u為速度,σ為總應力張量,F為外力,cp為定壓比熱容,T為溫度,κ為導熱系數．值得注意的是,本文僅考慮純受迫流動,浮升力相對于慣性力很小,因此忽略了自然對流[22]．

總應力張量σ可分解為各向同性壓力、溶劑貢獻和聚合物貢獻之和:

σ=-pI+2ηsd+τ,

(4)

其中,I為單位矩陣,ηs為溶劑黏度,τ表示彈性偏應力張量,d為形變率張量,

(5)

通過對總應力張量的分解,方程(2)可表示為

(6)

1.1 Oldroyd-B本構方程

對于彈性偏應力張量τ,考慮如下微分形式的Oldroyd-B本構方程:

(7)

式中,λ表示流體的松弛時間;ηp表示聚合物黏度．ηs和ηp滿足如下關系:

η=ηs+ηp,

(8)

(9)

綜合式(7)和(9),有

(10)

對于流體松弛時間和黏度的溫度依賴,選用Reynolds指數模型進行建模和計算:

ηs(T)=ηse-φ(T-T0),

(11)

ηp(T)=ηpe-φ(T-T0),

(12)

λ(T)=λe-φ(T-T0),

(13)

其中,φ為溫度依賴系數,T0為參考溫度．另外,為了更好地描述流動,引入以下無量綱數:Péclet數Pe=ρcpVL/κ、Reynolds數Re=ρVL/η、Weissenberg數Wi=λV/L、溶劑黏度比β=ηs/η,其中V和L分別表示流動的特征速度和特征長度．

1.2 狀態方程

SPH方法常采用狀態方程來模擬弱可壓縮流動．本文使用的狀態方程如下[23]:

(14)

其中,ρ0為參考密度,c為聲速,γ=7．為了保證弱可壓縮流動行為足夠接近不可壓縮流動行為,聲速c通常選擇比最大流體速度大10倍左右．

2 SPH方法

2.1 SPH離散

對于質量和動量守恒方程,本文采用的SPH離散為[24]

(15)

(16)

式中,i表示粒子編號,j為粒子i的相鄰粒子,Wij=W(|rij|,h)為核函數,rij=ri-rj,ri為粒子i的位置,h為光滑長度．由于五次樣條核函數具有連續的二階導數,其精度和穩定性均較好,故本文選用它作為核函數,其表達式如下:

(17)

其中,q=r/h,對于二維問題W0取值7/(478πh2)．

對于Oldroyd-B本構方程,其SPH離散取決于速度偏導數的離散．引入

(18)

于是,Oldroyd-B本構方程的SPH離散可表示為

(19)

傳統SPH方法的計算精度低,限制了其向更復雜流動問題的應用．為了提高計算精度,本文對SPH離散(15)—(19)中的核函數梯度進行修正,即利用修正矩陣[25]

(20)

將原核函數梯度修正為

(21)

利用新的核梯度(21)進行計算,可提高SPH方法的計算精度．關于核梯度修正SPH離散的更多內容,可參閱文獻[24]．

2.1.1 溫度方程的離散

若要實現非等溫黏彈性流動問題的SPH模擬,還需對溫度方程進行SPH離散．對于溫度方程(3),其等號右端部分涉及Laplace算子的離散化．通常,可采用SPH散度算子和梯度算子的離散對Laplace算子進行離散,但其所產生核函數的二階導數對粒子的無序性十分敏感,容易引起數值計算的不穩定．借鑒Shao和Lo[26]對Laplace算子離散的思想,本文對溫度方程采用如下的SPH離散格式:

(22)

2.2 邊界處理

2.2.1 固壁邊界

由于邊界或鄰近邊界處的粒子在積分時會被邊界截斷,這種單向影響會導致求解錯誤,因此有必要對固壁邊界進行特殊處理．

本文采用由固壁粒子和固壁外虛粒子組成的增強型處理技術[9]．首先,沿固壁邊界以粒子初始間距Δx布置一層固壁粒子．在計算過程中,固壁粒子參與控制方程中速度、壓力和彈性應力的計算,以有效地防止粒子穿透邊界．固壁粒子的密度和位置均保持不變,壓力和彈性應力通過其支持域內流體粒子壓力和彈性應力的正則化插值得到:

(23)

其中,B表示壓力或彈性應力,i表示固壁粒子,j表示與固壁粒子i相鄰的流體粒子．

其次,為彌補邊界處粒子在積分時被邊界截斷的問題,在固壁邊界外以粒子初始間距Δx額外布置3層固壁外虛粒子．與邊界粒子類似,固壁外虛粒子的密度和位置也保持不變．為滿足壓力和彈性應力的Neumann邊界條件,固壁外虛粒子的壓力和彈性應力設置為與邊界法向上相鏈接的固壁粒子值相同,如圖1所示．

圖1 固壁邊界處理示意圖Fig.1 Sketch of the wall boundary treatment

2.2.2 周期邊界

對于黏彈性Couette流動的模擬,采用周期邊界條件進行無限長平板的處理．在粒子運動過程中,沿x方向施加周期邊界條件,其實施步驟如下:在每一時間步,通過右側邊界離開流域的SPH粒子將從左側邊界在相同的y位置重新插入流域;左側邊界周圍的流體粒子被定義為右側邊界附近粒子的相鄰粒子,反之亦然．

2.3 時間積分

由于蛙跳積分法具有存儲需求量低、計算效率高的優點,因此本文選用該方法進行時間積分．關于該積分格式的詳細內容,可參閱文獻[8]．

3 黏彈性Couette流模擬

鑒于目前尚未見到運用SPH方法對基于Oldroyd-B模型的非等溫黏彈性Couette流動問題進行數值模擬的文獻報道,因此本文選用該問題作為數值算例,并討論非等溫流動相較于等溫流動的不同流動特征,分析各物理參數對流動過程的影響．

3.1 等溫流動

首先考慮等溫流動．圖2給出了該問題的幾何模型和初始狀態．當t=0時,流體是靜止的,且位于兩塊固定、垂直間距為L的無限長平板之間;當t>0時,上平板突然以平行于x軸正方向的恒定速度u0運動．

圖2 黏彈性Couette流的幾何區域Fig.2 The geometric region of the viscoelastic Couette flow

模擬中,沿x軸方向采用周期邊界約束,故計算區域可取為Lx×Ly=0.5×1．上平板速度u0=1,流體密度ρ=1,各物理參數取值為Re= 0.1,Wi=0.1,β= 0.1,其中特征長度L=1,特征速度V=1．粒子初始間距設置為Δx=0.02,對應于1 274個流體粒子,64個固壁粒子和192個固壁外虛粒子．核函數采用五次樣條函數,光滑長度采用與Ellero等[13]模擬黏彈性Poiseuille流相同的取值,即h=0.9Δx．聲速c=10V,時間步長Δt=1.0×10-5．選用3個空間點A,B和C記錄流體的瞬態流動行為,它們到下平板的高度分別為0.25L,0.5L和0.75L．對于該問題,本文用聯想計算機(Intel(R) Core(TM) i7-9700 CPU@3.60 GHz)模擬10萬個時間步,所耗時間約為983 s．

圖3給出了黏彈性Couette流動在t=1時刻的粒子分布和速度分布;而圖4給出了A,B和C三點處的速度u和彈性剪切應力τxy隨時間的變化情況,并將本文SPH數值解和解析解[27]進行了比較．可以看出,由于受彈性作用的影響,Couette流發生了明顯的速度過沖．A,B和C三點處的速度u均經歷了先增大、后減小,逐漸達到穩態的過程;離上平板越近,穩態時的速度值越大．點C處的穩態速度值明顯高于點A和B處的穩態速度值．對于剪切應力,離上平板越近,應力發生的時間越早,但穩態時的剪切應力值幾乎相同．無論對于速度u還是剪切應力τxy,SPH數值解與解析解均吻合,從而驗證了本文SPH方法求解黏彈性流體的有效性和準確性．

圖3 等溫黏彈性Couette流動在t=1時刻的粒子分布和速度分布Fig.3 The particle distribution and the velocity distribution of the isothermal viscoelastic Couette flow at t=1

(a) 速度u (b) 彈性剪切應力τxy(a) Velocity u (b) Elastic shear stress τxy圖4 等溫黏彈性Couette流動A,B,C三點處的速度u和彈性剪切應力τxy隨時間的變化:SPH數值解與解析解的對比Fig.4 Isothermal viscoelastic Couette flows at 3 points A,B and C:a comparison of SPH and analytical solutions for velocity u and elastic shear stress τxy

圖5給出了Re=1和10時點B處的速度u和彈性剪切應力τxy隨時間的變化情況,其他參數保持不變．可以看出,即便對于較大的Re,本文SPH方法也能夠精確穩定地模擬．對于Re=1和10,所得SPH數值解和解析解均吻合較好．另外,Re越大,速度過沖現象越不明顯,但穩態值幾乎相同．這是因為:Re越大,流動慣性力越占優,導致彈性力所占比重下降,因此速度的過沖現象減弱．Re越大,剪切應力穩態值越小,這是因為:Re越大,流體黏度越小,導致應力計算的穩態值也越小．

(a) 速度u (b) 彈性剪切應力τxy(a) Velocity u (b) Elastic shear stress τxy圖5 Re=1和10時,等溫黏彈性Couette流動的SPH模擬:點B處的速度u和彈性剪切應力τxy隨時間的變化Fig.5 SPH simulations of isothermal viscoelastic Couette flows with Re=1 and 10:the time changes of velocity u and elastic shear stress τxy at point B

3.2 非等溫流動

接下來,我們將黏彈性Couette流動擴展到非等溫情況下進行模擬．取流體溫度和下平板溫度Tf=313,上平板溫度Tw=443,溫度依賴系數φ=0.001,Péclet數Pe=100．模擬過程中,其他參數均設置為與等溫流動相同．

圖6給出了非等溫情況下點B處的速度u和彈性剪切應力τxy隨時間的變化情況,并將非等溫結果與等溫結果進行了比較．在流動早期階段,非等溫和等溫情況的流動特征幾乎相同;但在t>1之后,溫度的引入對流體的流動特征產生了一定的影響．等溫流動在t=1時刻已達到了穩態,但非等溫流動的速度和剪切應力值均出現輕微下降,這表明非等溫流動并未達到熱平衡．該下降趨勢一直持續下去,即便在t=2時刻,非等溫流動仍未達到穩態．為了分析溫度依賴系數φ對流動過程的影響,本文特別增加了φ=0.005和0.010的數值試驗,其他參數均保持不變．從圖 6可以看出,φ越大,流動后期速度和剪切應力的下降趨勢越明顯．這表明:溫度依賴系數對流動過程的影響隨著φ的增加而愈加明顯．但即便對于較大的φ= 0.010,流動在t=2仍未達到穩態．

(a) 速度u (b) 彈性剪切應力τxy(a) Velocity u (b) Elastic shear stress τxy圖6 非等溫黏彈性Couette流的SPH模擬:φ對點B處速度u和彈性剪切應力τxy隨時間變化的影響Fig.6 The SPH simulation of the non-isothermal viscoelastic Couette flow:the effect of temperature dependent parameter φ on the time changes of velocity u and elastic shear stress τxy at point B

由于φ越大,流動后期速度和彈性剪切應力的下降趨勢越明顯,因此接下來本文選取較大的φ= 0.01,研究Pe對流動過程的影響．圖7給出了3個不同Pe(Pe=1,20,100)下點B處的速度u和彈性剪切應力τxy隨時間的變化情況．可以看出,Pe越小,流動達到穩態所需的時間越短．對于Pe=1,速度u和彈性剪切應力τxy在t=1時刻已達到穩態,但Pe=20和100還沒有．這是因為,Péclet數表示對流速率與擴散速率之比,Pe越小,比熱容越小,流體熱擴散越快,流動會更快地達到熱平衡．為了檢驗Pe=1時的流動在t=1時刻是否已達到熱平衡,圖8給出了Couette流動在6個不同時刻的溫度關于y的分布情況．流動初期,熱量由上平板向下傳遞,流體溫度隨之升高;在t=1時刻,流體溫度呈線性分布,這表明Couette流動在此時刻已經達到熱平衡．值得注意的是,對于Pe=20和100的算例,流動在達到熱平衡后也可以達到穩態,只是所需的時間更長．

(a) 速度u (b) 彈性剪切應力τxy(a) Velocity u (b) Elastic shear stress τxy圖7 非等溫黏彈性Couette流的SPH模擬:Pe對點B處速度u和彈性剪切應力τxy隨時間變化的影響Fig.7 The SPH simulation of the non-isothermal viscoelastic Couette flow:the effect of Pe on the time changes of velocity u and elastic shear stress τxy at point B

圖8 非等溫黏彈性Couette流的SPH模擬(Pe=1,φ=0.01):6個不同時刻的溫度關于y的分布 圖9 利用SPH得到的溫度分布圖的收斂性分析及其與FVM解的比較Fig.8 The SPH simulation of the non-isothermal viscoelastic Couette flow (Pe=1,φ=0.01):the temperature distribution vs.y at 6 different moments Fig.9 Convergence analysis of the temperature distribution obtained with the SPH and the comparison with the FVM solution

為了驗證本文SPH方法求解非等溫黏彈性流動問題的有效性,圖9給出了Pe=1和φ=0.01時,利用SPH方法和FVM方法[28]得到的3個不同時刻溫度分布圖的比較．不難看出,SPH解和FVM解吻合較好．另外,為了驗證本文SPH方法求解非等溫黏彈性流動問題的數值收斂性,本文特別增加了粒子初始間距分別為Δx=0.01,0.012 5,0.04的數值試驗,其他參數均保持不變,得到的溫度分布圖如圖 9所示．可以看出,通過3組不同粒子初始間距和原始粒子初始間距(Δx=0.02)得到的溫度分布結果無顯著差異．上述結果表明,本文SPH方法對于非等溫黏彈性流動問題的模擬是準確、收斂的．

3.2.1β的影響

β表示溶劑黏度與流體總黏度之比．為了進一步研究β對流動過程的影響,圖10給出了不同β(β=0.05,0.10,0.30,0.50)取值下點B處的速度u和彈性剪切應力τxy隨時間的變化情況．圖10中的各物理參數取值如下:Re=0.1,Wi=0.1,Pe=1,φ=0.01,因為此時流動可較快地達到熱平衡．可以看出,β并不影響速度的穩態值,但顯著影響速度的瞬態行為．β越大,速度的過沖現象越不明顯．這一結果表明:β較小,即溶劑黏度比總黏度小得多時,此時彈性占優,因此對黏彈性瞬態流動的影響較大;反之,當β增大時,彈性作用對流動的影響減小．對于彈性剪切應力τxy,其穩態值隨著β的增大而減小,這是因為:β越大,流體越接近于Newton流體,彈性力越弱．

(a) 速度u (b) 彈性剪切應力τxy(a) Velocity u (b) Elastic shear stress τxy圖10 非等溫黏彈性Couette流的SPH模擬:β對點B處速度u和彈性剪切應力τxy隨時間變化的影響Fig.10 The SPH simulation of the non-isothermal viscoelastic Couette flow:the effect of β on the time changes of velocity u and elastic shear stress τxy at point B

3.2.2Wi的影響

Wi是黏彈性流動的另一個重要參數．它比較了彈性力與黏性力,通常由流體的應力松弛時間與具體的加工時間的關系給出．為了研究Wi對流動產生的影響,圖 11給出了Re=0.1,β= 0.1,Pe=1和φ=0.01時,不同Wi取值下(Wi=0.05,0.10,0.20,0.50)點B處的速度u和彈性剪切應力τxy隨時間的變化情況．可以看出,Wi對速度瞬態行為的影響較大,但并不影響穩態值的大小．Wi越大,速度過沖現象越明顯,達到最值所需時間也越長．對于Wi=0.05,0.10,0.20,0.50,速度過沖達到的最值分別為umax≈0.520,0.575,0.615,0.640,而達到最值所需的時間分別為t≈0.06,0.10,0.12,0.20．對于較大的Wi=0.5,還觀察到2次明顯的速度過沖．

(a) 速度u (b) 彈性剪切應力τxy(a) Velocity u (b) Elastic shear stress τxy圖11 非等溫黏彈性Couette流的SPH模擬:Wi對點B處速度u和彈性剪切應力τxy隨時間變化的影響Fig.11 The SPH simulation of the non-isothermal viscoelastic Couette flow:the effect of Wi on the time changes of velocity u and elastic shear stress τxy at point B

對于彈性剪切應力,我們觀察到隨著Wi的增大,流動前期的剪切應力值逐漸變小．特別是對于Wi=0.05的試驗,其剪切應力很快達到穩態;但當Wi=0.5時,其剪切應力比Wi=0.05的小,且在持續增長．Wi越大,剪切應力達到穩態值所需的時間越長,但它并不改變其穩態值的大小．對于Wi=0.05,0.10,0.20,0.50,剪切應力達到穩態值所需的時間分別為0.69,0.78,0.91,1.51．這是因為,Wi越大,流體的應力松弛時間越長．在這種情況下,流動需要更長的時間來達到穩態．

4 結 論

本文運用SPH方法對瞬態非等溫黏彈性流動問題進行了數值模擬研究．首先,模擬了等溫情況下基于Oldroyd-B模型的黏彈性Couette流動;隨后,將其擴展到非等溫情況下進行模擬．討論了非等溫流動相較于等溫流動的不同流動特征,分析了溫度依賴系數、Péclet數等對流動過程的影響．本文所得結論如下:

1) 通過與有限體積方法解的比較和對數值收斂性的評價,表明了本文SPH方法可準確有效地模擬非等溫黏彈性流動問題．

2) 溫度的引入對黏彈性流動過程有影響．溫度依賴系數φ越大,流動后期速度和剪切應力下降得越快．Pe越小,流動達到熱平衡所需的時間越短．

3)β不影響速度的穩態值,但顯著地影響其瞬態行為．β越大,速度過沖現象越不明顯．剪切應力的穩態值隨著β的增大而減小．

4)Wi不影響速度和剪切應力的穩態值,但影響其瞬態行為．Wi越大,速度過沖現象越明顯,速度和剪切應力達到穩態值所需的時間越長．