非等溫黏彈性聚合物流體圓柱繞流的高精度數(shù)值模擬*

莊 昕， 劉付軍， 孫艷萍， 王惠玲

（1. 河南工程學(xué)院 理學(xué)院，鄭州 451191；2. 山西財(cái)經(jīng)大學(xué) 應(yīng)用數(shù)學(xué)學(xué)院，太原 030006）

引 言

聚合物制品由于其優(yōu)良的力學(xué)性能，在航空航天、國防軍工、汽車制造以及生物醫(yī)藥等領(lǐng)域的應(yīng)用日益廣泛.聚合物流體(聚合物溶液)是一類典型的非Newton 黏彈性流體，在加工成型過程中會出現(xiàn)既非Hooke 彈性體，又非Newton 黏性流體的復(fù)雜奇異流變行為，這種特殊流變行為對其最終制品的綜合性能產(chǎn)生重要影響.因此，對各向異性黏彈性聚合物流體在不同流動區(qū)域流變行為的研究，一直是流變學(xué)研究領(lǐng)域的難點(diǎn)和熱點(diǎn)問題，具有非常重要的理論意義和實(shí)際工程應(yīng)用價(jià)值[1].然而，傳統(tǒng)的解析求解和實(shí)驗(yàn)分析都很難解決實(shí)際工程問題.隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值模擬成為研究聚合物流動過程的一種強(qiáng)有力的手段.與Newton 流體和非Newton 黏性流體相比，非等溫黏彈性流動問題的控制方程包括質(zhì)量方程、動量方程、能量方程和本構(gòu)方程.復(fù)雜非線性黏彈本構(gòu)方程的引入不僅增加了控制方程未知變量的個數(shù)，同時也使得數(shù)值模擬的計(jì)算難度顯著增大.

建立黏彈流變本構(gòu)方程的途徑主要有唯象性方法和分子論方法兩種[2].經(jīng)典的唯象性黏彈本構(gòu)方程有線性Maxwell 模型、非線性O(shè)ldroyd-B 模型等[3].和唯象性方法相比，分子論方法建立的本構(gòu)模型能從機(jī)理上闡明黏彈性材料的奇異流變現(xiàn)象，能更清晰地揭示黏彈性流體在流動過程中的復(fù)雜力學(xué)響應(yīng).對于聚合物稀溶液，基于分子論方法建立的有限伸展FENE 類啞鈴模型[4-5]能夠有效描述其非線性黏彈流變行為.對于聚合物濃厚體系，大分子鏈間相互纏結(jié)的機(jī)理和本質(zhì)非常復(fù)雜.基于Lodge 網(wǎng)絡(luò)理論的PTT (Phan-Thien Tanner)模型[6]雖然能夠預(yù)測濃溶液的剪切變稀和應(yīng)力過沖現(xiàn)象，但此模型在剪切流動初期時會出現(xiàn)虛假振蕩.Mcleish 和Larson[7]首次基于聚合物熔體蠕動和管道理論提出了非線性PP(Pom-Pom)模型，此模型可以描述支化聚合物熔體的奇異流變性質(zhì)，是聚合物濃厚體系流變本構(gòu)方程研究的一個重大突破.然而，PP 模型除了在穩(wěn)態(tài)拉伸流動中出現(xiàn)間斷現(xiàn)象外，其描述取向的方程在高剪切流動中還存在無界現(xiàn)象.Verbeeten 等[8]基于PP 模型提出了著名的XPP 黏彈本構(gòu)模型.研究結(jié)果表明[9-12]，XPP 模型不僅能夠正確地描述聚合物濃溶液的拉伸和剪切行為，還能有效克服應(yīng)力奇點(diǎn)問題，使可計(jì)算的We數(shù)更大.因此，近年來基于分子論方法建立的XPP 模型已被應(yīng)用于聚合物熔體注塑成型[13-14]和擠出成型[15]的數(shù)值預(yù)測中.然而，上述基于XPP 非線性黏彈本構(gòu)方程的模擬研究均在等溫條件下進(jìn)行.事實(shí)上，聚合物成型過程是非等溫非穩(wěn)態(tài)流動過程，溫度對流動狀態(tài)和受力狀態(tài)的影響不可忽略，一些非等溫條件下聚合物溶液流動問題的模擬也只是采用純黏性本構(gòu)模型[16-17].

