非等温黏弹性聚合物流体圆柱绕流的高精度数值模拟*

庄 昕， 刘付军， 孙艳萍， 王惠玲

（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））