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

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

引 言

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

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

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

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

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

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

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

1.2 數(shù)值方法

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

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

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

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

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

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

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

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

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

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

2 數(shù) 值 算 例

2.1 Newton 黏性流體圓柱繞流

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

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

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

圖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都呈對(duì)稱分布，本文計(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ì)算問(wèn)題是關(guān)鍵問(wèn)題之一.下面采用應(yīng)力積分法計(jì)算圓柱邊界表面總應(yīng)力，由曳力系數(shù)的定義得到曳力系數(shù)的公式為

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

圖4 t=100 時(shí)，不同 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

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

圖5 Re=100 時(shí)，半個(gè)周期內(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 黏彈性流體自左向右在兩個(gè)平行平板間流動(dòng)，圓柱位于平板間的對(duì)稱位置上.圓柱的半徑為R，兩個(gè)平板的間距為4R.平板的上游和下游長(zhǎng)度均取10R，R為1 個(gè)單位長(zhǎng)度.整個(gè)計(jì)算區(qū)域?yàn)??=[0,20] × [0,4]，其中圓柱的圓心位于(10,2)，半徑為1 個(gè)單位長(zhǎng)度.應(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è)置為無(wú)滑移邊界條件，即u=v=0.

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

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

圖8 當(dāng)y=2 時(shí)，We 對(duì)水平速度、拉升量及法向應(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)對(duì)流占優(yōu)溫度方程的有效性，下面繼續(xù)研究非等溫條件下XPP 流體圓柱繞流問(wèn)題.聚合物熔體選用低密度聚乙烯(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時(shí)相應(yīng)的XPP 模型和Arrhenius 模型參數(shù)如表3 所示[9].非等溫圓柱繞流問(wèn)題的計(jì)算區(qū)域如圖6 所示.速度和應(yīng)力邊界條件和等溫問(wèn)題相同，入口處溫度為T(mén)melt=443， 上下固壁和圓柱表面溫度為T(mén)wall=313.

表3 LDPE 熔體在T r =443時(shí)的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時(shí)，對(duì)流作用較弱，對(duì)流導(dǎo)致的熱傳遞比較弱，因此溫度在固壁附近和圓柱后面比較低.隨著Pe的增加，對(duì)流作用逐漸增加，對(duì)流導(dǎo)致的熱傳遞作用也逐漸增強(qiáng)，因此相應(yīng)區(qū)域溫度逐漸增大.當(dāng)Pe=100 000時(shí)(圖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é) 論

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

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

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

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

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

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

［1］ 楊旭, 梁英杰, 孫洪廣, 等. 空間分?jǐn)?shù)階非Newton流體本構(gòu)及圓管流動(dòng)規(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)滑移流動(dòng)與傳熱分析[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］ 曹偉, 王蕊, 申長(zhǎng)雨. 塑料熔體在注塑模中的三維流動(dòng)模擬[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］ 周文, 歐陽(yáng)潔, 楊斌鑫, 等. 三維非等溫非牛頓流體充模過(guò)程的建模與模擬[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］ 仇軼, 由長(zhǎng)福, 祁海鷹, 等. 用無(wú)網(wǎng)格法求解不同Re下圓柱繞流問(wèn)題[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））