目前，有限元法[9-11,18]、有限體積法[12-13]和有限差分法[15]在XPP 黏彈性流動問題的數(shù)值模擬方面已取得了一些研究成果.早在20 世紀(jì)70 年代，國外流變學(xué)者[19-20]就已經(jīng)采用有限差分法數(shù)值研究了黏彈性流動問題.有限體積法是在有限差分法的基礎(chǔ)上逐漸發(fā)展起來的數(shù)值計(jì)算方法.和有限元法相比，有限體積法對于守恒型方程的離散方程可以保持守恒，并且占用計(jì)算機(jī)內(nèi)存小、處理方程源項(xiàng)方便、數(shù)值穩(wěn)定性好[21].近年來，隨著流變學(xué)理論的成熟和計(jì)算機(jī)的快速發(fā)展，有限體積法不僅在黏性流體力學(xué)、電磁學(xué)、聲學(xué)以及傳熱學(xué)計(jì)算領(lǐng)域得到廣泛的應(yīng)用和推廣，而且在Oldroyd-B 黏彈性流動問題的數(shù)值模擬方面也顯示出了獨(dú)特的優(yōu)勢[22-24].

本文針對非等溫不可壓XPP 黏彈性流體流動特征建立控制方程，針對動量方程和XPP 本構(gòu)方程的對流項(xiàng)，采用基于延遲修正格式的高階AVLsmart 格式.另外，為了保證黏彈性流動迭代的快速收斂，采用基于壓力的同位網(wǎng)格有限體積(CLEAR)算法求解了Newton 黏性和非Newton XPP 黏彈性流動控制方程.最后，通過圓柱繞流問題進(jìn)行模擬，驗(yàn)證了本文數(shù)學(xué)模型及數(shù)值方法的有效性和準(zhǔn)確性.

1 數(shù)學(xué)模型和數(shù)值方法

1.1 數(shù)學(xué)模型

二維非等溫不可壓縮黏彈性流體的控制方程由連續(xù)性方程、動量方程、能量方程和本構(gòu)方程組成.連續(xù)性方程、動量方程和能量方程可表示為

1.2 數(shù)值方法

無量綱形式的控制方程(8) ～ (11)可以寫成如下統(tǒng)一的守恒形式：

在直角坐標(biāo)系中，Φ(可表示u， τ 和T)為待求解變量， θ 和 δ 為 常數(shù)，SΦ為依賴于所考察方程的源項(xiàng)(見表1). 本文采用CLEAR 法對流場、溫度場和應(yīng)力場的控制方程進(jìn)行離散，將所有求解變量都存儲于控制單元的中心.

表1 方程(18)中函數(shù)和常數(shù)的表達(dá)式Table 1 Definition of functions and constants in eq. (18)

為了提高數(shù)值計(jì)算的精度，本文采用高精度AVLsmart 格式[25]來離散動量和本構(gòu)方程的對流項(xiàng)，文獻(xiàn)[25]采用此格式求解了可壓縮黏彈性流動問題.為了保證待求解代數(shù)方程組滿足對角占優(yōu)的條件，增加動量方程和本構(gòu)方程的對流項(xiàng)高精度格式計(jì)算的穩(wěn)定性，本文采用延遲修正的方法在源項(xiàng)中加入高階AVLsmart 格式的修正項(xiàng).為了增強(qiáng)數(shù)值計(jì)算的穩(wěn)定性，提高迭代求解過程的收斂速度，本文對控制方程的離散方程組采用CLEAR 全隱式算法[26]進(jìn)行求解.充分保證了速度與壓力以及速度和應(yīng)力的耦合，解決了SIMPLE (semiimplicit method for pressure linked equation)算法中的第二個基本假設(shè)（壓力修正方程的推導(dǎo)中忽略了鄰點(diǎn)速度修正的影響）.

