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

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

（1. 河南工程學院 理學院，鄭州 451191；2. 山西財經大學 應用數學學院，太原 030006）

引 言

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

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

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

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

1 數學模型和數值方法

1.1 數學模型

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

1.2 數值方法

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

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

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

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

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

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

⑦ 求解本構方程以及能量方程.

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

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

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

2 數 值 算 例

2.1 Newton 黏性流體圓柱繞流

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

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

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

圖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都呈對稱分布，本文計算結果與文獻[27]非常吻合，充分驗證了本文方法的有效性.

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

圓柱繞流中圓柱邊界受力的計算問題是關鍵問題之一.下面采用應力積分法計算圓柱邊界表面總應力，由曳力系數的定義得到曳力系數的公式為

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

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

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

圖4 t=100 時，不同 Reynolds 數下圓柱繞流的流線圖：(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

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

圖5 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/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 黏彈性流體圓柱繞流

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

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

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

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

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

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

圖8 當y=2 時，We 對水平速度、拉升量及法向應力的影響：(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

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

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

圖9 給出了不同Pe數下圓柱繞流區域的溫度分布圖. 從圖9(a)中可以看出當Pe=100時，對流作用較弱，對流導致的熱傳遞比較弱，因此溫度在固壁附近和圓柱后面比較低.隨著Pe的增加，對流作用逐漸增加，對流導致的熱傳遞作用也逐漸增強，因此相應區域溫度逐漸增大.當Pe=100 000時(圖9(d))，壁面低溫區域基本消失，圓柱繞流區域的溫度分布和圖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 結 論

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

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

2）XPP 黏彈性流體圓柱繞流中速度、應力分量及拉升量等值線光滑無震蕩，且均在圓柱周圍達到最大值.圓柱附近的速度值隨We的 增大而變小；而拉升量和第一法向應力隨著We的增大而變大. 事實上，AVLsmart 格式求解XPP 流動既具有較高的數值精度，又兼具延遲修正方法求解方便、魯棒性強、數值穩定性好等特點.

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

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

參考文獻（ References ） ：

［1］ 楊旭, 梁英杰, 孫洪廣, 等. 空間分數階非Newton流體本構及圓管流動規律研究[J]. 應用數學和力學, 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］ 吳其曄, 巫靜安. 高分子材料流變學[M]. 北京: 高等教育出版社, 2002. （WU Qiye, WU Jing’an.Polymer Rheology[M]. Beijing: Higher Education Press, 2002. （in Chinese））

［3］ 白羽, 方慧靈, 張艷. Oldroyd-B流體繞拉伸楔形體的非穩態滑移流動與傳熱分析[J]. 應用數學和力學, 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]. 化工學報, 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]. 化工學報, 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］ 陶文銓. 數值傳熱學[M]. 2版. 西安: 西安交通大學出版社, 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］ 仇軼, 由長福, 祁海鷹, 等. 用無網格法求解不同Re下圓柱繞流問題[J]. 清華大學學報(自然科學版), 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））