用无网格法求解不同Re下圆柱绕流问题_仇轶

ISSN 100020054 CN 1122223öN 　 清华大学学报 (自然科学版)J T singhua U niv (Sci & Tech) , 2005 年 第 45 卷 第 2 期2005, V o l. 45, N o. 2 20ö362202223 　 用无网格法求解不同 R e 下圆柱绕流问 仇　轶, 　由长福, 　祁海鹰, 　徐旭常 (清华大学 热能工程系, 煤的清洁燃烧技术国家重点实验室, 北京 100084) 收稿日期: 2004202227 基金项目: 国家自然科学基金资助项目 (10102009) 作者简介: 仇轶 (19782) , 男 (汉) , 湖南, 博士研究生。 通讯联系人: 由长福, 副教授, E2m ail: youcf@ te. tsinghua. edu. cn 摘　要: 采用无网格伽辽金法对不同 R e 下二维不可压粘性 流场圆柱绕流问题进行了数值模拟。该计算不仅预测了不同 条件下的流动特性, 而且预测了流体作用在圆柱上的曳力和 升力。这些力通过积分圆柱表面的压力和粘性力获得。结果 显示, 当R e 大于 80 时, 圆柱所受曳力和升力开始振荡。计算 所得的脱落频率 (S r)与前人的结果吻合良好。 关键词: 流体力学; 不可压缩粘性流体; 无网格 Galerk in 法 (EFG) ; 分步法 中图分类号: O 357. 1 文献标识码: A 文章编号: 100020054 (2005) 0220220204 Num er ica l s im ula tion of f low around a cyl inder for d ifferen t R e using a m eshless m ethod Q IU Yi1, YOU Cha ngfu, Q I Ha iy ing , XU Xucha ng (State Key Laboratory of Clean Combustion of Coal, D epartmen t of Thermal Eng ineer ing, Tsinghua Un iversity, Be ij ing 100084, China) Abstract: T he elem en t free Galerk in m ethod, EFG, w as used to sim ulate viscous, incomp ressib le 2D flow s around a cylinder fo r differen t Reyno lds num bers. T h is analysis p redicted no t on ly the flow patterns fo r differen t conditions bu t also the drag and lift fo rces exerted by the flu id on the cylinder as w ell. T hese fo rces w ere calcu lated by integrating the p ressure and shear fo rce over the cylinder surface. T he resu lts show that the drag and lift fo rces o scillate fo r Reyno lds num bers greater than 80. T he calculated Strouhal num ber, a param eter characterizing the o scillat ions, agrees w ell w ith established resu lts. Key words: f lu id m echan ics; viscous incomp ressib le flow s; elem ent free Galerk in (EFG) m ethod; fractional step m ethod 高 R e 流动计算关键问题之一是解决由于强对 流作用带来的数值波动问题。目前公认精度较高的 有流线迎风格式 (st ream line upw ind2Petrov2 Galerk in, SU PG) [1 ] , T aylo r2Galerk in (T 2G) 格式[2 ] 等。T 2G 格式对于线性对流方程具有三阶精度, 但 该格式有高阶空间导数项, 应用于非线性问题比较 困难。文[ 3 ]提出了一种三阶有限元分步格式, 没有 新的高阶导数项出现, 适用于非线性问题的求解。 无网格的目的在于消除一部分或者全部的计算 网格, 而完全基于点的近似, 不需要将结点连成单 元。无网格方法目前主要应用于固体力学领域[4 ] , 在 流体力学研究领域该方法一直没有得到很大的发 展, 应用也比较少。 本 文采用无 网 格 伽 辽 金 法 ( elem en t free Galerk in m ethod, EFG) [4, 5 ]计算不同 R e 下二维粘 性不可压流场圆柱绕流问题。计算中采用文[ 3 ]提出 的速度压力分离的格式, 能够模拟出圆柱体尾流区 中交错排列的 Karm an 涡街, 且得到的涡街脱落的 频率与前人的计算结果吻合得较好。关于EFG 基本 原理及数值离散化的实现可参考文[ 4 6 ]。 1　控制方程及离散化的实现 二维不可压粘性流动控制方程为9u i9 t + u j u i, j = - p , iΘ + v (u i, j + u j , i) , j + f i, u i, i = 0. (1) 其中 f i 为体积力。 采用文[ 3 ]提出的三步有限元格式, 可得到N 2S 方程的时间离散格式: u n+ 1ö3 i - u n i∃ tö3 = - unj u ni, j - p n, iΘ +Μ(u ni, j + u nj , i) , j + f ni , u n+ 1ö2 i - u n i∃ tö3 = - un+ 1ö3j un+ 1ö3i, j - p n, iΘ + Μ(un+ 1ö3i, j + un+ 1ö3j , i ) , j + f n+ 1ö3i , u n+ 1 i - u n i∃ tö3 = - un+ 1ö2j un+ 1ö2i, j - p n+ 1, iΘ +Μ(un+ 1ö2i, j + un+ 1ö2j , i ) , j + f n+ 1ö2i . (2) 采用的 Galerk in 法对方程 (2) 进行空间离散, 可得最终的无网格方程如下: M ΑΒ un+ 1ö3iΒ - uniΒ∃ tö3 = - N nΑΒuniΒ - L iΑΒp nΒΘ - ΜS niΑ + M ΑΒf niΒ +∫# 5 ΑΜ(uni, j + unj , i) n id# , M ΑΒ un+ 1ö2iΒ - uniΒ∃ tö2 = - N n+ 1ö3ΑΒ un+ 1ö3iΒ - L iΑΒp nΒΘ - ΜS n+ 1ö3iΑ + M ΑΒf n+ 1ö3iΒ + ∫# 5 ΑΜ(un+ 1ö3i, j + un+ 1ö3j , i ) n id# , M ΑΒ un+ 1iΒ - uniΒ∃ t = - N n+ 1ö2ΑΒ u n+ 1ö2iΒ - L iΑΒp n+ 1ΒΘ - ΜS n+ 1ö2iΑ + M ΑΒf n+ 1ö2iΒ + ∫# 5 ΑΜ(un+ 1ö2i, j + un+ 1ö2j , i ) n id#. (3) 其中: M ΑΒ = κ8 e 5 Α5 Βd8 , L iΑΒ = κ8 e 5 Α5 Β, id8 , N nΑΒ = κ8 e 5 Α5 Μ5 Β, ju njΜd8 , S niΑ = κ8 e 5 Α, i (5 Β, ju niΒ + 5 Β, iu njΒ) d8. (4)5 Α是形函数。由方程 (3)可知, 在求解 un+ 1i 前必须先 求出 p n+ 1。对方程 (3) 最后一式两边同时取散度, 并 引入不可压条件 un+ 1i, i = 0, 可导出压力 Po isson 方 程。同样, 对该方程采用 Galerk in 离散, 可得到最终 的无网格压力方程 G ΑΒp n+ 1ΒΘ = - 1∃ tH iΑΒu niΒ - K jΑΒun+ 1ö2jΒ - S Α, (5) 其中: G ΑΒ = κ8 e 5 Α, i5 Β, id8 , L iΑΒ = κ8 e 5 Α5 Β, id8 , K jΑΒ = κ8 e 5 Α, i5 Β5 Μ, j un+ 1ö2iΜ d8 ,S Α =∫# 5 Α un+ 1i - uni∃ t d#. (6)因此, 由式 (5)求出 p n+ 1后就能由方程 (3)求出 un+ 1i 。这种三阶有限元分步格式数值稳定性好, 并且没有引入人工参数, 适合用于高 R e 下不可压流动的数值模拟。值得注意的是, 本文采用L agrange 乘子法引入速度和压力本质边界条件, 具体方法请参看文[ 5 ]。2　计算结果及讨论为了将无网格法的计算结果和前人的结果进行比较, 本文采用了文[ 7 ]的算例。计算区域见图 1。图 1　二维不可压粘性流体圆柱绕流计算区域取D = 1m , 流体密度为 Θ= 1. 0 kg öm 3。在管道入口处, x 方向流速呈抛物线分布, 形式如下:u x = 6UH 2 [ (y - yB )H - (y - y B ) 2 ],u y = 0. (7)其中: U = 1. 0m ös, 为入口平均流速; H = 4. 1D ,为管道的高度。以圆柱圆心为中心, 向周围辐射式分布 4 907 个结点, 本文分别对 R ep = 20、100、200、1 000 四种情况下的流场分布进行了数值模拟。得到的计算结果如图 2 所示。由图 2 可知, 在R e 较小时 (R e= 20, t= 100 s) ,流动在圆柱背后形成对称的旋涡, 此时为定常流动。随着 R e 增加, 圆柱后面的旋涡逐渐趋于非对称, 当R e 超过某个临界值时, 流动变为非定常流动, 而且涡开始从圆柱上脱落, 并最终演变为周期性的流动。如图 2 所示, 在R e= 100 时, 圆柱后面形成交错排列的 Karm an 涡街。在 R e= 1 000 时, 对流作用强, 流态变化复杂, 用本文的方法也能得到稳定的数值解。由于圆柱上有旋涡脱离的现象, 圆柱受到的流体作用的曳力和升力处于脉动状态, 定义曳力系数 122仇　轶, 等: 　用无网格法求解不同 R e 下圆柱绕流问题 和升力系数分别为: C d = f xΘU 2ö2, C l = f yΘU 2ö2. (8) 其中 f x 和 f y 分别为圆柱体所受的 x 方向和 y 方向 的作用力, 由积分圆柱表面所受的压力和粘性力 求得。 图 2　不同Re 下圆柱周围的流线 图 3 所示为图 1 计算域下R e= 100 时曳力系数 和升力系数随时间变化曲线。如果流场计算域是完 全对称的话, 那得到的升力系数应该在 0 左右脉动。 而本文采用的计算域是非对称的, 相应的升力系数 则在C 1,m in= - 1. 004 和C 1,m ax= 1. 085 之间脉动, 曳 力系数则在 C d,m in = 3. 111 和 C d,m ax = 3. 180 之间 脉动。 由图 3 还可以看出, 曳力系数脉动频率为升力 系数脉动频率的两倍。这是由于在每一个旋涡的生 长和脱落过程中, 曳力总有一个最大值和最小值, 而 升力的变化则取决于这个涡是位于圆柱的上方还是 下方。定义 S r = DU T . (9) 式中 T 为升力系数 C l 的振荡周期。由图 3 可得, R e= 100 时, C l 的振荡周期约为 3. 4 s。相应的 S r 为 0. 294, 与文[ 7 ]中的结果 0. 301 8 非常接近。 图 3　Re= 100 时, 曳力系数和升力系数 随时间变化曲线 上述采用的算例中, 圆柱受壁面影响较大。本文还计 算了圆柱在不受壁面影响的情况下, S r 为 0. 204, 与文[ 7 ]中 0. 18～ 0. 2 的结果也相当接近。因为在管 道中壁面的限制加速了旋涡脱离的过程, 从而导致 圆柱在受壁面影响条件下的 S r 比不受壁面影响条 件下的 S r 大得多。 为了形象地考察圆柱后面旋涡生成与脱落的过 程, 本文计算了R e= 100 的流动。图 4 给出了不同时 刻圆柱附近的流线。如图所示, 在 t= 78. 2 s 时, 圆柱 右上方有旋涡生成。随着时间的推移, 这个旋涡不断 扩大并向圆柱右下方移动。在 t= 79. 8 s 时, 圆柱右 下方也有一个小旋涡生成。这一过程记为旋涡脱落 的半周期, 相应的 S r 为 0. 312 5, 与文[ 7 ]的结果很 接近。 3　结　论 本文在时间上用三阶分步有限元法, 在空间上 用无网格Ga lerk in 法进行离散, 对二维圆柱绕流进 行了模拟。结果表明, 在R e 较小时, 圆柱后面形成 对称的旋涡; 随着R e 的增加, 圆柱背后的旋涡逐渐 趋于非对称, 并形成交错排列的 K a rm an 涡街。此 外, 本文还对圆柱所受流体作用力进行了研究, 由计 算得到的作用力脉动规律与文献结果符合得很好, 表明本文所构造的无网格离散格式可用于流体流动 的数值模拟, 丰富了流体力学计算方法。 222 清 华 大 学 学 报 (自 然 科 学 版) 2005, 45 (2) 图 4　不同时刻圆柱附近流线 (Re= 100) 参考文献　 (References) [ 1 ] T ezduyar T E, Ganjoo D K. Petrov2Galerk in fo rm ulations w ith w eigh ting funct ions dependent upon spatial and tempo ral discretion: app lication to transient convect ion2diffusion p roblem s [J ]. Com p u t M ethod s A p p l M ech E ng , 1985, 59: 49 71. [ 2 ] Donea J , Q uartapelle L , Selm in V. A n analysis of tim e discretization in the fin ite elem en t so lution of hyperbo lic p roblem s [J ]. J Com p u t P hy s, 1978, 70: 463 499. [ 3 ] 江春波, 徐照明, 李秀丽. 求解不可压缩流动的分步有限元 格式 [J ]. 清华大学学报 (自然科学版) , 2002, 42 (2) : 278 280. J IAN G Chunbo, XU Zhaom ing, L I X iuli. F ractional step fin ite elem ent fo rm ulation fo r so lving incomp ressib le flow s [J ]. J T sing hua U n iv (S ci & T ech) , 2002, 42 (2) : 278 280. ( in Chinese) [ 4 ] 宋康祖, 陆明万, 张雄. 固体力学中的无网格方法 [J ]. 力学 进展, 2000, 30 (1) : 55 65. SON G Kangzu, LU M ingw an, ZHAN G X iong. M esh less m ethod fo r so lid m echan ics [J ]. A d vances in M echan ics, 2000, 30 (1) : 55 65. ( in Ch inese) [ 5 ] Belytschko T , Gu L , L iu Y Y. F racture and crack grow th by elem en t2free Galerk in m ethods [J ]. M odel S im u l M aster S ci E ng , 1994, 2: 519 534. [ 6 ] Q iu Y. Imp rovem ent of search ing nodes arithm etic in m eshfree m ethod [A ]. P roceedings of the A SM EöJSM E Jo in t F luids Engineering Conference [C ]. Hono lulu, H I, U SA : Am erican Society of M echan ical Engineers, 2003. 2023 2026. [ 7 ] Joel H F, M ilovan P. Computational M ethods fo r F lu id D ynam ics [M ]. Germ any: Sp ringer, 2002. 322仇　轶, 等: 　用无网格法求解不同 R e 下圆柱绕流问题