CLEAR 算法中每一層改進(jìn)后的速度和壓力直接由動量方程和壓力方程得出，不需要引入修正值.CLEAR 算法求解二維非等溫黏彈性流動問題的實(shí)施步驟如下：

⑥ 判斷是否滿足收斂條件，若不滿足，則返回步驟②繼續(xù)執(zhí)行下一時間步的迭代.

⑦ 求解本構(gòu)方程以及能量方程.

⑧ 重復(fù)步驟②～⑦，直到所有變量都滿足收斂條件.

上述迭代求解過程中，引入Φ(可表示u， τ 和T)的L2范數(shù)誤差：

其中，k和k?1 分別表示所求變量相鄰兩次的迭代次數(shù)，當(dāng)?shù)鷿M足收斂條件E(Φ)≤10?5時，迭代終止.

2 數(shù) 值 算 例

2.1 Newton 黏性流體圓柱繞流

圖1 Newton 流體圓柱繞流示意圖Fig. 1 Schematic diagram of the Newtonian flow past a cylinder

區(qū)域擴(kuò)充法是一種求解不規(guī)則區(qū)域流動和傳熱問題的有效方法[21].在實(shí)際計(jì)算過程中，將不規(guī)則區(qū)域擴(kuò)充為規(guī)則的計(jì)算區(qū)域（矩形區(qū)域），使得算法在程序中易于實(shí)施.以二維圓柱繞流問題為例，將擴(kuò)充區(qū)域（圓面內(nèi)）流體的動力黏度設(shè)為非常大的數(shù)（如1 020），這樣做相當(dāng)于把擴(kuò)充區(qū)域看作是黏度無限大的流體，使得該區(qū)域內(nèi)“流體的速度”比其他區(qū)域流體速度小許多數(shù)量級（擴(kuò)充區(qū)域內(nèi)流體幾乎不流動）.

圖2 給出了Re=10，20 和40 時，流場達(dá)到穩(wěn)態(tài)時圓柱附近的壓力和流線分布.由圖2 可知，在3 種Re下，壓力等值線光滑且呈上下對稱分布；對流線而言，邊界層發(fā)生分離，在圓柱的下游有一對穩(wěn)定的對稱漩渦，漩渦的中心逐漸遠(yuǎn)離圓柱表面，漩渦尺寸隨著Re的增大而不斷拉長，并向下游發(fā)展，流動為定常流動.

圖2 不同Re 下，圓柱附近的壓力(上)和流線分布(下)Fig. 2 Distributions of pressure (top) and streamlines (bottom) near the cylinder for different Re numbers

圖3 分別給出了Re= 10，20 和40 的水平速度u和垂直速度v的等值線圖.如圖所示在各種Re下，水平速度u和垂直速度v都呈對稱分布，本文計(jì)算結(jié)果與文獻(xiàn)[27]非常吻合，充分驗(yàn)證了本文方法的有效性.

圖3 不同Re 下，圓柱繞流的速度u(上)和速度v(下)分布Fig. 3 Distributions of velocity u (top) and velocity v (bottom) near the cylinder for different Re numbers

圓柱繞流中圓柱邊界受力的計(jì)算問題是關(guān)鍵問題之一.下面采用應(yīng)力積分法計(jì)算圓柱邊界表面總應(yīng)力，由曳力系數(shù)的定義得到曳力系數(shù)的公式為

其中，D為圓柱直徑，上述所有的值均是基于當(dāng)前時刻計(jì)算獲得的.表2 給出了用外插方法計(jì)算出的低Reynolds 數(shù)下黏性流體圓柱繞流的曳力系數(shù)，從表中可以看出隨著Reynolds 數(shù)的增大曳力系數(shù)逐漸變小，本文計(jì)算結(jié)果和文獻(xiàn)[28-30]結(jié)果一致.

