許曉陽, 趙雨婷

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

0 引 言

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

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

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

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

1 控 制 方 程

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

(1)

(2)

(3)

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

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

σ=-pI+2ηsd+τ,

(4)

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

(5)

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

(6)

1.1 Oldroyd-B本構(gòu)方程

對于彈性偏應(yīng)力張量τ,考慮如下微分形式的Oldroyd-B本構(gòu)方程:

(7)

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

η=ηs+ηp,

(8)

(9)

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

(10)

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

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

(11)

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

(12)

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

(13)

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

1.2 狀態(tài)方程

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

(14)

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

2 SPH方法

2.1 SPH離散

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

(15)

(16)

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

(17)

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

對于Oldroyd-B本構(gòu)方程,其SPH離散取決于速度偏導數(shù)的離散．引入

(18)

于是,Oldroyd-B本構(gòu)方程的SPH離散可表示為

(19)

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

(20)

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

(21)

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

2.1.1 溫度方程的離散

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

(22)

2.2 邊界處理

2.2.1 固壁邊界

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

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

(23)

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

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

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

2.2.2 周期邊界

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

2.3 時間積分

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

3 黏彈性Couette流模擬

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

3.1 等溫流動

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

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

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

圖3給出了黏彈性Couette流動在t=1時刻的粒子分布和速度分布;而圖4給出了A,B和C三點處的速度u和彈性剪切應(yīng)力τxy隨時間的變化情況,并將本文SPH數(shù)值解和解析解[27]進行了比較．可以看出,由于受彈性作用的影響,Couette流發(fā)生了明顯的速度過沖．A,B和C三點處的速度u均經(jīng)歷了先增大、后減小,逐漸達到穩(wěn)態(tài)的過程;離上平板越近,穩(wěn)態(tài)時的速度值越大．點C處的穩(wěn)態(tài)速度值明顯高于點A和B處的穩(wěn)態(tài)速度值．對于剪切應(yīng)力,離上平板越近,應(yīng)力發(fā)生的時間越早,但穩(wěn)態(tài)時的剪切應(yīng)力值幾乎相同．無論對于速度u還是剪切應(yīng)力τxy,SPH數(shù)值解與解析解均吻合,從而驗證了本文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) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖4 等溫黏彈性Couette流動A,B,C三點處的速度u和彈性剪切應(yīng)力τxy隨時間的變化:SPH數(shù)值解與解析解的對比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和彈性剪切應(yīng)力τxy隨時間的變化情況,其他參數(shù)保持不變．可以看出,即便對于較大的Re,本文SPH方法也能夠精確穩(wěn)定地模擬．對于Re=1和10,所得SPH數(shù)值解和解析解均吻合較好．另外,Re越大,速度過沖現(xiàn)象越不明顯,但穩(wěn)態(tài)值幾乎相同．這是因為:Re越大,流動慣性力越占優(yōu),導致彈性力所占比重下降,因此速度的過沖現(xiàn)象減弱．Re越大,剪切應(yīng)力穩(wěn)態(tài)值越小,這是因為:Re越大,流體黏度越小,導致應(yīng)力計算的穩(wěn)態(tài)值也越?。?/p>

(a) 速度u (b) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖5 Re=1和10時,等溫黏彈性Couette流動的SPH模擬:點B處的速度u和彈性剪切應(yīng)力τ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,溫度依賴系數(shù)φ=0.001,Péclet數(shù)Pe=100．模擬過程中,其他參數(shù)均設(shè)置為與等溫流動相同．

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

(a) 速度u (b) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖6 非等溫黏彈性Couette流的SPH模擬:φ對點B處速度u和彈性剪切應(yīng)力τ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

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

(a) 速度u (b) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖7 非等溫黏彈性Couette流的SPH模擬:Pe對點B處速度u和彈性剪切應(yīng)力τ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個不同時刻的溫度關(guān)于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方法求解非等溫黏彈性流動問題的數(shù)值收斂性,本文特別增加了粒子初始間距分別為Δx=0.01,0.012 5,0.04的數(shù)值試驗,其他參數(shù)均保持不變,得到的溫度分布圖如圖 9所示．可以看出,通過3組不同粒子初始間距和原始粒子初始間距(Δx=0.02)得到的溫度分布結(jié)果無顯著差異．上述結(jié)果表明,本文SPH方法對于非等溫黏彈性流動問題的模擬是準確、收斂的．

3.2.1β的影響

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

(a) 速度u (b) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖10 非等溫黏彈性Couette流的SPH模擬:β對點B處速度u和彈性剪切應(yīng)力τ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是黏彈性流動的另一個重要參數(shù)．它比較了彈性力與黏性力,通常由流體的應(yīng)力松弛時間與具體的加工時間的關(guān)系給出．為了研究Wi對流動產(chǎn)生的影響,圖 11給出了Re=0.1,β= 0.1,Pe=1和φ=0.01時,不同Wi取值下(Wi=0.05,0.10,0.20,0.50)點B處的速度u和彈性剪切應(yīng)力τxy隨時間的變化情況．可以看出,Wi對速度瞬態(tài)行為的影響較大,但并不影響穩(wěn)態(tài)值的大?。甒i越大,速度過沖現(xiàn)象越明顯,達到最值所需時間也越長．對于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) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖11 非等溫黏彈性Couette流的SPH模擬:Wi對點B處速度u和彈性剪切應(yīng)力τ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

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

4 結(jié) 論

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

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

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

3)β不影響速度的穩(wěn)態(tài)值,但顯著地影響其瞬態(tài)行為．β越大,速度過沖現(xiàn)象越不明顯．剪切應(yīng)力的穩(wěn)態(tài)值隨著β的增大而減小．

4)Wi不影響速度和剪切應(yīng)力的穩(wěn)態(tài)值,但影響其瞬態(tài)行為．Wi越大,速度過沖現(xiàn)象越明顯,速度和剪切應(yīng)力達到穩(wěn)態(tài)值所需的時間越長．