（流体力学专业论文）圆柱涡致振动的高精度数值计算.pdf

浙江人学硕l ：学位论文 摘要 本文采用任意拉格朗日欧拉( a l e ) 方法和区域分块法( d d m ) 数值求解了均 匀来流下圆柱的涡致振动问题。a l e 方法用于求解移动边界问题，而d d m 方法 适合于复杂流场的求解。流体运动通过原始变量n a v i e r - s t o k e s 方程描述，压力 采用p o i s s o n 方程形式，而柱体运动简化为弹簧质量一阻尼系统。计算网格采 用h 0 型非交错网格系统，易于采用严格的入流和出流边界条件。通过在每一 时间步分别求解n s 方程、p o i s s o n 方程和柱体运动方程，数值模拟流体和柱体 之间的非线性耦合运动。n s 方程的对流项和扩散项分别采用三阶迎风紧致格式 和四阶中心紧致格式离散，柱体运动采用经典龙格一库塔方法求解。 通过大量数据的计算，成功地捕捉到了“锁定”、“拍”和相位开关等现象， 详细分析了小雷诺数条件下圆柱尾涡形态、升阻力系数以及圆柱位移等在不同的 频率比、质量比、阻尼系数等条件下的变化规律。对单自由度和两自由度模型计 算结果的比较分析证明，尽管横向振动起主导作用，单流向振动对涡致振动有一 定的影响。另外，本文也对三维圆球绕流问题进行了数值模拟，为三维涡致振动 的研究做了铺垫。 关键词：a l e ；涡致振动；分块耦合方法；圆柱绕流；紧致格式；非交错网格 圆球绕流 t 本文受国家自然科学基金1 0 2 7 2 0 9 4 资助 3 浙江人学硕士学位论文 a b s t r a c t an u m e r i c a ls t u d yh a sb e e nc o n d u c t e di n t ov o r t e x i n d u c e dv i b r a t i o n s o fa ne l a s t i cc i r c u l a rc y l i n d e r ，u s i n gt h ea r b i t r a r yl a g r a n g i a n e u l e r i a n f o r m u l a t i o n ( a l e ) a n d t h ed o m a i n d e c o m p o s i t i o nm e t h o d ( d c ，哪t h ea l e f o r m u l a t i o ni sag e n e r a l i z e dc o n c e p tt od e a lw i t ht h ep r o b le mo fm o v in g b o u n d a r ie sa r o u n dac o n t i n u u m ，a n dt h ed u mm e t h o disa d v a n t a g e o u st os o l v e s u c hc o m p l i c a t e df l o wp r o b l e m t h ef l u i dm o t i o ni sd e s c r i b e db yt h e i n c o m p r e s s i b l en a v i e r s t o k e s ( n s ) e q u a t i o n s ，a n dt h ec y l i n d e rm o t i o ni s m o d e l e d b y a s p r i n g d a m p e r m a s ss y s t e m t h ep r i m i t i r ev a r i a b l e a n d p r e s s u r e p o i s s o ne q u a t i o nf o r m u l a t i o ni s e m p o y e d f o rt h en u m e r i c a l s o l u t i o no ft h ens e q u a t i o n s ，i nt e n s o rf o r m s ，o n a nh 一0 t y p e o f n o n s t a g g e r e dg r i d st h a tm a k e si tc o n v e n i e n tt oa p p l yas t r i c ti n f l o wa n d o u t f l o wb o u n d a r yc o n d i t i o n t h en o n 一1 i n e a rf l u i d s t r u c t u r ei n t e r a c t i o n i ss i m u l a t e dn u m e r i c a l l yb ys o l v i n gt h en - se q u a t i o n s ，t h e p r e s s u r e e q u a t i o na n dt h ec y l i n d e rm o t i o n e q u a t i o ni nt u r nf o re v e r yt i m es t e p t h en - se q u a t i o n sc o n v e c t i o nt e r ma n dd i s s i p a t i o nt e r ma r ed i s c r e t i z e d u s i n gt h et h i r d o r d e ru p w i n dc o m p a c ts c h e m ea n dt h ef o u r t h o r d e rc e n t r a l c o m p a c ts c h e m e ，r e s p e c t i v e l y t h ec y l i n d e rm o t i o n e q u a t i o ni ss o l v e d u s i n gt h er u n g e k u t t am e t h o d b y n u m e r o u s c o m p u t a t i o n a ld a t a ， w es u c c e s s f u ll y c a p t u r e d t h e “l o c k e d i n ”，“b e a t ”a n d “p h a s es w i t c h ”p h e n o m e n a ，a n dt h ep r e s e n t d i s c u s s e si nd e t a i lt h ev o r t e xs t r u c t u r e ，t h eu n s t e a d y1 i f ta n dd r a go n t h ec y li n d e r ，a n dt h ec y l i n d e r sd i s p l a c e m e n ta tv a r i o u sc o n d i t i o n so f t h ef r e q u e n c yr a t i o ，t h em a s sr a t i oa n dt h ed a m pc o e f f i c l e n ti nal o w r e y n o l d s r a n g e t h er e s u l t sa b o u tt h eo n e d e g r e e o f f r e e d o mm o d e la n d t w o d e g r e e o f f r e e d o m m o d e la r e c o m p a r e d ：t h ec o m p a r i s o n s h o w st h e s t r e a m w i s em o t i o no ft h e b o d y d o e si n f l u e n c et h e v o r t e x i n d u c e d v a r i a t i o n s ，e v e nt h o u g ht h et r a n s v e r s em o t i o ni sd o m i n a n tf o rab l u f fb o d y i nac r o s s f l o w b e s i d e s ，a3 - dn u m e r i c a ls i m u l a t i o no fu n i f o r mf l o wp a s t 4 浙江人学硕l 学位论文 a s p h e r e i s p r e s e n t e d ，t h a tp r e p a r e s f o r t h e i n v e s t i g a t i o n o f3 - d v o r t e x i n d u e e dv i b r a t i o n s k e yw o r d ：a l e ：v o r t e x i n d u c e dv i b r a t i o n s ：d o m a i nd e e o m p o s i t i o nm e t h o d ： f l o wo v e rac i r c u l a rc y l i n d e r ，c o m p a c ts c h e m e s ：n o n s t a g g e r e dg r i d s ：f l e w o v e ras p h e r e 4 t h e p r e s e n t i ss u p p o s e db y n a t u r a ls c i e n c ef o u n d a t i o no f c h i n a1 0 2 7 2 0 9 4 5 浙江火学硕一l 学位论文 第一章绪论 1 1 研究背景和前人成果 流体力学问题中流场的复杂性与其控制方程的高度非线性导致在许多带有 一般性的问题上，直接求解解析解存在相当大的数学困难，甚至有时几乎不可能。 随着现代电子计算机的飞速发展以及计算方法的不断进步，对一些特定流场下的 模型方程的直接数值模拟得到了巨大发展，并逐渐形成了一门完整的学科一计算 流体力学。近年来，随着计算方法和相关技术的不断改进，数值模拟在流体力学 研究中发挥了越来越重要的作用，计算流体力学已经成为研究流体力学的理论、 实验和计算三大支柱领域之一。这其中，有限差分法( f d m ) 是其中具有代表 性，也是应用广泛的一个分支。 流固耦合问题具有很强烈的应用背景和重要的学术价值，特别是近年来海洋 工程和风工程等实际应用的需要，加之实验技术的提高，使得流固耦合问题的研 究成为了现代流体力学发展的一个挑战性问题。其中，一个典型的流固耦合问题 就是由于结构涡脱落导致的涡致振动。在绕流钝体的尾部通常有周期性的涡脱落 发生，对于固定钝体的涡脱落现象，通常使用无量纲化斯特劳哈尔( s t r o u h a l ) 数来描述：s t = 工d u ，其中z 为涡脱落频率，u 为无穷远处来流速度，d 为 钝体尺度。伴随着涡的脱落，将产生周期性的作用力作用于该钝体，而当钝体为 弹性支撑，且其自振频率接近自然涡脱落频率的时候，将产生钝体结构的自激振 动，导致钝体结构的振动幅度剧烈增大，甚至造成结构的破坏。另一方面，钝体 和周围流体的耦合作用在某种条件下也可以抑制结构的振动，这在实际工程领域 可用于设备振动的控制。 鉴于涡致振动在工程中的强烈应用背景，前人已经在实验和数值方面做了一 定的工作。早期的实验研究主要集中在刚性结构的绕流( r i c h t e r & n a u d a s c h e r 1 9 7 6 【1 1 1 ；s o & s a v k a r1 9 8 1 t 1 羽；b a b a ne ta 1 1 9 8 9 1 1 3 1 ) ，后来则丌始了对弹性支撑结 构的绕流实验研究( a n a g n o s t o p o u l o s & b e a r m a n1 9 9 2 1 1 4 1 ，f e i r e i s e ne ta 1 1 9 9 4 1 1 5 1 ： w e s t & a p e l t1 9 9 7 1 6 1 ) 。对涡致振动的实验研究大致可以分为两个方向。一部分 实验研究主要关心结构振动为尾涡行为的影响，而使结构以指定的振幅和频率强 6 浙江大学硕士学位论立 迫振动，具有代表性的工作有k o o p m a n ( 1 9 6 7 ) 17 1 ，g r i f f i n & v o t a w ( 19 7 2 ) b 8 1 ，g r i f f i n & r a m b e r g ( 1 9 7 4 ) 1 9 】等。研究发现，结构的振动将改变尾涡的模式和形态，以及 涡与涡之间的间距。由于这些实验研究都没有考虑流体对结构的反馈作用，因此 所获得的数据和结果，并不是都能够准确刻画流固耦合的涡致振动现象。同时， 也有相当部分的实验研究直接对流固耦合的涡致振动进行了模拟，以确保实验结 构能够真实地刻画流体和结构之间复杂的非线性耦合问题，主要代表性的研究是 对弹性支撑圆柱的涡致振动研究，包括有f e n g ( 1 9 6 8 ) 睇，g r i f f i n & k o o p m a n n ( 1 9 7 7 ) 1 2 “，g r i f f i n ( 1 9 8 0 ) ”1 ，g r i f f i n & r a m b e r g ( 1 9 8 2 ) 2 3 】等。 另一方面，对涡致振动的数值研究则相对要少的多。与实验研究类似，一部 分研究人员主要侧重于研究振动结构对流场的影响。陆夕云等( 1 9 9 3 ) 1 2 4 】通过求解 涡量流函数形式的n s 方程，对圆柱横向强迫振动下的圆柱近迹尾涡结构进行 了数值研究，成功捕获了“锁定”、“相位开关”等现象，并对不同振动频率、振 幅下的复杂尾涡结构进行了分类，但由于涡量流函数法本身的限制，很难推广到 三维问题。刘松& 符松( 2 0 0 0 ) 2 5 通过动网格和静网格结合求解n s 方程数值研究 圆柱的流向强迫振荡，紧贴圆柱的内层网格跟随圆柱刚体运动，远离圆柱的外层 网格点则固定于惯性坐标系，中间一层过渡层可如弹性体做收张变形，这样，圆 柱的振动范围被限制在中间过渡层上，也使得其应用范围受到了很大限制。显然， 上述的研究工作并没有把流体和结构之间的耦合作用考虑在内，而实际上，在很 多的工程实际中，这种耦合作用经常起着非常重要的作用而不可忽略。具有代表 性的对流固耦合问题的研究有j a d i c e ta 1 1 9 9 8 2 剐，他们通过迭代求解流体控制方 程和结构运动方程，模拟流体和结构之间的耦合作用，成功求解了弹性支撑机翼 的绕流问题，揭示了流固耦合的非线性来源于机翼表面的移动边界条件； n e w m a n & k a r n i a d a k i s ( 1 9 9 7 ) 口9 】进一步改进了该方法，引进贴体坐标，对弹性电 缆的涡致振动现象进行了研究；c y z h o ue ta 1 1 9 9 8 【3 0 】则采用离散涡法计算了圆 柱的涡致振动，将流域用离散分布的涡来近似，流体和结构的耦合采用对每一时 间步的迭代进行模拟，在固定r e = 2 0 0 的条件下，揭示了决定圆柱涡致振动现象 的各种因素，发现结构的振动幅度与升力和圆柱位移的相位差有很大的关系，进 一步证明，振动幅度依赖于一个等效s g = 8 x 2 s t 2 a m + 阻尼参数，s t = 五d 乞为 s t r o u l a l 数，m = 聊p d 2 为无量纲化质量比，口为结构振动方程中的阻尼系数， 浙江人学硕j 学位论文 d 为圆柱直径，u 为均匀来流速度，f 为涡脱落频率，m 为无量纲化质量，p 为 流体密度。在国内，曹丰产& 项海帆【2 6 】吲采用动坐标系下的n s 方程，数值模 拟了弹性支撑圆柱的涡致振动现象，研究了锁定和非锁定区域圆柱位移、升力、 尾涡模态的变化，该方法将参照系固定在和圆柱一起做刚性运动的计算区域上， 计算网格随圆柱运动从而回避了移动边界导致的网格畸变问题。 在流固耦合问题的数值研究中，最难处理的是流体和刚性体之间的移动边界 问题，同时，这也是模拟非线性流固耦合的关键。任意拉格朗日一欧拉方法( a l e ) 是目前比较成功的用于解决移动边界和自由面问题的方法，a l e 方法在移动边 界上采用拉格朗r 描述很好地解决了流体和刚性体之间的一致性和平衡条件，同 时在远离移动边界的流体区域采用欧拉描述克服了拉格朗日描述不易解决的出 流外边界问题，也使流场的涡脱落等现象能够得到很好的描述。a l e 方法最早 由c w h i r t ( 1 9 7 4 ) 3 1 】推导，后来t j r h u g h e se t a l l 9 8 1 3 2 1 将其应用于解决自由面 问题，tb e l y t s c h k oe ta l1 9 8 2 【3 弱j d o n e ae ta l1 9 8 2 t 3 6 】，tn o m u r ae ta l1 9 9 2 【3 钔r ， w e ie ta l1 9 9 5 ( 3 3 1 等将其应用于解决流固耦合问题，但都主要应用在有限单元法。 1 2 本文研究内容与方法 本文的任务是对均匀来流下弹性支撑圆柱的涡致振动现象进行数值模拟，剖 析各种复杂的涡致振动现象和影响涡致振动的各种因素，以期达到对实际工程应 用的指导作用。 本文采用基于一般曲线坐标系的高精度紧致差分格式和分块耦合方法求解 任意拉格朗日一欧拉方式描述的原始变量二维不可压粘性流体的n s 方程，压 力采用p o i s s o n 方程形式求解；计算网格为h o 型非交错网格系统，该网格的优 点易于实施严格的进出口边界条件，且应用于a l e 方法时有利于节约计算机内 存和机时资源；柱体运动采用弹簧一阻尼一质量系统模拟，并通过经典龙格一库 塔方法数值求解求解。本文通过对流体和柱体之间相互耦合的数值模拟，观察柱 体的受力、运动特性和柱体尾部的涡脱落特性等，揭示各种复杂的涡脱落现象， 并详细讨论产生和影响这些现象的各种参数。 本文的第二章将详细介绍求解圆柱涡致振动问题的控制方程和计算方法，包 括a l e 方法、高精度紧致差分格式、区域分解算法、网格分块和生成、边界条 8 浙江人学坝j 。学位论文 件和初始条件以及具体计算过程等。第三章讨论刚性支撑钝体的绕流问题，包括 圆柱绕流问题和圆球绕流问题，其计算结果用于验证本文所使用方法和计算程序 的正确性，并为后面章节的涡致振动数值模拟提供比较的基础，圆球绕流的计算 是二维绕流问题的延伸，同时也为日后三维涡致振动问题的求解打下基础。由于 在大部分的涡致振动问题中，横向振动要比流向振动更为重要，目前大部分文献 都将横向振动作为研究重点，因此本文第四章着重讨论单自由度的涡致振动问 题，即圆柱仅允许在横向作来回运动；通过对不同雷诺数、不同频率比、不同质 量比的几组数据的计算，观察和分析一些著名的涡致振动现象，并对各种参数的 影响作详细剖析。第五章采用两自由度模型进行数值模拟，即在单自由度模型的 基础上考虑柱体的流向振动，对两种模型的计算结果进行比较分析，以验证流向 振动对涡致振动现象的影响。 9 浙江大学硕l 学位论文 第二章控制方程和计算方法 2 1 任意拉格朗日一欧拉法( a l e l 时间相关的不可压缩粘性流体的控制方程包括n a v i e r - s t o k e s 方程和连续性 方程，写成无量纲形式如下： p 憎m 善 _ 鹏+ 堡0 x s c z ， 婺= o ( 2 1 2 ) 僦 方程( 2 1 1 ) 为动量方程，方程( 2 1 2 ) 为连续性方程，“为应力张量，可定义如下： 伊即睁瓦o u j ( 2 1 3 ) 方程( 2 1 1 ) 一( 2 1 3 ) 仅在流体为不可压缩流体，密度为常量的情况下有效。珥为速 度矢量在i 方向的分量，t 为时间，6 j 为单位体力矢量在i 方向的分量，p 为流体 密度，p 表示压力，民是克罗内克尔( k r o n e c k e r ) 符号，为流体的粘性系数。 图2 1 a l e 方法描述 方程( 2 1 1 ) 仅在固定欧拉坐标系下有效，但不能应用于移动边界问题。我们 假设存在一个移动的参考坐标系x ，它可以随移动边界按某种规律变形，如图( 2 ，1 ) 浙江大学硕十学位论文 所示，假设一个流体质点在时问血内，从点a 经过位移a x 到达b ，而参考坐 标x 从点a 经过位移4 y 到c 点，则： x 一6 y = ( u v ) a t( 2 1 4 ) 这早，u 表示流体质点的速度，v 表示参考坐标x 的移动速度。此时，从点c 的角度看b 点，可以得出函数f ( x ，r 1 的全微分定义如下： o f = f ( x ，) + ，( x ，f ) = f ( x + ( 6 x 一6 y ) ，+ a t )( 2 1 5 ) 对上次在时间上进行二阶精度的泰勒展丌 o f 。f ( x , t ) - i - a f ( x , t ) ：f ( x , t ) + 丛鲨一a y ) + d f ( x ，, t ) a t + o ( 6 t 2 ) ( 2 1 6 ) o xd f 当f 趋向于零，将式( 2 1 4 ) 代入式( 2 1 6 ) ，得 要。笪些l + 堂些( u ( 2 1 - 7 ) d f 以i ，办 、 、7 从式( 2 1 7 ) 中可以看到，在任意拉格朗日一欧拉( a l e ) 描述下，( u v ) 的存在改 变了n s 方程的对流项，而由于应力张量、体力矢量和连续性方程对没有时间 的微分，它们在a l e 描述中跟欧拉描述一致。同时，我们可以注意到，当参考 坐标移动速度v = 0 时，式( 2 1 7 ) 与欧拉描述下的形式一致，而当坐标移动速度与 流体质点移动速度相等时( v = u ) ，式( 2 1 7 ) 与拉格朗日描述下的形式一致，这也 就是我们把该方法称为任意拉格朗日一欧拉方法的原因。综上所述，任意拉格朗 日一欧拉描述下的基本流体控制方程如下 p 陪+ ( u j - v j ) 针p b , + 百i 8 t i j b m ， 婴：0( 2 1 9 ) 2 2 控制方程及其差分离散 2 2 1 流体动量方程及其差分离散 时间相关的不可压缩粘性流体控制方程包括n a v i e r - s t o k e s 方程和连续性方 浙江大学填 学位论文 稗，以任意拉格朗日欧拉描述，可表示为如下的无量纲彤式： 百o u + ( 【f v 坶u = 一v p + 去v 2 u ( 2 2 1 ) v u = 0f 2 2 2 ) 其中u 为速度矢量，v 为网格速度矢量； 在二维任意曲线坐标系下，n s 方程又可表示为 上式中，u , v 为笛卡儿坐标系r 的速度分量，“v 为任意曲线罐标系下逆变速度 在# ，r l 方向的分量，为网格速度在任意曲线坐标系下逆变速度分量，其 定义为 ：妻鹄” 2 擘量 口：棚 【r = 仇u + r l j v【= 仉”。+ 叩，v 。 r 7 2 2 2 压力p o i s s o n 方程及其差分离散 对动量方程2 2 1 的等式两边取散度得： v 售+ ( u - v 妒i j ) + v2 p = 去vm ) ( 2 2 5 ) 交换散度和时间偏微分符号的位置，( 2 25 ) 式化为： 甲2 p 书( 1 v 2 u - ( uv ) 删) _ 掣 ( 2 26 ) 将续性方程( 2 2 2 ) 式代入l 式，则得a l e 描述下的压力p o i s s o n 方程： 俨肚v 去1 9 2 u - - ( u - v ) v u ) 阻7 ) 在不可压缩流动的问题中，连续性约束条件的满足是必须解决的一个问题，如果 连续性条件满足的不好，整个计算过程会g q t 计算瀑差的扩散而失败，即使得到 稳定解也是不精确的。根据本文采用的h a r l o w 的连续性余量修正 去，假设在第 n 步的连续性余量为d = v u 。，而第n + l 步的连续性余量为零，则经过修正的 压力p o i s s o n 方程为 驴 弘 。一如。一k 圳 协 铴 = 一一 朋 腑 圪 一 一 缈 心 峨 u 吒 砌一i毪加一西 塑垩盔兰堡土兰皇笙兰一 v ：脚( 1 v 2 u - ( u v u 卜詈 ( 2 2 8 ) 同样，二维压力p o i s s o n 方程的离散式可从方程组( 2 2 3 ) 推导： 首先，定义雅克比( j a c o b i ) 矩阵为 j ( j a c o b i ) = 嬲 ( 2 2 9 ) 专鲁+ 专( u 训妒专c 矿吲= 一专( 虽既+ 删+ 去扣亿：加、 专害+ 专( u 也咐专( 矿v 。= 一；( 钮w 一去扣 将( 2 2 1o ) 点乘v 孝再对善微分，加上将( 2 2 1 0 ) 点乘v 叩再对7 7 微分，可得压力泊 ar - 7 六w 抛+ 号( 【，一地+ 等( 矿一心+ 白詈+ 等( u 一乩也+ 号一) 】+ 石a 玎i 。g jo 研u 斗予一u 埘地+ 予( 矿一h + 予詈+ 号( u u m ) v ；+ 予( 矿一n = 一睦【等( 虽砖+ 仉岛) + 等嗨既+ 巩岛) 】+ 南号( 颤p ；+ 仇岛) + 号( 与p ；+ 吼岛) + 面1 晓3 【- 1 候v 2 甜+ 与v 2 v ) 】+ i ( g v 2 u + r y v 2 v ) ) ( 2 2 1 1 ) 考虑到粘性项肯定消失，则压力方程应为 丧b ( 量岳p + 最仇岛+ 与与姥+ 每仉岛) 1 + 茜b ( 玑最+ 玑仇岛+ 巩专, p ；+ r l y q y 岛h + 杀号c 最詈+ 每詈，+ 南哆魄詈+ 十仉詈，+ 妻号肛u m ) u ；+ ( v - 跏小号w 埘挑删吲咽+ + 杀号【( u 一心+ ( 矿一心心 + 号【一砜璩+ ( y 一心 ：o 佗2 1 2 ) 榍椐仵煮曲线举标变换理论，可知 1 3 浙江大学硕1 等位论文 9 1 = 占x 十 y y g u = 圆。+ ，hy g = 矾仉+ 叩。叩。 g “= 仉戋+ r y 乞, 定义： q = 鲁 号【( u 一) + ( y 一) 叫+ 号 ( u 一乩) + 唯一圪) k c r 2 = 南 予【( u 一砜) + ( 矿一圪) 叫+ 予 ( u 一) k + 缈一圪) y = 琵0 【_ 1 【文百o u + 彭争+ 丽0 了1 ( 仉百o u + 仉o 甜r ) = 毒虽粤+ 孚卅南虽c 芋+ 芋肛西a 嘴o ( i f ) + 南号 ：o d a ， 将( 2 2 1 3 ) 和( 2 2 1 4 ) 式代入( 2 2 1 2 ) 式，得 ( 2 2 ，1 3 ) r 2 2 1 4 ) 琵0 1g 既+ g m p q ) + 旦o t i 三j ( 毋+ 9 2 2 岛) 】+ 喜h + 南h + = 。( 2 2 1 5 ) 再定义： q = g “j = v 善田善，= 专( 茧量+ 岛每 哎= g ”j = v 孝o t t l j = ( 丢仉+ 乞巩) 8 t = 9 2 1 j = v , 7 田4 j = l ( r m x + q 占 岛= 9 2 2 ，= v 唧叩，= ( 仉仉+ 巩巩) r 2 2 1 6 ) 则任意拉格朗日一欧拉描述下，采用h a r l o w 的连续性余量修正项，压力p o i s s o n 方程在任意曲线坐标系下的离散表达式为： 专心职+ 口z 岛+ 吼) + 南( 届砖+ 屈岛+ 吧) + 詈= 。 1 4 r 2 2 1 7 ) 浙江大学倾j 学位论文 2 2 3 柱体运动方程 通常，一个刚性体在二维平面的运动可以由三个位移分量柬刻画： d = 西，疋，钟。( 2 2 1 8 ) 4 ，民分别表示刚性体在两个方向的平动分量，臼为刚性体的转动分量。在研究 柱体的涡致振动现象时，通常我们将柱体结构简化为如图2 2 所示的的弹簧一阻 尼一柱体模型，柱体在三个自由度上分别由弹簧和阻尼支撑，不考虑三个自由度 之间的耦合，则柱体的运动方程为 图2 2 埘口- i - c v + 胎= x ( 2 2 1 9 ) v 和口分别为柱体的速度矢量和加速度矢量 v = 万，口= d ( 2 2 2 0 ) m , c , k 分别为质量矩阵、阻尼矩阵和刚度矩阵，且均为常系数对角线矩阵。方程 ( 2 2 1 5 ) 的右端项x 为作用在刚性体g 的质心的牵引力矢量： x = x i ，x 2 ，m ( 2 2 2 1 ) 五，x ：分别为刚性体在两个方向所受到的作用力，m 为受到的转矩。 由于时间和篇幅的限制，本文仅讨论圆柱单自由度( 横向平动) 和两自由度 ( 两个方向的平动) 的涡致振动现象，由于上述模型的运动在各个自由度上是相 互独立的，因此方程( 2 2 1 5 ) 对本文研究问题依然有效。方程( 2 2 1 5 ) 在两自由 度下，可表示为如下形式： 其中，孝，分别为弹簧系统阻尼比和自振频率；和巧分别为柱体受到的阻力 冯 2 q 伽伽 巧 = = 锄 前面 十 + 锄鼬鲰甄 佗 佗 浙江人学硕j 学位论文 和升力；m 为柱体质量；。* 、t 嘶、膏嘶和y 嘶、j = ，嘶、j j 。为圆柱在两个方向的 位移、速度和加速度，各个量都已无量纲化。 2 3 高精度紧致差分格式 众所周知，流动问题是一个多尺度的物理问题，对于一般的流动函数，可以 用富氏展开或富氏变换来描述。为了能够模拟这种复杂的多尺度结构，使离散方 法能正确逼近富氏级数中的高频分量，从而使离散函数能够最大程度逼近真实 解，必须采用精细的网格或者高精度格式。1 。近些年随着计算流体力学的长足发 展，在格式研究和算法研究方面都有了很大进步。国内的专家教授作出了很大贡 献，如在高分辨率格式方面，张涵信教授的n n d 方法。7 1 就作出了很大贡献；在高 精度计算方法方面，付德薰、马延文教授的方法”“，陈国谦教授方法“”以及邓 4 , 冈t j 教授方法“”“”等都获得了广泛承认。国外的高精度方法主要有l e l e 法“、 p a d e 法、d e n n i s 法1 等等。 高精度差分格式从发展流派讲基本上分为二种类型：一类是对一阶导数孚， c 宜 二阶导数箕作高精度离散，这方面的文章主要有l e l e 、付德薰等人的文章； 戚 另一类是对模型方程一“庐= r ( d e n n i s 方法) 或模型方程2 a # 。= + r ( 陈 国谦方法) 进行高精度离散。有人称之为总体高精度方法( g l o b a l ) 。在时间精度 方面大都采用变化的j a m e s o n 四阶龙格一库塔法“”( r k ) 。中科院院士钟万勰教 授提出时间精确积分法，可惜在流动计算中还应用不多。下面，我们就几种常用 的高精度紧致差分格式作简单介绍。 2 3 1 中心差分紧致格式( c 0 8 ) 我们以非线性、带粘性及无粘性的b u r g e r s 方程作为模型方程，即： 丝+ 笪：上塑 西氖r e 缸2 厂：昙“2 ( 有粘性) 1 6 浙江火学顺 学位论文 甜( 0 ，f ) 2 1“0 , 0 2 o7 刈 ( 2 3 1 ) 塑+ 笪：0，：三“z( 无粘性) o to x 。 2 u ( o ，f ) = 1 ( 2 3 2 ) ( 觌= f ( 斟越 。， 其中f ，墨的内点格式为： 扭+ 詈f + 担= 去( 4 “+ 1 ) 壶& ，+ 岩墨+ 击乳。= 吉( 一z “，+ 1 ) ( z 。a ) 显式及四阶隐式等四种边界格式。模仿c a r p e n t e r m l 的做法在左边界0 点的物理 阻】k 】一k 】k 】 i = l ，n 阻川2 嘉阪m i - 1 ，n ( 2 3 5 ) 警= 一去【盯1 m m 】+ 矿1 面1 盯1 马m ， i - 1 ，n ( 2 3 6 ) 【一。 一1 【一：】【厂f 】= 【一。r 1 【一2 】l w ? i = - 。r l b ：】b ? ：m l ? ：l 。m ，l l 。u 。m 。1 w j 。u ，j ： ：； 。村。l ， 2 3 - 7 l 州h 1 ”i历可“j卅 n “ ll “刈 1 7 浙江人学硕1 ：学位论文 警= 一去川引+ i 1 ：r e l 吲- b 一忸】【引 在线法方程的系数矩阵m 。中含有变量u ，可考虑非线性影响。 同样，线性化b u r g e r s 方程( a u 新+ c 尝= 三k e 昙) 也可以离散为：优础积。 等= 一御- 1 缈f 】+ 吉去盯1 剐叫 = 一三h2 帅f 】+ 击h 上r e 【b 1 】1 【占2 酬 z - 其中，m ；b 。】_ 1 【4 ：】：l ：m 。l 口j 。m 加l ni l “n l ”n j “i n j ( 2 3 8 ) 2 3 2 傅德薰的三阶上风紧致格式( u c d a 3 ) 其内点格式为：詈+ 三3f j 一= i 1 咭1 乃。+ j 2 乃一i 5 乃一，) 采用的边界点格式为三阶显式： _ = 去( 一詈 + 詈，2 一；如+ ；，4 晶= ；卜：钿+ ：抽一芸钿+ 罢厶， 2 3 3 傅德薰的五阶上风格式( u c d a 5 ) ；巧十2 5 f 6 鬲1 ( 奶- 2 - - 4 4 f j _ + 3 6 f j + l 砺。一乃+ ：) 边界用四阶显式格式： e = 面1 ( 一2 5 f , + 4 8 f 2 3 6 十1 6 a 3 a ) e = 一- 2 5 f o + 4 毪，：t 一，一s 6 f 一：+ 16 _ 一，一玑。) 边界相邻的点上用四阶隐式格式： + 3 疋= 去( _ 1 7 z + 9 a 十9 a 一，4 ) ( 2 3 1 0 ) 浙江大学硕i 学位论文 只+ 3 只一= 一去( 一1 7 兀+ 9 六一+ 9 兀一：一f - 3 ) 2 4 区域分解算法与网格分块 2 4 1 区域分解方法概述 区域分解算法( d o m a i nd e c o m p o s i t i o nm e t h o d ，d 删) 是八十年代崛起的新 方向，由于该方法能将大型问题分解为小型问题、复杂边值问题分解为简单边值 问题、串行问题分解为并行问题，因此1 9 8 5 年以后研究渐趋活跃。进入九十年 代，区域分解算法已成为计算数学的热门领域，其趋势方兴未艾。任安禄发展了 分块三维贴体坐标数值生成的方法，并由此进行了复杂流场的分块耦合求解 4 7 j 4 8 4 9 o 为了适应流场的区域分块，用非交叉网格系统，a b d a l l a h 等发表了一 系列的文章”“”。1 ，论述了非交叉网格下的方程形式和离散格式等，这为我们研 究复杂流场的分块耦合求解提供了有用的工具，采用f i b d a l l a h 所提出的非交叉 网格下的方程形式和离散格式，流场变量定义在同一网格节点上，使得解在各子 块间的传递十分容易，促使我们成功地发展了复杂流场的分块耦合求解方法。 区域分解算法特别受关注是因为它具有其它方法无以比拟的优越性： 1 它把大问题化为若干小问题，缩小计算规模。 2 子区域形状如果规则( 如长方形) ，其上或者允许使用熟知的快速算法， 如快速f o u r i e r 变换( f f t ) 、谱方法等；或者已经有解这类规则问题的 高效率软件备用。 3 允许使用局部拟一致网格，无需用整体拟一致网格。甚至各子域可以用 不同离散方法进行计算。这对于形态极不规则的问题，如锅炉燃烧问题： 炉体部分与烟筒部分几何尺寸相差很大，整体计算为了对付烟筒部分， 不得不把网格加得很密，而区域分解算法可以把这两部分分别处理，具 有很大得灵活性。 4 允许在不同子域选用不同的数学模型，以便整体模型更适合于工程物理 实际情况。如气体绕飞行体流动，在边界层附近为粘性流，在边界层外 为无粘流，二者有不同的数学模型，使用区域分解算法易于在不同子域 选用更合于实际的模型。 1 9 浙江人学硕l 学位论文 5 算法是高度并行的，即计算的主要步骤是在各子域内独立进行的。 6 对对称区域问题有更简单的区域分解算法。 在计算流体力学领域，使用分块求解方法有利于边界条件的实现和贴体坐 标系统的生成。而这两点在整个数值模拟中起着至关重要的作用。大量的计算实 践告诉我们，使用不适当的网格系统会破坏质量、动量或动能的守恒，甚至导致 解的中断。网格生成的成功有时可理解为整个数值模拟的成功。然而，在大多数 的情况下，因为无法形成总体的贴体坐标系统，人们不得不回避这些复杂的三维 流动问题，或者作一些近似处理来进行研究。另一方面可以使得计算使用内存大 大节省，使得在小机器上能够解决大型的三维复杂流场的计算。 简而言之，d d m 的基本思路是：将求解区域q 分解成为若干个子区域： 五：u 豆，每个子域面，都具有相对简单规则的几何形状，于是原问题的求解转 t = l 化为在子域上求解，再通过一定的方法使子域间的信息得以传递，实现解的耦合， 最终得到全区域上的解。对于子域之间的耦合一般可分为重叠对接耦合与非重叠 对接耦合。如图1 ( a ) 为重叠对接耦合，块i 和块i i 间重叠一层网格，块i 的拟边 界为块i i 的内点，在完成每块的单独计算之后，只需将块i 的拟边界上的值用块 i i 内点的值替代即可，同样，块i i 的拟边界为块1 的内点，也只需用同样方法处 理即完成这两块之间的耦合。图1 ( b ) 为非重叠对接耦合，块i 与块i i 间仅在交界 处对接，无重叠网格，这时我们采用d i r i c h l e t - n e u m a n n 耦合，在块的拟边界上 采用导数相等( n e u m a n n 问题) 或插值( d i r i c h l e t 问题) 来得到。 l l i l 1 l 8 重叠耦合 图2 3 区域耦合方式 有关i ) d m 的理论分析、收敛性证明等已有不少论述， 2 4 2 圆柱涡致振动的区域分块 b 非重叠耩台 可参见文“。 为解决圆柱涡致振动中的移动边界问题，根据任意拉格朗日一欧拉方法的思 渐江大学硕l 学位论文 路，a l e 区域的网格在每个时间步都将随圆柱体的运动更新，因此，若将繁个 计算区域都作为a l e 区域，势必造成比较低的计算效率。实际上，相对于整个 计算区域来讲，圆柱振动的区域要相对小的多，因此，我们可以仅对圆柱周围的 一个小区域应用a l e 方法，而对其他区域仍然采用欧拉描述的方法，这样，在 整个计算过程中，仅圆柱周围的小区域需要更新网格，而其他区域网格固定不动， 使得计算效率得到一定的提高。如图2 4 所示，我们的计算区域取圆柱上游6 倍 直径，下游2 5 倍直径，圆柱上下方各取6 倍直径；整个计算区域划分为9 个子 区域( 块) ，仅在包围圆柱的第9 块应用a l e 方法，而其他8 块采用欧拉描述。 块1 与块4 、块4 与块6 、块2 与块1 、块2 与块3 、块5 与块3 、块5 与块8 、 块7 与块6 、块7 与块8 之间采用重叠一层的耦合方式，方便块间的数据传递； 而块9 与块2 、4 、7 、6 之间由于很难重叠网格，可采用非重叠d n 耦合方式， 也可在计算过程中，将第9 块网格向周围延伸一层网格，起到类似于重叠网格的 作用，实践证明，后者的计算结果更为理想。 一 j16 j ：( u t 、_ 一 358 _ i 2 4 3 网格生成方法 图2 4 计算区域分块 计算网格的生成( g r i dg e n e r a t i o n ) 是数值模拟的前处理过程，对整个计 算过程起着十分重要的作用。一个与边界形状相适应的网格，或更进一步与解的 分辨相适应的网格能大大简化微分方程的计算，起到计算快、分辨率高的效果。 而一种不合适的网格系统则会造成收敛慢、分辨率低，甚至导致计算不稳定。因 此如何来划分计算网格，或者说如何生成网格是偏微分方程组数值求解过程中的 一个重要问题。对于流体力学方程组，各种现象十分复杂，如边界层、激波、压 缩性等等，网格的好坏对计算问题更起到十分重要的作用。 尽管网格生成技术( n u m e r i c a lg r i dg e n e r a t i o n ) 日益成熟，但当流动区 2 l 浙江人学硕l 学位论义 域的形状较复杂时，网格的生成会遇到困难。例如当流动区域是一个多连通域时， 就很难生成统一的贴体网格。同时，在几何形状变化较大的区域附近，如壁而的 凸、凹角处，生成的网格往往具有较大的曲率，正交性不好，贴体性也较差，这 些将严重影响到数值解的精度以及算法的稳定性。而另一方面，这些形状变化较 大的区域却往往是流动更复杂的区域，也是人们更感兴趣的部分。区域分解法能 够( d d m ) 很好地解决这些问题。同过合适的区域分解可以在每个子域上独立地 生成性态良好的贴体网格。 网格生成方法主要有复变量方法、代数方法和微分方程方法。微分方程方法 中，t h o m p s o n 法。”比较成熟并且适用范围广泛，而且已被成功应用于多块耦合 网格生成中。因此本文中的网格生成采用t h o