表2 不同Re 下的曳力系數(shù)Table 2 Drag coefficients at different Re numbers

圖4 給出?t=0.01，t=100 時，在較高的Reynolds 數(shù)下，圓柱繞流周圍的流線圖.從圖4 明顯可以看出當(dāng)Re=60時，圓柱下游呈現(xiàn)不對稱性，流動為非定長流動；當(dāng)Re=80，100 和200 時，尾流中出現(xiàn)周期交替脫落的旋渦，并最終演變?yōu)橹芷谛缘牧鲃?，圓柱后面形成交錯排列的Karman 渦街，尾流保持層流狀態(tài).盡管圓柱的幾何形狀十分簡單，但流動卻十分復(fù)雜，因此圓柱繞流常用來檢驗(yàn)數(shù)值方法的有效性，上述模擬結(jié)果和已有文獻(xiàn)[31]結(jié)果一致.因此，基于區(qū)域擴(kuò)充的CLEAR 算法能夠有效模擬黏性流體圓柱繞流問題.

圖4 t=100 時，不同 Reynolds 數(shù)下圓柱繞流的流線圖：(a) Re=60；(b) Re=80；(c) Re=100；(d) Re=200Fig. 4 Streamlines of the flow near the cylinder at t=100 for different Re numbers: (a) Re=60; (b) Re=80; (c) Re=100; (d) Re=200

為了形象地考察圓柱后面旋渦生成與脫落的過程，本文給出了當(dāng) ?t=0.01，Re=100 時，半個周期內(nèi)不同時刻的流線圖.如圖5 所示，在t=T/12 時，圓柱右上方有旋渦生成.隨著時間的推移，這個旋渦不斷擴(kuò)大并向圓柱右下方移動.在t=6T/12 時，圓柱右下方也有一個小旋渦生成，這一過程記為旋渦脫落的半周期.

圖5 Re=100 時，半個周期內(nèi)圓柱附近流線圖：(a) t=T/12；(b) t=2T/12；(c) t=3T/12；(d) t=4T/12；(e) t=5T/12；(f) t=6T/12Fig. 5 Streamlines of the flow near the cylinder in a half cycle for Re=100: (a) t=T/12; (b) t=2T/12; (c) t=3T/12; (d) t=4T/12; (e) t=5T/12; (f) t=6T/12

2.2 非Newton XPP 黏彈性流體圓柱繞流

圓柱繞流是模擬黏彈性流體的基準(zhǔn)算例之一.如圖6 所示，XPP 黏彈性流體自左向右在兩個平行平板間流動，圓柱位于平板間的對稱位置上.圓柱的半徑為R，兩個平板的間距為4R.平板的上游和下游長度均取10R，R為1 個單位長度.整個計(jì)算區(qū)域?yàn)??=[0,20] × [0,4]，其中圓柱的圓心位于(10,2)，半徑為1 個單位長度.應(yīng)力邊界采用文獻(xiàn)[15]的方法進(jìn)行處理，速度和壓力邊界條件分別為：

圖6 XPP 流體圓柱繞流示意圖Fig. 6 Schematic diagram of the XPP flow past a cylinder

① 入口處速度為拋物型速度，umax=1.0，v= 0；

③ 圓柱表面和固壁面上設(shè)置為無滑移邊界條件，即u=v=0.

