基于区域分解, 有限差分法和涡法杂交的数值方法, 结合高阶隐式差分格式, 和以修正的不完全LU分解为预处理器的共轭梯度法作求解器. 研究了雷诺数Re＝1000, 旋转速度比α∈(0.5,3.25)范围内, 绕旋转圆柱从突然起动到充分发展, 长时间内尾流旋涡结构和阻力, 升力系数的变化规律. 计算所得流动图案与实验流场显示符合很好, 数值试验证实了临界状态的存在

凌国平

凌国灿

Institute Of Mechanics,Chinese Academy of Sciences

第　29　卷　第　1　期1997　年　1　月　　　力 学 学 报A CTA M ECHAN ICA S IN ICAV o l. 29,N o. 1　　Jan. , 19971)江苏省自然科学基金和中国科学院力研究所LNM 开放实验室资助课题.1990- 02- 09 收到第一稿, 1996- 07- 30 收到修改稿.绕旋转圆柱流动涡尾流结构和临界状态特性1)凌国平(苏州大学数学系, 苏州 215006)凌国灿(中国科学院力学研究所, 北京 100080)摘要　采用作者提出的基于区域分解、有限差分法与涡法杂交的数值方法, 结合高阶隐式差分格式, 和以修正的不完全LU 分解为预处理器的共轭梯度法作求解器. 系统地研究了雷诺数 R e= 1 000, 旋转速度比 a∈ (0. 5, 3. 25) 范围内, 绕旋转圆柱从突然起动到充分发展, 长时间内尾流旋涡结构和阻力、升力系数的变化规律. 计算所得流动图案与实验流场显示符合很好. 数值试验证实了临界状态的存在, 并首次给出了临界状态时的旋涡结构特性, 指出最佳升阻比也在该状态附近得到.关键词　旋转圆柱, 涡尾流, 有限差分法, 涡法, 预处理共轭梯度法, 不完全LU 分解引　言绕旋转圆柱的流动是一种典型的非定常分离流动. 它在改变旋涡的释放规律, 提高升力, 减小阻力, 降低流动引起的振动等一系列现代流动控制研究中, 有着重要的理论和实际应用价值,一直受到流体力学家们的重视.圆柱旋转产生的一个最显著的效应, 就是减弱和抑制在圆柱某一边流动的分离和旋涡的释放, 同时加强和发展在圆柱另一边流动的分离和旋涡的释放, 从而使整个流动在垂直来流方向成为不对称的, 将有一个升力作用在圆柱上, 这就是著名的M agnu s 效应. 影响流动的参数除了雷诺数 R e (= 2U ∞aöv , U ∞为来流速度, a 为圆柱半径, v 为流体的运动粘性系数) 外, 最重要的是旋转速度比 a (= 8 aöU ∞, 8 为圆柱旋转角速度). 　早期的研究集中在小 a 和低 R e 的范围. 80 年代中、后期开始对较大 a 和较高 R e 范围的流动进行了一系列实验[ 1～ 3 ]和数值[ 3～ 6 ]的研究, 但也只限于突然起动后的初始阶段 ( t≤ 24) , 和柱后的近尾流区, 而对长时间的充分发展流动和柱后涡尾流的整体结构特性仍缺乏研究. 多数数值结果与实验结果相符, 但也有例外, 其中最突出的是文[ 5 ]. 它指出, 在 R e= 200, a = 3. 25 时, 有多于一个主涡的所谓单边释放涡现象发生, 这是与文[1 ]的实验结果相矛盾的. 文[1 ]实验观察指出, 对如此大的 a , 在第一个主涡释放后, 流动即趋于稳定, 以后再无涡释放. 最近文[7 ]对 R e= 1 000, 较大 a 范围内绕旋转圆柱长时间的流动特性进行了系统的数值研究, 支持了文[1 ]的结论, 但未对流动从周期性双边涡结构过渡到稳定的单边涡结构之间临界状态的流动特性作进一步的探讨.© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd.   All rights reserved.　　本文采用作者 (1992) 提出的基于区域分解、有限差分法与涡法杂交的数值方法[ 8～ 10 ] , 并作了改进: 差分离散采用高阶格式, 时间步全隐式推进, 由此产生一个宽带大型稀疏矩阵方程, 采用目前认为属于最优一类的以修正的不完全 LU 分解作预处理器的共轭梯度法 (M ILU -CG) [ 11, 12 ]进行求解. 系统地研究了 R e= 1 000, a∈ (0. 5, 3. 25) 范围内, 绕旋转圆柱从突然起动到充分发展, 长时间内整个尾流旋涡结构和阻力、升力系数的变化规律. 计算所得典型的流动图案与实验流场显示符合很好. 计算证实了临界状态的存在, 并首次指出此时的涡尾流呈藕节状单边释放涡的结构, 阻力、升力呈微幅周期波动. 最后还讨论了控制优化问题, 指出最佳升阻比就发生在临界状态附近.1　控制方程和定解条件设不可压、粘性流体由无穷远处以均匀速度U ∞, 沿 X 正方向流向半径为 a , 且以角速度 8按逆时针方向旋转的无限长圆柱 (如图 1 所示).图 1　流动示意图F ig. 1　F low difin it ion sketch　　取原点在圆柱轴心, 但不随圆柱一起旋转的极坐标系 (r, Η). 　在该坐标系下, 以涡量 Ξ- 流函数 7 的形式表示, 并将变量 t, r, (V r,V Η) , 7 , Ξ 分别用 aöU ∞, a ,U ∞,U ∞, a 和U ∞öa 无量纲化后的流动控制方程为 5Ξ5t + 1r 55r (rV rΞ) + 1r 55Η(V ΗΞ) = 2R eý 2Ξ (1)ý 27 = - Ξ (2)其中, (V r,V a)分别为径向与周向速度, 它们与流函数 7 的关系为V r = 1r575Η,V Η = - 575r为加密柱面附近的网格分布, 引入变换 r= exp (2ΠΝ) , Η= 2ΠΓ. 流动控制方程 (1)、(2)变成E 5Ξ5t + 55Ν(U Ξ) + 55Γ(V Ξ) = 2R e (52Ξ5Ν2 + 52Ξ5Γ2) (3)5275Ν2 + 5275Γ2 = - E Ξ (4)其中, E = 4Π2e4ΠΝ,U = 575Γ= E 1ö2V r,V = - 575Ν= E 1ö2V Η.9第　1　期 　 凌国平等: 绕旋转圆柱涡尾流结构和临界状态特性 　© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd.   All rights reserved.　　边界条件为, 在柱面上 (Ν= 0)满足7 = 0, 575Ν = - E 1ö2aΞ = - 1E 5275Ν2 (5)在无穷远处 (Ν→∞)满足 5757 = E 1ö2 õ sin (2ΠΓ)5Ξ5Ν= 0 (6)初始条件为 Ξû t= 0 = 0, Ν> 0 (7)由方程 (3)～ (7)求得流场中的涡量分布后, 就可推出沿圆柱表面的压力和切应力分布, 进而求出作用在圆柱上的阻力系数C d 和升力系数C ΙC d = 2R e10(5Ξ5Ν- 2ΠΞ) û Ν= 0 õ sin (2ΠΓ) d ΓΘC Ι = 2R e 10(2ΠΞ - 5Ξ5Ν) û Ν= 0 õ cos (2ΠΓ) d ΓΘ (8)2　数值求解方法数值求解的基本思想是作者 (1992) 提出的基于区域分解、有限差分法和涡法杂交的数值模型, 并在差分格式和求解方法上作了较大改进.将流场分解为靠近柱面范围为O (a)量级的内区 8 1 和余下的外区 8 2 两部分. 内区采用有限差分解, 外区采用网格涡解, 内外区流动通过交界面 I 互相耦合.对方程 (3)中的对流项采用三阶偏心格式, 以 55Ν(U Ξ)为例, 其离散形式为55Ν(U Ξ) = [ (U Ξ) i+ 2, j - 2 (U Ξ) i+ 1, j + 9 (U Ξ) i, j - 10 (U Ξ) i- 1, j + 2 (U Ξ) i- 2, j ]ö6△Ν, 　　U i, j > 0[ - 2 (U Ξ) i+ 2, j + 10 (U Ξ) i+ 1, j - 9 (U Ξ) i, j + 2 (U Ξ) i- 1, j - (U Ξ) i- 2, j ]ö6△Ν, U i, j < 0方程 (3)中的粘性项和方程 (4)采用二阶中心差分格式. 方程 (3)的时间推进为全隐式的.这样得到的差分方程, 其系数矩阵分别为 9 对角和 5 对角的, 原方法中的求解器对 9 对角矩阵已无效. 这里, 我们对所有的差分方程统一采用对宽带 (多于三对角) 的大型稀疏矩阵都适用的修正的不完全LU 分解的预处理共轭梯度法 (M ILU - CG)进行求解.预处理共轭梯度法 (PCG)是近年来计算数学领域内发展十分迅速的一类求解大型稀疏矩阵A 的极其高效的算法[ 11 ] , 其计算效率依赖于预处理器的选择. 目前, 比较先进的预处理器有不完全的LU 分解 ( ILU ). 定义下标集S A = { ( i, j ) : a ij: ≠ 0}. 若仅对矩阵A 中下标属于S A 的 a ij进行LU 分解, 这就是 ILU 分解. 由于A 的稀疏性, S A 中的元素很少, 因此 ILU 分解比LU 分解的计算量要小的多. 若将 ILU 分解中略去的量- a rira rrj , ( i, j )∈S A 加到主对角线元上, 就形成最简01 　 力 学 学 报 1997　年　第　29　卷© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd.   All rights reserved.单的一种更有效的所谓修正的不完全LU 分解方法 (M ILU ) [ 11, 12 ]. M ILU 分解的算法, 可简述如下A 0: = Afo r r: = 1 to n dobegin　　 fo r j Ε r and (r, j )∈ S A do　　　a rrj: = a r- 1rj　　fo r i> r and ( i, r)∈ S A do　　　a rir: = a r- 1ir öa rrr　　 fo r i, j > rand ( i, r)∈ S A and (r, j )∈ S A do　　　 begin　　　　q: = - a rira rrj　　　　if ( i, j )∈ S A then a rij: = a r- 1ij + q　　　　else a rii: = a rii+ q　　　endend　　在同样的精度控制下, M ILU - CG 方法的求解速度是传统的线松弛L SOR 结合LU 分解方法的 6～ 8 倍[ 11, 12 ] , 且这种计算速度上的优势, 随着矩阵规模和复杂性的增加, 会变得更加明显. 目前, 在大型矩阵方程的计算中, 后者已被前者所取代.涡法部分的计算仍与原方法一样, 详细过程请参阅文献[8 ]～ [10 ].3　计算结果分析我们计算了 R e= 1 000 和 a∈ (0. 5, 3. 25) 范围内, 均匀来流绕旋转圆柱的流动. 内区细网格, 其网格数为 144×240, 对应物理平面区域为 aΦ rΦ 3a , 0Φ ΗΦ 2Π. 外区粗网格, 内外区粗网格总数为 300×240, 对应物理平面区域为 aΦ rΦ 110a , 0Φ ΗΕΦ 2Π. 时间步长为△t= 0. 01.计算的时间范围为 t= 0→80.为验证本文方法的准确性, 我们将部分计算所得的流线分布图与已有的实验流场显示照片作了比较. 图 2 是 (a) a= 0. 5, t= 3, (b) = 3. 0, t= 4 时两者的对比 (实验资料均取自文献[ 3 ]).可见, 计算结果与实验结果是完全吻合的, 与其他计算结果[ 3～ 5, 7 ]相比, 也相当接近.随着旋转速度比 a 变化, 柱后尾流将呈现不同的流动结构.图 3, 4 分别是 a= 0. 5 和 1. 0 时, 在各个时刻的流线与涡量分布图. 与绕固定圆柱流动时的情形相仿, 柱后可见明显的从上、下柱面周期交替释放涡形成的类似 Karm an 涡街的旋涡结构,但其中心线随 a 的增加, 逐渐沿圆柱旋转方向偏移. 图 5 是 a = 2. 0 时, 在各个时刻的流线与涡量分布图. 此时, 下方释放的涡与上方释放的涡相比, 强度已大大减弱, 自 t = 30. 0 后, 流线的波动也已相当小了. 图 6 是 a = 3. 0 时, 在各个时刻的流线与涡量分布图. 此时, 下方涡的发展完全被抑制, 上方自第一个主涡释放后, 流动即趋于稳定, 流线呈平直状, 尾流中除了一个与上柱面相连的单边附着涡外, 再无其他涡释放. 我们的数值结果是支持文献[1, 3 ]的实验结论, 并与文[7 ]的数值结果相一致的. 文[5 ]的计算时间较短 (Φ 24) , 它的结论还有待商榷.11第　1　期 凌国平等: 绕旋转圆柱流动涡尾流结构和临界状态特性 　© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd.   All rights reserved.图 3　R e= 1 000, a= 0. 5, 在各时刻的流线与涡量分布图F ig. 3　Stream line and vo rtocity contours at differen t t im e instan ts fo rR e= 1 000, a= 0. 05　　可以推测在周期性交替双边释放涡结构, 和稳定的单边附着涡结构之间, 应该存在一种过渡性的所谓临界状态. 对 R e= 1 000, 我们在临界值 ac≈ 2. 2 处找到了这种状态. 图 7 就是此时在各个时刻的涡量分布图. 下方释放的涡已弱到可不显示出来的地步, 上方自第一个主涡释放后, 附着涡呈明显的藕节状, 而不是稳定状态时单边附着涡的那种几乎平直的形状. 原因是此时下方附着涡仍存在, 它对上方涡的结构还起着周期性的影响. 所以, 上方涡就成为象藕一样一节节的了, 随着时间的发展它会一节节地脱落. 当然, 这些脱落涡的强度比起主涡来要弱得多了.这很象文[ 5 ]说的单边释放涡的情形, 只是这种结构应该发生在临界状态, 而不是 a 很大的稳定状态.21 　 力 学 学 报 1997　年　第　29　卷© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd.   All rights reserved.图 6　R e= 1 000, a = 3. 0, 在各时刻的流线与涡量分布图F ig. 6　Stream line and vo rtocity contours at differen t t im e instan ts fo r R e= 1 000, a = 3. 031第　1　期 　 凌国平等: 绕旋转圆柱流动涡尾流结构和临界状态特性© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd.   All rights reserved.图 7　R e= 1 000, a = 2. 2, 在各时刻的涡量分布图F ig. 7　V o rticity contours at differen t t im einstan ts fo r R e= 1 000, a= 2. 2　　图 8, 图 9 分别是对不同的 a (= 0. 5, 1. 0, 2. 0, 3. 0) 阻力系数C d 和升力系数C Ι (× (- 1) )随时间变化的曲线. 对于 a= 0. 5, 1. 0 小于临界值的情形, C d , C Ι 呈现较大幅度的同频周期波动(注意: 与绕固定圆柱流动时C d 变化频率是 C Ι 的 2 倍的情形不同). 对于 a = 3. 0＄大于临界值的情形, C d , C Ι 趋于稳定值. 对于 a = 2. 0 接近临界值的情形, C d , C Ι 有微幅的同频周期波动.　　　　图 8　阻力系数随时间的变化　　　　　　　　　　　　 图 9　升力系数随时间的变化　　　R e= 1 000; a = 0. 5, 1. 0, 2. 0, 3. 0　　　　　　　　　R e= 1 000; a = 0. 5, 1. 0, 2. 0, 3. 0　　F ig. 8　T he variat ion of the drag fo rce coefficien t 　　　 F ig. 9　T he variat ion of the lift fo rce coefficien t　　w ith tim e fo r R e= 1 000; a = 0. 5, 1. 0, 2. 0, 3. 0　　　　　w ith tim e fo rR e= 1 000; a = 0. 5, 1. 0, 2. 0, 3. 0图 10　平均升阻比C ΙöC d 随 a 的变化F ig. 10　T he variat ion of C ΙöC d w ith a 　　研究绕旋转圆柱流动的一个重要目的, 是实现用旋转圆柱的方法进行流动控制优化的问题.从图 8, 图 9 可以看出, 随着 a 的增加, 升力、阻力均呈增加趋势, 但增加的幅度不一样. 因此,有可能在某一 a 处, 升阻比可获得极大值. 图 10是R e 分别为 1 000, 200 时, 平均升力系数C Ι (×(- 1) ) 和平均阻力系数 C d 之比 C ΙöC d 随 a 的变化曲线图. 可见, 平均升阻比 C ΙöC d 都将在临界值附近取得极大值. 这再一次表明了研究临界状态的流动特性, 对流动控制优化问题的重要性.41 　 力 学 学 报 　 1997　年　第　29　卷© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd.   All rights reserved.4　结　论1) 旋转速度比 a 对柱后尾流的旋涡结构和升力、阻力系数的变化, 有着决定性的影响.2) 存在临界的 ac. 当 a< ac 时, 尾流呈类似 Karm an 涡街的周期交替双边释放涡结构, 阻力、升力呈较大幅度同频周期波动. 当 a > ac 时, 尾流呈稳定的单边附着涡结构, 阻力、升力呈稳定值. 当 a≈ ac 时, 尾流呈节状单边释放涡结构, 阻力、升力呈微幅同频周期波动.3) 最大升阻比可在临界状态附近得到.4) 对高阶差分格式形成的宽带矩阵方程, 推荐采用高效的带预处理器的共轭梯度法 (如本文的M ILU - CG 方法)求解.参　考　文　献1　Coutanceau M , M enard C. Influence of ro tation on the near- w akedevelopm ent beh ind an impulsively started circu lar cylin2der. J . F lu idM ech. 1985, 158: 399～ 4462　Badr HM , D ennis SCR. T im e- dependent viscous flow past an impulsively started ro tating and translat ing circu lar cylinder.J . F lu id M ech. 1985, 158: 447～ 4883　Badr HM , Coutanceau M , D ennis SCR , et al. U nsteady flow past aro tating circu lar cylinder at Reyno lds num bers 103 and104. J . F lu id M ech. 1990, 220: 459～ 4844　Chang CC, Chern R - L. V o rtex shedding from an impulsively started ro tating and translat ing circu lar cylinder. J . F lu idM ech. 1991, 233: 265～ 2985　Chen Y- M , O u Y- R , Pearlstein A J. D evelopm ent of the w ake beh ind a circu lar cylinder impulsively started in to ro tatingand rectilinear mo tion. J . F lu id M ech. 1993, 253: 449～ 4846　O u Y- R. M athem atical modeling and num erical sim ulation inex ternal flow contro l. F low Contro l, Sp ringer- V erlag, 19957　Chew YT , Cheng M , L uo SC. A num erical study of flow past a ro tating circu lar cylinder using a hybrid vo rtex schem e. J.F lu id M ech. 1995, 299: 35～ 718　凌国灿, 凌国平, 王运平. 钝体分离旋涡流动的区域分解, 杂交数值模拟— I. 理论方法及其应用. 中国科学, A 辑, 1992(1) : 70～ 799　凌国灿, 凌国平, 顾启扬. 高雷诺数下绕钝体分离流动的一种新数值算法. 华中理工大学学报, 1992, 20 (4) : 87～ 9310　L ing Guo- p ing, L ing Guo- can. A num erical study of the vo rtex mo tion in o scillat ing flow around a circu lar cylinder atlow and m iddleKc num bers. A cta M echan ica S in ica, 1994, 10 (3) : 212～ 21911　吕涛, 石济民, 林振宝. 区域分解算法. 北京: 科学出版社, 199212　石济民, 林振宝, 岳兴业. M ILU - CG 方法在三维三相黑油模型中的应用. 石油勘探与开发, 1994, 21 (6) : 52～ 5851第　1　期 　 凌国平等: 绕旋转圆柱流动涡尾流结构和临界状态特性 　© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd.   All rights reserved.NUM ER ICAL STUDY ON VORTEX PATTERNS AND CR ITICALSTATE CHARACTERS IN THE W AKE OF A ROTATINGC IRCULAR CYL IND ERL ing Guop ing(D ep artm en t of M athem atics, S uz hou U niversity , S uz hou 215006, Ch ina)L ing Guocan( Institu te of M echan ics, Ch inese A cad em y of S ciences,B eij ing 100080, Ch ina)Abstract　A hyb rid fin ite d ifference and vo rtex m ethod, w h ich is based on dom ain decom po si2t ion and p ropo sed by the au tho rs (1992) , is u sed fo rca lcu la t ing the flow around a ro ta t ing cir2cu lar cylinder a t R eyno ldsnum ber R e= 1 000 and the ro ta t ion velocity ra t io a∈ (0. 5, 3. 25) re2spect ively. Som e im p rovm en ts are m ade in so lu t ion p rocess: ah igh o rder fu ll im p licit d ifferencea lgo rithm is adop ted, and the So lvers is rep laced by a very efficien t m odified incom p lete LU de2com po sit ion con juga te grad ien t m ethod (M ILU - CG). T he long- t im e fu llydeveloped fea tu resabou t the varia t ion s of the vo rtex pat tern s in thew ake , and drag, lif t fo rces on the cylinder areg iven. T he ca lcu la ted stream line con tou r p lo ts agreed w ellw ith the experim en ta lvisua lized flowp ictu res. T he ex istence of the crit ica l sta te isnum erica lly show n fo r the first t im e, and the vo r2tex pa t tern s a t co rresponding sta tes are g iven, too . T he m ax im um lif t to drag fo rcera t io can beob ta ined near the crit ica l sta te.Key words　 ro ta t ing circu la r cylinder, vo rtex w ake flow , fin ite d ifference m ethod, vo rtexm ethod, p recondit ioned con juga te grad ien t m ethod, incom p lete LU decom po sit ion61 　 力 学 学 报 　 1997　年　第　29　卷© 1995-2004 Tsinghua Tongfang Optical Disc Co., Ltd.   All rights reserved.

绕旋转圆柱流动涡尾流结构和临界状态特性

http://dspace.imech.ac.cn/bitstream/311007/35410/1/Paper0199.pdf

绕旋转圆柱流动涡尾流结构和临界状态特性

Abstract

Similar works

Full text

Available Versions

Institute Of Mechanics,Chinese Academy of Sciences