許曉陽, 趙雨婷

(西安科技大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,西安 710054)

0 引 言

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

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

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

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

1 控 制 方 程

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

(1)

(2)

(3)

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

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

σ=-pI+2ηsd+τ,

(4)

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

(5)

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

(6)

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

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

(7)

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

η=ηs+ηp,

(8)

(9)

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

(10)

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

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

(11)

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

(12)

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

(13)

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

1.2 狀態(tài)方程

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

(14)

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

2 SPH方法

2.1 SPH離散

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

(15)

(16)

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

(17)

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

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

(18)

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

(19)

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

(20)

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

(21)

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

2.1.1 溫度方程的離散

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

(22)

2.2 邊界處理

2.2.1 固壁邊界

由于邊界或鄰近邊界處的粒子在積分時(shí)會(huì)被邊界截?cái)?這種單向影響會(huì)導(dǎo)致求解錯(cuò)誤,因此有必要對(duì)固壁邊界進(jìn)行特殊處理．

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

(23)

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

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

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

2.2.2 周期邊界

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

2.3 時(shí)間積分

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

3 黏彈性Couette流模擬

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

3.1 等溫流動(dòng)

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

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

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

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

圖3 等溫黏彈性Couette流動(dòng)在t=1時(shí)刻的粒子分布和速度分布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流動(dòng)A,B,C三點(diǎn)處的速度u和彈性剪切應(yīng)力τxy隨時(shí)間的變化:SPH數(shù)值解與解析解的對(duì)比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時(shí)點(diǎn)B處的速度u和彈性剪切應(yīng)力τxy隨時(shí)間的變化情況,其他參數(shù)保持不變．可以看出,即便對(duì)于較大的Re,本文SPH方法也能夠精確穩(wěn)定地模擬．對(duì)于Re=1和10,所得SPH數(shù)值解和解析解均吻合較好．另外,Re越大,速度過沖現(xiàn)象越不明顯,但穩(wěn)態(tài)值幾乎相同．這是因?yàn)?Re越大,流動(dòng)慣性力越占優(yōu),導(dǎo)致彈性力所占比重下降,因此速度的過沖現(xiàn)象減弱．Re越大,剪切應(yīng)力穩(wěn)態(tài)值越小,這是因?yàn)?Re越大,流體黏度越小,導(dǎo)致應(yīng)力計(jì)算的穩(wěn)態(tài)值也越?。?/p>

(a) 速度u (b) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖5 Re=1和10時(shí),等溫黏彈性Couette流動(dòng)的SPH模擬:點(diǎn)B處的速度u和彈性剪切應(yīng)力τxy隨時(shí)間的變化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 非等溫流動(dòng)

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

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

(a) 速度u (b) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖6 非等溫黏彈性Couette流的SPH模擬:φ對(duì)點(diǎn)B處速度u和彈性剪切應(yīng)力τxy隨時(shí)間變化的影響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

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

(a) 速度u (b) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖7 非等溫黏彈性Couette流的SPH模擬:Pe對(duì)點(diǎn)B處速度u和彈性剪切應(yīng)力τxy隨時(shí)間變化的影響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個(gè)不同時(shí)刻的溫度關(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

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

3.2.1β的影響

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

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

(a) 速度u (b) 彈性剪切應(yīng)力τxy(a) Velocity u (b) Elastic shear stress τxy圖11 非等溫黏彈性Couette流的SPH模擬:Wi對(duì)點(diǎn)B處速度u和彈性剪切應(yīng)力τxy隨時(shí)間變化的影響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

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

4 結(jié) 論

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

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

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

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

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