圖7 分 別 給 出 了 當(dāng)We=1.0，Re=1.0 ， β=1/9， α=0.15，q=2 ， λ0b/λ0s=3時 的 速 度u,v， 應(yīng) 力 分 量τxx,τxy, τyy和拉伸量 λ的等值線圖.從圖7 可以看出XPP 流體圓柱繞流時，流體在圓柱周圍區(qū)域經(jīng)歷了十分復(fù)雜的變形，所有的應(yīng)力分量和拉伸量在圓柱周圍產(chǎn)生較大的梯度，因此高速度、高剪切和高拉伸現(xiàn)象均出現(xiàn)在圓柱周圍. τxx的最小值出現(xiàn)在圓柱前面，而局部最大值出現(xiàn)在圓柱的上下壁面附近，靠近圓柱的上下兩個管道壁面附近以及圓柱后面下流區(qū)域的位置.正好相反，第二法向應(yīng)力 τyy的最大值出現(xiàn)在圓柱前沿區(qū)域.由于圓柱周圍存在高的剪切率，因此剪切應(yīng)力 τxy極值的絕對值出現(xiàn)在圓柱周圍.特別地，分子拉伸量 λ在圓柱前沿區(qū)域形成一個邊界，在圓柱周圍和管道壁面λ 存在其局部極大值.另外，從圖7 可以明顯看出，法向應(yīng)力 τxx， τyy和拉升量λ 關(guān)于y=2 對稱，而 τxy關(guān)于y=2 反對稱.總之，本文黏彈性流體圓柱繞流的數(shù)值模擬能夠提供圓柱周圍光滑并且無振蕩的應(yīng)力和拉升量結(jié)果，說明了本文高精度數(shù)值模擬方法模擬黏彈性流體圓柱繞流問題的有效性.

圖7 We=1，Re=1 時，速度、應(yīng)力分量及拉升量的分布：(a) u；(b) v；(c) τxx；(d) τyy；(e) τxy；(f) λFig. 7 Distributions of velocity vectors, stress components and stretches at We=1, Re=1: (a) u; (b) v; (c) τxx; (d) τyy; (e) τxy; (f) λ

We數(shù)是黏彈性流體流動中表征黏彈效應(yīng)的一個重要參數(shù)，圖8 給出了當(dāng)y=2.0(即在y方向中軸線上)，Re=1.0， β =1/9， α =0.15，q=2， ε =1/3時 ，We對速度、拉升量、第一法向應(yīng)力和第二法應(yīng)力的影響.由圖8 可知，We對表征黏彈性流體宏觀信息的量產(chǎn)生極大的影響，在We較 小時，彈性表現(xiàn)得不夠明顯；隨著We的增大，彈性才得以表現(xiàn).這主要是由于We很小時，流體彈性與黏性相比要小的多，即此時黏彈性流體與Newton 流體表現(xiàn)出相似的性質(zhì)；隨著We增大，彈性表現(xiàn)的越來越顯著，流體的彈性得以充分表現(xiàn).從圖8 可以清晰地看出，圓柱附近的速度值隨We的 增大而變小，而拉伸量隨著We的增加明顯增大，中軸線上靠近圓柱處的拉升量和法向應(yīng)力的絕對值明顯大于其他區(qū)域.

圖8 當(dāng)y=2 時，We 對水平速度、拉升量及法向應(yīng)力的影響：(a) u；(b) λ；(c) τxx；(d) τyyFig. 8 Effects of We on horizontal velocities, stretches and normal stresses at y=2: (a) u; (b) λ; (c) τxx; (d) τyy

為進(jìn)一步驗(yàn)證本文方法求解強(qiáng)對流占優(yōu)溫度方程的有效性，下面繼續(xù)研究非等溫條件下XPP 流體圓柱繞流問題.聚合物熔體選用低密度聚乙烯(low density polyethylene, LDPE)，流體密度為ρ =780 kg·m?3，比熱容為c=2 540 J·kg?1·K?1， 熱傳導(dǎo)率為 κ =0.241 W·m?1·K?1. 支鏈聚合物L(fēng)DPE 在Tr=443時相應(yīng)的XPP 模型和Arrhenius 模型參數(shù)如表3 所示[9].非等溫圓柱繞流問題的計(jì)算區(qū)域如圖6 所示.速度和應(yīng)力邊界條件和等溫問題相同，入口處溫度為Tmelt=443， 上下固壁和圓柱表面溫度為Twall=313.

表3 LDPE 熔體在T r =443時的XPP 和Arrhenius 流變參數(shù)[9]Table 3 Rheology parameters for LDPE melt for the XPP and Arrhenius models at T r =443[9]

圖9 給出了不同Pe數(shù)下圓柱繞流區(qū)域的溫度分布圖. 從圖9(a)中可以看出當(dāng)Pe=100時，對流作用較弱，對流導(dǎo)致的熱傳遞比較弱，因此溫度在固壁附近和圓柱后面比較低.隨著Pe的增加，對流作用逐漸增加，對流導(dǎo)致的熱傳遞作用也逐漸增強(qiáng)，因此相應(yīng)區(qū)域溫度逐漸增大.當(dāng)Pe=100 000時(圖9(d))，壁面低溫區(qū)域基本消失，圓柱繞流區(qū)域的溫度分布和圖9(a)相差較大.

圖9 不同Pe 下的溫度分布：(a) Pe=100；(b) Pe=1 000；(c) Pe=10 000；(d) Pe=100 000Fig. 9 Distributions of temperatures at different Pe numbers: (a) Pe=100; (b) Pe=1 000; (c) Pe=10 000; (d) Pe=100 000

3 結(jié) 論

本文針對不可壓Newton 流體、XPP 黏彈性流體圓柱繞流問題，發(fā)展了一種耦合延遲修正技術(shù)的高階AVLsmart 格式，采用CLEAR 算法結(jié)合區(qū)域擴(kuò)充方法求解控制方程.本文研究工作能為聚合物及纖維增強(qiáng)聚合物加工成型過程的研究奠定堅(jiān)實(shí)的理論基礎(chǔ). 所得主要結(jié)論如下：

1）Newton 流體在穩(wěn)態(tài)情形下壓力、流線、速度以及曳力系數(shù)和已有結(jié)果一致.在較高的Re下，尾渦的變化趨勢以及半個周期內(nèi)圓柱后面旋渦生成與脫落過程都和其他研究結(jié)果相近.數(shù)值結(jié)果表明本文高精度方法能夠穩(wěn)定、高效、精確地模擬Newton 黏性流動問題.

2）XPP 黏彈性流體圓柱繞流中速度、應(yīng)力分量及拉升量等值線光滑無震蕩，且均在圓柱周圍達(dá)到最大值.圓柱附近的速度值隨We的 增大而變?。欢亢偷谝环ㄏ驊?yīng)力隨著We的增大而變大. 事實(shí)上，AVLsmart 格式求解XPP 流動既具有較高的數(shù)值精度，又兼具延遲修正方法求解方便、魯棒性強(qiáng)、數(shù)值穩(wěn)定性好等特點(diǎn).

3）Pe對溫度場的分布影響非常大，在低Pe下，高溫區(qū)位于入口附近, 管道上下壁面以及圓柱后方溫度較低，隨著Pe的增大，高溫區(qū)逐漸向低溫區(qū)域推進(jìn)，壁面低溫區(qū)域基本消失.

致謝 本文作者衷心感謝河南工程學(xué)院博士培育基金(DKJ2019013)對本文的資助.

參考文獻(xiàn)（ References ） ：

［1］ 楊旭, 梁英杰, 孫洪廣, 等. 空間分?jǐn)?shù)階非Newton流體本構(gòu)及圓管流動規(guī)律研究[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2018, 39（11）:1213-1226. （YANG Xu, LIANG Yingjie, SUN Hongguang, et al. A study on the constitutive relation and the flow of spatial fractional non-Newtonian fluid in circular pipes[J].Applied Mathematics and Mechanics, 2018, 39（11）:1213-1226.（in Chinese））

［2］ 吳其曄, 巫靜安. 高分子材料流變學(xué)[M]. 北京: 高等教育出版社, 2002. （WU Qiye, WU Jing’an.Polymer Rheology[M]. Beijing: Higher Education Press, 2002. （in Chinese））

［3］ 白羽, 方慧靈, 張艷. Oldroyd-B流體繞拉伸楔形體的非穩(wěn)態(tài)滑移流動與傳熱分析[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2022, 43（3）:272-280. （BAI Yu, FANG Huiling, ZHANG Yan. Unsteady slip flow and heat transfer analysis of Oldroyd-B fluid over the stretching wedge[J].Applied Mathematics and Mechanics, 2022, 43（3）: 272-280.（in Chinese））

［4］LIELENS G, KEUNINGS R, LEGAT V. The FENE-L and FENE-LS closure approximations to the kinetic theory of finitely extensible dumbbells[J].Journal of Non-Newtonian Fluid Mechanics, 1999, 87（2/3）: 179-196.

［5］HYON Y K, DU Q, LIU C. An enhanced macroscopic closure approximation to the micro-macro FENE models for polymeric materials[J].Multiscale Modeling Simulation, 2008, 7（2）: 978-1002.

［6］PHAN-THIEN N, TANNER R I. A new constitutive equation derived from network theory[J].Journal of Non-Newtonian Fluid Mechanics, 1977, 2（4）: 353-365.

［7］MCLEISH T C B, LARSON R G. Molecular constitutive equations for a class of branched polymers: the Pom-Pom polymer[J].Journal of Rheology, 1998, 42（1）: 81-110.

［8］VERBEETEN W M H, PETERS G W M, BAAIJENS F T P. Differential constitutive equations for polymer melts: the extended Pom-Pom model[J].Journal of Rheology, 2001, 45（4）: 823-843.

［9］VERBEETEN W M H, PETERS G W M, BAAIJENS F P T. Viscoelastic analysis of complex polymer melt flows using the eXtended Pom-Pom model[J].Journal of Non-Newtonian Fluid Mechanics, 2002, 108（1/3）: 301-326.

［10］VERBEETEN W M H, PETERS G W M, BAAIJENS F P T. Numerical simulations of the planar contraction flow for a polyethylene melt using the XPP model[J].Journal of Non-Newtonian Fluid Mechanics, 2004, 117（2/3）:73-84.

［11］TAMADDON-JAHROMI H R, WEBSTER M F. Transient behaviour of branched polymer melts through planar abrupt and rounded contractions using Pom-Pom models[J].Mechanics of Time-Dependent Materials, 2011, 15:181-211.

［12］AGUAYOJP, PHILLIPSPM, PHILLIPSTN, et al. Thenumerical prediction of planar viscoelastic flows using the Pom-Pom model and high-order finite volumeschemes[J].Journal of Computational Physics, 2007, 220（2）: 586-611.

［13］YANG B X, OUYANG J, LI Q, et al. Modeling and simulation of the viscoelastic fluid mold filling process by level set method[J].Journal of Non-Newtonian Fluid Mechanics, 2010, 165（19/20）: 1275-1293.

［14］LI X J, HE J H. Variational multi-scale finite element method for the two-phase flow of polymer melt filling process[J].International Journal of Numerical Methods for Heat & Fluid Flow, 2020, 30（3）: 1407-1426.

［15］OISHI C M, MARTINS F P, TOME M F, et al. Numerical solution of the eXtended Pom-Pom model for viscoelastic free surface flows[J].Journal of Non-Newtonian Fluid Mechanics, 2011, 166（3/4）: 165-179.

［16］ 曹偉, 王蕊, 申長雨. 塑料熔體在注塑模中的三維流動模擬[J]. 化工學(xué)報(bào), 2004, 55（9）: 1493-1498. （CAO Wei,WANG Rui, SHEN Changyu. Three dimensional flow simulation of plastic melt in injection mold[J].CIESC Journal, 2004, 55（9）: 1493-1498.（in Chinese））

［17］ 周文, 歐陽潔, 楊斌鑫, 等. 三維非等溫非牛頓流體充模過程的建模與模擬[J]. 化工學(xué)報(bào), 2011, 62（3）: 618-627.（ZHOU Wen, OUYANG Jie, YANG Binxin, et al. Modeling and simulation of 3D non-isothermal non-Newtonian fluid filling process[J].CIESC Journal, 2011, 62（3）: 618-627.（in Chinese））

［18］GAO P Y, ZHAO Z T, YANG Y. The numerical modeling and study of gas entrapment phenomenon in non-isothermal polymer filling process[J].Journal of Non-Newtonian Fluid Mechanics, 2021, 294: 104575.

［19］CROCHET M J, PILATE G. Plane flow of a fluid of second grade through a contraction[J].Journal of Non-Newtonian Fluid Mechanics, 1976, 1（3）: 247-258.

［20］PERERA M G N, WALTERS K. Long-range memory effects in flows involving abrupt changes in geometry, partⅠ: flows associated with L-shaped and T-shaped geometries[J].Journal of Non-Newtonian Fluid Mechanics,1977, 2（1）: 49-81.

［21］ 陶文銓. 數(shù)值傳熱學(xué)[M]. 2版. 西安: 西安交通大學(xué)出版社, 2001. （TAO Wenquan.Numerical Heat Transfer[M].2nd ed. Xi’an: Xi’an Jiaotong University Press, 2001.（in Chinese））

［22］HABLA F, TAN M W, HABLBERGER J, et al. Numerical simulation of the viscoelastic flow in a three-dimensional lid-driven cavity using the log-conformation reformulation in OpenFOAM[J].Journal of Non-Newtonian Fluid Mechanics, 2014, 212: 47-62.

［23］PIMENTA F, ALVES M A. Stabilization of an open-source finite-volume solver for viscoelastic fluid flows[J].Journal of Non-Newtonian Fluid Mechanics, 2017, 239: 85-104.

［24］ZHOU W, OUYANG J, WANG X D, et al. Numerical simulation of viscoelastic fluid flows using a robust FVM framework on triangular grid[J].Journal of Non-Newtonian Fluid Mechanics, 2016, 236: 18-34.

［25］ZHUANG X, OUYANG J, LI Y G, et al. A three-dimensional thermal model for viscoelastic polymer melt packing process in injection molding[J].Applied Thermal Engineering, 2018, 128: 1391-1403.

［26］TAO W Q, QU Z G, HE Y L. A novel segregated algorithm for incompressible fluid flow and heat transfer problems-CLEAR (coupled and linked equations algorithm revised) part Ⅰ: mathematical formulation and solution procedure[J].Numerical Heat Transfer(Part B):Fundamentals, 2004, 45（1）: 1-17.

［27］GAO P Y, OUYANG J, DAI P F, et al. A coupled continuous and discontinuous finite element method for the incompressible flows[J].International Journal for Numerical Methods in Fluids, 2017, 84（8）: 477-493.

［28］DING H, SHU C, YEO K S, et al. Simulation of incompressible viscous flow past a circular cylinder by hybrid FD scheme and meshless least squarebased finite difference method[J].Computer Methods in Applied Mechanics and Engineering, 2004, 193（9/11）: 727-744.

［29］NITHIARASU P, ZIENKIEWICZ O C. Analysis of an explicit and matrix free fractional stepmethod for incompressible flows[J].Computer Methods in Applied Mechanics and Engineering, 2006, 195（41/43）: 5537-5551.

［30］TUANN S Y, OLSON M D. Numerical studies of the flow around a circular cylinder by a finite element method[J].Computers and Fluids, 1978, 6（4）: 219-240.

［31］ 仇軼, 由長福, 祁海鷹, 等. 用無網(wǎng)格法求解不同Re下圓柱繞流問題[J]. 清華大學(xué)學(xué)報(bào)(自然科學(xué)版), 2005, 45（2）:220-223. （QIU Yi, YOU Changfu, QI Haiying, et al. Numerical simulation of flow around a cylinder for differentReusing a meshless method[J].Journal Tsinghua University(Science and Technology), 2005, 45（2）: 220-223.（in Chinese））