（机械电子工程专业论文）拉管结构绕流的数值模拟及理论分析.pdf

大连理工大学硕士学位论文 摘要 圆截面拉管构件广泛应用于岸边集装箱起重机、履带起重机等工程机械产品上，这 种结构件往往长细比较大，因此极易发生风致振动，强烈的振动能很快引起杆件两端固 定处发生疲劳破坏。在风雨共同作用下，拉管表面会形成周期振荡的水线，进而激发拉 管的大幅振动，这是一种比较特殊的风致振动现象，它引起的拉管振幅远大于其它类型 的风致振动，学术上称之为风雨激振。本文应用数值模拟和理论分析的方法对拉管风雨 激振现象进行了较详细的研究。研究内容如下： ( 1 ) 应用数值模拟方法研究水线位置角p 对拉管绕流的影响。首先建立了数值模拟 的计算模型并设定边界条件，然后选取雷诺应力( r s m ) 模型，应用a n s y sc f x1 1 0 计 算流体力学软件对水线位于拉管不同位置的情况做数值模拟，得出拉管表面压力分布、 升力系数、阻力系数随水线位置角0 的变化规律，最后与前人的试验数据进行对比，验 证数值模拟的正确性； ( 2 ) 选取雷诺应力( r s m ) 模型，对四种不同水线尺寸下的拉管绕流进行数值模拟，得 出绕流特性，并分析水线尺寸对拉管表面气动力系数的影响； ( 3 ) 分别选用标准k 一占模型和雷诺应力( r s m ) 模型对拉管绕流现象做数值模拟，根 据分析结果证明雷诺应力( r s m ) 模型更适合拉管风雨激振的研究，为今后此类研究中湍 流模型的适当选取提供依据； ( 4 ) 建立拉管风雨激振的单自由度模型，并使用m a t l a b 软件编制基于四阶龙格库塔 法的拉管运动微分方程求解程序。通过此程序计算四组算例，分别用来讨论水线初始位 置角包、拉管自振频率厂、入口风速砜、拉管空间姿态( 倾角a ，风向角夕) 对拉管 振动的影响，并与已有的试验结果进行对比，验证理论分析的正确性。 关键词：拉管；数值模拟；风雨激振；理论分析 拉管结构绕流的数值模拟及理论分析 n u m e r i c a ls t i m u l a t i o na n dt h e o r e t i c a l a n a l y s i so f f l o wa r o u n ds t a yt u b e a b s t r a c t t h es t a yt u b ew i t hc i r c u l a rs e c t i o ni sw i d e l yu s e di nc o n t a i n e rc r a n e ，c r a w l e rc r a n ea n d o t h e rc o n s t r u c t i o nm a c h i n e r i e s t h i sk i n do fc o m p o n e n tu s u a l l yh a sab i gs l e n d e r n e s sr a t i o ， w h i c hc o u l dp r o b a b l ym a k ew i n d i n d u c e d v i b r a t i o no c c u r s t h es t r o n gv i b r a t i o n ，h o w e v e r ， c a nc a u s et h ef a t i g u ec r a c k i n go fr o d sb o t h - e n d u n d e rt h ea c t i o no fw i n da n dr a i n ，t h e r ew i l l f o r mr i v u l e t 、i t hc y c l i c a lm o t i o ni nt h es u r f a c eo fs t a yt u b e ，w h i c hc a nc a u s el a r g ea m p l i t u d e o fs t a yt u b e t h i si sas p e c i a lt y p eo fw i n d i n d u c e dv i b r a t i o nw h i c hi sc a l l e dw i n d r a i n i n d u c e dv i b r a t i o n i tc o u l dc a u s em u c hl a r g e ra m p l i t u d et h a no t h e rt y p e so fv i b r a t i o n i nt h i s p a p e r , i tm a i n l ys t u d i e ss t a yt u b e sw i n d - r a i n - i n d u c e dv i b r a t i o nb ym e a n so fn u m e r i c a l s i m u l a t i o na n dt h e o r e t i c a la n a l y s i s t h em a i n w o r ki n c l u d e sa sf o l l o w s ： ( 1 ) t h es t u d ym a i n l yr e s e a r c h e st h ei n f l u e n c eo fa n g l ep o s i t i o no fr i v u l e to nf l o w a r o u n dt h es t a yt u b eb ym e a n so fn u m e r i c a ls i m u l a t i o n f i r s t l y ，e s t a b l i s h i n gt h ec o m p u t a t i o n a l m o d e la n ds e t t i n gt h eb o u n d a r yc o n d i t i o n t h e n ，a p p l y i n gt h er e y n o l d ss t r e s sm o d e l ( r s m ) t om a k et h en u m e r i c a ls i m u l a t i o no nd i f f e r e n tr i v u l e tp o s i t i o n sb ym e a n so fa n s y sc f x 11 0c o m p u t a t i o n a lf l u i dd y n a m i c ss o f t w a r e ，a n dg e t t i n gt h ec h a n g el a wo fd i s t r i b u t i o no f s u r f a c ep r e s s u r e ，l i f tc o e f f i c i e n ta n dd r a gc o e f f i c i e n ta l o n gw i t ht h ea n g l ep o s i t i o no fr i v u l e t f i n a l l y ，c o m p a r i n gt h ed a t aw i t ht h ef o r m e r st e s td a t at ov e r i f yt h ev a l i d i t yo ft h en u m e r i c a l s i m u l a t i o n ； ( 2 ) c h o o s i n gt h er e y n o l d ss t r e s sm o d e l ( r s m ) t om a k en u m e r i c a ls i m u l a t i o no ff l o w a r o u n ds t a yt u b e 、析t hf o u rd i f f e r e n ts i z eo fr i v u l e t t h e ng e t t i n gt h ef l o wc h a r a c t e r i s t i c sa n d a n a l y z i n gt h ei n f l u e n c eo fr i v u l e ts i z eo nt h es u r f a c ec o e f f i c i e n to fa e r o d y n a m i co fs t a yt u b e ； ( 3 )a n a l y z i n gt h es l i p s t r e a ma r o u n dt h es t a yt u b eb yn u m e r i c a ls i m u l a t i o nt 1 1 a ta d o p t s s t a n d a r dk - - em o d e la n dr e y n o l d ss t r e s sm o d e l t h e n ，a n a l y s i sc o n t r a s t i n gt h er e s u l t so ft w o m o d e l sd e m o n s t r a t e st h a tr e y n o l d ss t r e s sm o d e li sb e t t e r m e a n w h i l e ，t h ea n a l y z i n gr e s u l t p r o v i d e sr e f e r e n c et oa d o p tt u r b u l e n c em o d e l ； ( 4 ) s e t t i n gu pt h es o d fm o d e lo fs t a yt u b ea n du s i n gm a t l a bs o f t w a r et ow r i t e p r o g r a mf o rs o l v i n gt h ed i f f e r e n t i a le q u a t i o no fs t a yt u b em o v e m e n tb a s i n go nf o u r t ho r d e r r u n g e - k u t t am e t h o d c o m p u t i n gt h ep r o g r a mt oa n a l y z et h ee f f e c to ft h ei n i t i a la n g l e p o s i t i o no fr i v u l e t ，n a t u r a lf r e q u e n c yo fs t a yt u b e ，i n l e tw i n ds p e e du 0 ，s p a t i a la t t i t u d eo fs t a y t u b e ( i n c l i n a t i o n 口，w i n da n g l e 夕) ，w h i c hi n f l u e n c et h ev i b r a t i o no ft h es t a yt u b e t h e nt h e i i 大连理工大学硕士学位论文 c o m p a r i s o nb e t w e e nt h ea n a l y t i c a lr e s u l t sa n dt h ee x p e r i m e n t a ld a t aw a sc a r r i e do u tt o d e m o n s t r a t et h ev a l i d i t yo ft h e o r e t i c a la n a l y s i s k e yw o r d s ：s t a yt u b e ；n u m e r i c a ls i m u l a t i o n ；w i n d - r a i n i n d u c e dv i b r a t i o n ； t h e o r e t i c a la n a l y s i s 大连理工大学学位论文独创性声明 作者郑重声明：所呈交的学位论文，是本人在导师的指导下进行研究 工作所取得的成果。尽我所知，除文中已经注明引用内容和致谢的地方外， 本论文不包含其他个人或集体已经发表的研究成果，也不包含其他已申请 学位或其他用途使用过的成果。与我一同工作的同志对本研究所做的贡献 均已在论文中做了明确的说明并表示了谢意。 若有不实之处，本人愿意承担相关法律责任。 学位论文题目：撞篁结捡缝速鲍数焦搓拯丞堡诠佥堑 作者签名：盘壁。至日期：竺! 乞年j l 月生日 大连理工大学硕士学位论文 大连理工大学学位论文版权使用授权书 本人完全了解学校有关学位论文知识产权的规定，在校攻读学位期间 论文工作的知识产权属于大连理工大学，允许论文被查阅和借阅。学校有 权保留论文并向国家有关部门或机构送交论文的复印件和电子版，可以将 本学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、 缩印、或扫描等复制手段保存和汇编本学位论文。 学位论文题目：撞笪结捡缝速鲍熬焦搓拯区堡诠佥堑 一一7z一、 作者签名：益竺兰日期：三! 12 年l 月二生日 导师签名：兰乒递丛一 日期：三生年生月l 日 大连理丁大学硕士学位论文 1 绪论 11 引言 111 风灾害 风是空气从气压大的地方向气压小的地方流动形成的。气流一旦遇到结构的阻塞， 就形成高压气幕。风速愈大，对结构产生的压力也愈大，从而导致结构产生大的变形和 振动。如果结构物的抗风设计不当，或者产生过大的变形，就会使结构不能正常工作， 或者产生局部破坏，甚至整体破坏。随着经济的飞速发展，高、大、细、长的柔性工程 结构越柬越多，对于起重机结构亦是如此。在国内外，风载荷作用下的起重机事故每年 都时有发生。 2 0 0 3 年9 月1 2 号到1 3 号夜间袭击韩国的第1 4 号台风“鸣蝉”是韩国近一个世纪 内所遭遇慑强的台风，其风速达到每小时2 1 6 公里。给韩国约四分之一的地区造成了 严重的损害，死亡及失踪人数已经达到1 万人，直接经济损失将近1 3 亿美元。瞬间风 速达4 2 7 米，秒的狂风刮倒了釜山港区集装箱码头的1 0 多台重达近千吨的货物起重机。 1 9 9 1 年7 月1 9 日凌晨3 时3 0 分左右，青岛市区上空因强对流天气形成飓风，晟大风速 达3 07 米，秒，雨量1 8 0 毫米，码头断电，海水涌上岸壁，两台岸边集装箱起重机被吹 倒，两机相撞倒塌，直接经济损失1 5 0 0 余万元，并造成1 人死亡、1 人重伤、停产4 个 多月，成为青岛港最严重的一次灾祸。1 9 9 8 年7 月天津港集装箱公司因突发阵风使一台 进口的岸边集装箱起重机倾倒入海，直接损失数千万元。如图11 为一起风灾事故发生 后的现场图片。 图11 风灾事故 f i g1 1w i n dd i s a s t e r a c c i d e n t 拉管结构绕流的数值模拟及理论分析 对于大型的起重机，结构比较轻柔，而且钢材的阻尼小，对风的作用相对敏感，动 力特性更加突出。大型起重机又多安装在常年有大风的江边和沿海，更应考虑风载荷对 结构造成的不利影响。为了满足起重机设备的安全性、工作性能、及操作者的舒适性要 求，必须对起重机的风振问题有深入的研究，指导起重机的结构设计等工程实践活动。 我国现行的起重机设计规范制定于1 9 8 3 年，其中风载荷部分主要参考建筑结构的 风载荷规范，多年来一直未做修订。对起重机的设计来讲，单纯套用一般建筑结构的风 振研究成果显然是不合适的。所以对起重机进行风振研究不但具有重大的理论价值而且 能推动起重机械的技术进步，具有很大的经济效益和社会效益。 1 1 2 起重机拉管振动的主要类型 目前，圆截面拉管结构广泛应用与大型工程起重机中，例如大吨位履带起重机、大 吨位岸边集装箱起重机等。在环境因素及结构本身的固有特性达到一定条件时，细长结 构的拉管构件极易发生风致振动，拉管受风激作用有可能发生的振动形式有涡激振动、 抖振、驰振、尾流驰振和风雨激振等【l ，2 】。目前前四种振动的产生机理已经清楚了，但是 风雨激振的产生机理尚未得出定论，本文主要针对拉管风雨激振现象进行了研究。图1 2 为拉管风致振动类型的详细分类。 图1 2 拉管风致振动类型 f i g 1 2t y p e so fw i n d - i n d u c e dv i b r a t i o no fs t a yt u b e ( 1 ) 涡激振动( v o r t e xs h e d d i n g ) ：当流体绕过拉管时，在满足特定因素的条件下( 如 雷诺数、拉管几何尺寸等) ，流体会在拉管后方形成交替脱落的漩涡，从而引起拉管表 面受力状况的周期变化，当漩涡脱落频率与拉管的固有频率一致时，就会引起拉管的大 大连理工大学硕士学位论文 幅振动。漩涡脱落的频率( 又称斯超海尔频率) 由下式表示： c ， f ，= 兰兰 ( 1 1 ) d 式中，u 为风速；d 为拉管直径；s t 为无量纲参数，称为斯超海尔( s t r o u h a l ) 数， 对于圆截面构件，当雷诺数在1 0 x1 0 43 0 x1 0 5 之间时，s t 大致等于o 2 。当流体为湍 流时，s t 的大小还与流体的湍流强度有关，本文第三章做了详细说明。 当漩涡脱落频率正好等于拉索的某阶自振频率f i 【时，则会激发拉管的涡激共振。将 式1 1 做适当变形，可得出发生涡激共振的风速为： ，n u 一= 坐兰 ( 1 2 ) “ s 对于常见的拉管，设管的k 阶模态频率氏为2 h z ，索直径d 为o 1 6 m ，则可得k 阶 模态下发生涡激共振的风速u v 。等于1 5 m s 。它表明涡激共振能在很低的风速下发生， 并且对于高阶模态其发生涡激共振的风速也越高。 涡激共振引起的拉管振幅可由下式近似计算： 抛o o s q ( 射1u c p d 2 di 2 j 式中，m 为拉管每米质量；孝为拉管模态阻尼比；p 为空气密度； 系数。 ( 1 3 ) c l 是气动升力 由上式可知，增大拉管的质量和阻尼比可以降低拉管的振幅。参数m 善p d 2 称为 s c r u t o n 数，s c r u t o n 数越大涡激振动幅度越小。 ( 2 ) 抖振( b u f f e t i n gv i b r a t i o n ) ：当一根拉管处于另一根拉管的卡门涡列之中时可 发生抖振。两根并列拉管背风向拉管与顺风向拉管频率接近时可能发生抖振现象，抖振 实际上是一种顺风向共振，属于强迫振动。 ( 3 ) 驰振( g l l o p i n g ) - 驰振是发生在具有特殊截面形状的长细结构物的一种典型的 单自由度气动不稳定现象。驰振通常发生在截面形状特殊的结构上，在风的作用下具有 气动不稳定的性质，产生了负阻尼导致横风向运动失稳。驰振是一种自激振动。 有关驰振临界的风速问题，常运用d e nh a r t o g 判据分析： f 堕+ c n1 0 ( 1 4 ) l 砌“ 这是驰振不稳定性的必要条件。它的充分条件是： 2 m 蝣+ i 1p 凇f 孕+ c o1 雨是发生风雨激振的必要条件，一般发生在雨量为小到中雨的天气； 通常只在沿风向向下倾斜的拉索上发生风雨激振； 发生振动的拉索直径一般为1 2 0 m m 2 0 0 m m ，拉索表面材料性质，灰尘等对振动有 重要影响。 ( 2 ) 风雨激振振动特性 振动主要发生在拉索平面内，但也存在一个面外振动分量； 风雨激振发生时拉管的最大振幅一般为3 5 倍拉管直径； 风雨激振属于稳定的单频振动，振动频率主要发生在0 6 3 h z ，振动常以“拍 的形 式出现； 非常不稳定的多振型振动也经常发生，振动中不稳定响应有一个很长的周期 当风雨激振发生时，拉索的上下表面各有个一条雨水形成的水线，沿倾斜的拉索向 下流动。水线随着拉索的振动在拉索表面振荡。 ( 3 ) 风雨激振抑制方法 机械控制措施是一种通过在拉索上附加阻尼器来提高拉索结构阻尼的减振方法。目 前，附加阻尼器的方法广泛的应用于大型工程起重机的拉管结构中，如德马格的 c c 8 8 0 0 型号履带起重机的变幅拉管就在根部固定端布置了阻尼块。 结构控制措施是采用辅助索将不同拉索联结。研究表明柔性辅助索能有效提高拉索 的阻尼，并有一定的能量耗散；而刚性索几乎没有能量耗削9 j 。 空气动力学控制措施主要通过改变拉索的截面形状及对拉索表面进行特殊处理来改 善拉索在风雨作用下的空气动力学性能【1 1 1 。空气动力学控制措施具有费用低、几乎 大连理工大学硕士学位论文 不需维护等有点。但需注意的是某种气动措施对某类拉索风致振动有效，而对另外 类型的风致振动可能产生更不利的后果。杜晓庆【1 3 】在人工降雨试验中研究了在拉索 表面缠绕螺旋线对抑制拉索风雨激振的作用，发现螺旋线的直径、缠绕方向、和螺 旋线之间的间距对风雨激振的影响并不是单一变化的，要慎重选择，避免对风雨激 振形成放大作用。 1 5 本文工作 本文研究工作以计算流体力学知识为基础，应用i c e m 网格划分软件及a n s y s c f x 计算流体力学软件对带水线的拉管结构进行数值模拟。同时建立拉管及水线耦合状 态下的运动微分方程，并应用四阶龙格库塔法求解运动微分方程，讨论影响拉管风雨 激振的参数。 第一章主要介绍了风灾害对国民经济造成的重大损失，及研究风工程的必要性与重 要性。对比实际情况总结了起重机拉管结构振动的主要类型，并且重点总结了拉索风雨 激振的研究手段及现状，用于以下章节对拉管结构绕流的研究。 第二章简要介绍了计算流体动力学方法的理论基础及实际应用。描述了方程的离散 方法、计算流体力学控制方程。 第三章基于a n s y sc f x1 1 0 计算流体力学软件对圆截面拉管构件进行数值模拟分 析，详细研究了在风雨共同作用的情况下，带水线拉管在高雷诺数下的表面压力分布、 阻力系数和升力系数随水线位置、水线形状的变化规律，与前人现场实测及试验结果进 行了对比，验证了数值模拟结果的正确性。对比了标准k 一占模型和r s m 模型这两种湍 流模型下的计算结果，得出适合此类问题计算的湍流模型。 第四章建立了拉管风雨激振的单自由度模型，推导出拉管的运动微分方程，并且应 用四阶龙格库塔法编制m a t l a b 计算程序求解二阶微分方程。选取四组算例分别讨论了 水线初始位置角、拉管自振频率、入口风速、拉管空间姿态四种参数对拉管振动的影响， 并与已有的试验结果进行对比，验证了本文结果的正确性。 最后，给出了本文研究工作的主要结论和对将来研究工作的展望。 拉管结构绕流的数值模拟及理论分析 2 计算流体动力学( c f d ) 数值模拟方法 2 1计算流体动力学简介 计算流体力动力学( c o m p u t a t i o n a lf l u i dd y n a m i c s ) 是通过计算机数值计算和图像 显示，对包含有流体流动和热传导等相关物理现象的系统所做的分析，它所涉及的学科 有流体力学、偏微分方程的数学理论、计算几何、数值分析、计算机科学等，它的发展 促进了这些学科的进一步发展。它以电子计算机为工具，应用各种离散化的计算方法， 对流体力学的各类问题进行数值试验、计算机模拟和分析研究，以解决各种实际问题。 计算流体力学的基本特征是数值模拟和计算机试验，它从基本物理定理出发，在很大程 度上替代了耗资巨大的流体动力学实验设备。在科学研究和工程技术中应用广泛。 任何流体运动的规律都是由以下三个定律为基础的：质量守恒定律；动量守恒定律 和能量守恒定律。这些基本定律可由数学方程组来描述，如欧拉( e u l e r ) 方程，n s 方 程。采用数值计算方法，通过计算机求解这些数学方程，研究流体运动特性，给出流体 运动空间定常或非定常流动规律，是应用计算流体动力学方法解决问题的主要途径。 计算流体动力学的兴起促进了实验研究和理论分析方法的发展，为简化流动模型的 建立提供了更多的依据，使很多分析方法得到发展和完善，例如目前在飞机工业中应用 广泛的面元法就是一个很好的例子。然而，更重要的是计算流体动力学采用它独有的新 的研究方法，用数值模拟方法研究流体运动的基本物理特性。这种方法的特点如下【4 1 ： ( 1 ) 给出流体运动区域内的离散解，而不是解析解。这区别于一般理论分析方法； ( 2 ) 它的发展与计算机技术的发展直接相关。这是因为可能模拟的流体运动的复杂程度、 解决问题的广度和所能模拟的物理尺度以及给出解的精度，都与计算机速度、内存、运 算及输出图形的能力直接相关； ( 3 ) 若物理问题的提法( 包括数学方程及其相应的边界条件) 是正确的，则可在较广泛 的流动参数( 如马赫数、雷诺数、气体性质、模型尺度等) 范围内研究流体力学问题， 且能给出流场参数的定量结果。这常常是风洞实验和理论分析难以做到的，然而，要建 立正确的数学方程还必须与实验研究相结合。更重要的是实际问题中所求解的多维非线 性偏微分方程组十分复杂，其数值解的现有数学理论尚不够充分。严格的稳定性分析， 误差估计和收敛性理论的发展还跟不上数值模拟的发展。虽然关于广义解唯一性存在性 等问题的严格数学理论已取得了长足的进展，但还不足以对一些感兴趣的具体的复杂问 题给出明确的回答。所以在计算流体动力学中，仍必须依靠一些较简单的、线性化的、 与原问题有密切关系的模型方程的严格数学分析，以及依靠启发性的推理给出所求解问 大连理工大学硕士学位论文 题的数值解的理论依据。然后再依靠数值实验，地面实验和物理特性分析，验证计算方 法的可靠性，从而进一步改进计算方法。 事实上，实验研究、理论分析方法和数值模拟是研究流体运动规律的三种基本方法， 它们的发展是相互依赖，相互促进的。无论采用何种方法，研究流体运动的理论基础都 是计算流体动力力学，它的实现步骤如下【4 】： ( 1 ) 建立反映工程问题或物理问题本质的数学模型。即建立反映问题的各个量之间关系 的微分方程及相应的定解条件，这是数值模拟的出发点。流体的基本控制方程通常包括 质量守恒方程、动量守恒方程、能量守恒方程，以及这些方程相应的定解条件。 ( 2 ) 寻求高效、高精度的求解方法，即建立针对控制方程的数值离散化方法。如有限差 分法、有限元法、有限体积法等。这里的计算方法不仅包括微分方程的离散化方法及求 解方法，还包括贴体坐标的建立，边界条件的处理等。以上内容是c f d 的核心。 ( 3 ) 应用商业软件或编制程序进行计算。包括计算网格的划分、初始条件和边界条件的 输入、控制参数的设定等。 ( 4 ) 计算结果后处理。提取计算结果，必要时以图表等直观方式输出结果。 2 2 计算区域离散方法 随着计算机技术和计算方法的发展，许多复杂的工程问题都可以采用区域离散化的 数值计算并借助计算机得到满足工程要求的数值解。数值模拟技术是现代工程学形成和 发展的主要动力之一【m 1 5 】。 区域离散化就是用一组有限个离散的点来代替原来连续的空间。具体操作过程是把 所计算的区域划分成许多互不重叠的子区域，确定每个子区域的节点位置和该节点所代 表的控制体积。节点代表了需要求解的未知物理量的几何位置，控制体积是应用控制方 程或守恒定律的最小几何单位。一般把节点看成控制体积的代表。控制体积和子区域并 不总是重合的。在区域离散化的过程中，由一系列与坐标轴相应的直线或曲线簇所划分 出来的小区域称为子区域。网格是离散的基础，网格节点是离散化物理量的存储位置。 常用的离散化方法有有限差分法、有限元法和有限体积法。对这三种方法分别介绍 如下。 有限差分法( f d m ) 是数值解法中最经典的方法。它是将求解区域划分为差分网格， 用有限个网格节点代替连续的求解域，然后将偏微分方程( 控制方程) 的导数用差商代 替，推导出含有离散点上有限个未知数的差分方程组。这种方法的产生和发展比较早， 也比较成熟，较多用于求解双曲线和抛物线型问题。用它求解边界条件复杂，尤其是椭 圆形问题不如有限元法或有限体积法方便。构造差分的方法有多种形式，目前主要采用 拉管结构绕流的数值模拟及理论分析 的是泰勒级数展开法。其基本的差分表达式主要有四种形式：一阶向前差分、一阶向后 差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式 为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差 分计算格式。 有限元方法( f e m ) 的基础是变分原理和加权余量法，其基本求解思想是把计算域 划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插 值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成 的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函 数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来 随着计算机的发展逐渐用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分 为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线性 组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成 的，则整个计算域内的解可以看作是由所有单元上的近似解构成。根据所采用的权函数 和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、 矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边 形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数 等。不同的组合同样构成不同的有限元计算格式。对于有限元方法，其基本思路和解题 步骤可归纳为： ( 1 ) 建立积分方程：根据变分原理或方程余量与权函数正交化原理，建立与微分方程初 边值问题等价的积分表达式，这是有限元法的出发点。 ( 2 ) 区域单元剖分：根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相 互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作 量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点 的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。 ( 3 ) 确定单元基函数：根据单元中节点数目及对近似解精度的要求，选择满足一定插值 条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单 元具有规则的几何形状，在选取基函数时可遵循一定的法则。 ( 4 ) 单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再 将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数，即单元中各 节点的参数值的代数方程组，称为单元有限元方程。 ( 5 ) 总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则 进行累加，形成总体有限元方程。 大连理工大学硕士学位论文 ( 6 ) 边界条件的处理：一般边界条件有三种形式，分为本质边界条件( 狄里克雷边界条 件) 、自然边界条件( 黎曼边界条件) 、混合边界条件( 柯西边界条件) 。对于自然边界 条件，一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按 一定法则对总体有限元方程进行修正满足。 ( 7 ) 解有限元方程：根据边界条件修正的总体有限元方程组，是含所有待定未知量的封 闭方程组，采用适当的数值计算方法求解，可求得各节点的函数值。有限元法( f e m ) 是将一个连续的求解域任意分成适当形状的许多微小单元，并于各小单元分片构造插值 函数，然后根据极值原理( 变分或加权余量法) ，将问题的控制方程转化为所有单元上 的有限元方程，把总体的极值作为各单元极值之和，即将局部单元总体合成，形成嵌入 了指定边界条件的代数方程组，求解该方程组就得到各节点上待求的函数值。对椭圆形 问题有更好的适应性。有限元法求解的速度比有限差分法和有限体积法慢，在商用c f d 软件中应用并不广泛。 有限体积法( f v m ) ，又称为控制体积法，是将计算区域划分为网格，并使每个网 格点周围有一个互不重复的控制体积，将待解的微分方程对每个控制体积积分，从而得 到一组离散方程。其中的未知数是网格节点上的因变量。子域法加离散，就是有限体积 法的基本思路。有限体积法的基本思路易于理解，并能得出直接的物理解释。离散方程 的物理意义就是因变量在有限大小的控制体积中的守恒原理，如同微分方程表示因变量 在无限小的控制体积中的守恒原理一样。有限体积法得出的离散方程，要求因变量的积 分守恒对任意一组控制体积都得到满足，对整个计算区域自然也得到满足，这是有限体 积法吸引人的优点。有一些离散方法，如有限差分法，仅当网格极其细密时，离散方程 才满足积分守恒，而有限体积法即使在粗网格情况下，也显示出准确的积分守恒。 就离散方法而言，有限体积法可视作有限元法和有限差分法的中间产物。三者各有 所长。有限差分法：直观，理论成熟，精度可选，但是不规则区域处理繁琐，虽然网格 生成可以使f d m 应用于不规则区域，但是对区域的连续性等要求较严。使用f d m 的 好处在于易于编程，易于并行。有限元方法：适合处理复杂区域，精度可选。缺陷在于 内存和计算量巨大，并行不如f d m 和f v m 直观。不过f e m 的并行是当前和将来应用 的一个不错的方向。有限体积法：适于流体计算，可以应用于不规则网格，适于并行。 但是精度基本上只能是二阶了。f v m 的优势正逐渐显现出来，目前的主流c f d 商业软 件都是基于f v m 方法开发的。 拉管结构绕流的数值模拟及理论分析 2 3 流场基本控制方程 ( 1 ) 质量守恒方程 任何流动问题都必须满足质量守恒定律。该定律可表述为：单位时间内流体中质量 的增加，等于同一时间间隔内流入该微元体的净质量。按照这一定律，可以得出质量守 恒方程【1 9 1 ： 望- i - 塑尘+ 皇盟- t a ( 一p w ) ：0 ( 2 1 ) 一 一十一= 、厶上， 8 t瓠巩 。z 引入矢量符号挑( 口) ：誓+ 要+ 冬，式( 2 1 ) 写成： 蹴卯化 粤+ 挑( ) ：0 ( 2 2 ) a ” 有的文献使用符号v 表示散度，即v 口：d i v ( 口) ：誓+ 要+ 誓，这样，式( 2 2 ) o xo z 写成： 鲁棚伽) = 0 ( 2 3 ) 在式( 2 1 ) 至( 2 3 ) 中，p 是密度，是时间，, 是速度矢量，“、v 和w 是速度矢量甜 在x 、y 和z 方向的分量。 上面给出的是瞬态三维可压流体的质量守恒方程。若流体不可压，密度p 为常数， 式( 2 1 ) 变为： 丝+ 塑+ 塑：o ( 2 4 ) + + 一= u l 么q ， 瓠 硒 砚 若流动处于稳态，则密度p 不随时间变化，式( 2 1 ) 变为： a ( p u ) + 旦边一t a ( p w ) ：0 ( 2 5 ) a瓮8y o z 质量守恒方程( 2 1 ) 或( 2 2 ) 常称作连续方程。 大连理工大学硕士学位论文 ( 2 ) 动量守恒方程 动量守恒定律也是任何流动系统都必须满足的基本定律。该定律可表述为：微元体 s e 流动的动量对时间的变化率等于外界作用在该微元体上的各种力之和。该定律实际上 是牛顿第二定律。按照这一定律，可导出x 、y 和z 三个方向的动量守恒方程1 1 9 1 ： 百a ( p u ) + d i v ( p u 炉一罢+ 等+ 鲁+ 誓+ 只 ( 2 6 a ) 优 呶o xo va z 百a ( p v ) + d i v ( p v 炉一考+ 鲁+ 鲁+ 誓+ 亿6 b ， 掣d t 油( 炉一挈o z + 孕o x + 肇o y + 鲁+ t ( 2 6 c ) 宓 其中，p 是流体微元体上的压力；k 、z x y 和k 等是因分子粘性作用而产生的作用 在微元体表面上的粘性应力f 的分量；足、毋和f z 是微元体上的体力，若体力只有重力， 且z 轴竖直向上，则f x - - 0 ，f y - - 0 ，f z = - p g 。 式( 2 6 ) 是对任何类型的流体( 包括非牛顿流体) 均成立的动量守恒方程。对于牛顿 流体，粘性应力f 与流体的变形率成比例，有： f 。：2 _ 9 u + 2 d i v ( 甜) f 眇：2 譬+ 2 d i l ，( “) f 荔= 2 掣+ 磁v ( 甜) 龙 当当 ( 2 7 ) f 砂2 f f2 ( 瓦+ 瓦) 得： f 2f “ 气牡2 气斜 其中，是动力粘度，a 是第二粘度，一般可取允= 一；。将式( 2 7 ) 代入式( 2 6 ) ， _ a ( p _ u ) + 破v ( 甜) ：旃，( g r a d u ) 一字+ s 。 ( 2 8 a ) d t 0 ) ) 伽一锄伽一砂 + + 抛一瑟加一瑟 l，k 工 乒 拉管结构绕流的数值模拟及理论分析 掣+ 访v ( p v u ) ：d i v ( , u , g r 础) 一窭+ s ， ( 2 8 b ) 掣+ d i v ( p w 甜) ：咖( g r a d w ) 一罢+ s 。 ( 2 8 c ) o t o z 其中，g r a d o ：票+ 票+ 罂，符号墨、& 和& 是动量守恒方程的广义源项， s 。= c + s ，s ，= e + j ，s 。= e + s ：，而其中的s x 、s y 和s ：的表达式如下： 以：昙( 罢) + 昙。尝) + 昙罢) + 拿一) ( 2 9 a ( 2 ( j u ( 2 d i v u 9 a ) 、 s ，2 i 【i ) + ii ) + ii ) + _ )仫 优劣卯劣必缎o x s ，：昙( 罢) + 晏( 宴) + 昙( 娑) + 晏( 胁甜) ( 2 9 b ) s y2 瓦( 瓦) + 瓦( 瓦) + 瓦( 瓦) + 荀( 黝v 甜) ( 2 s ：昙( 昙) + 晏( 刍+ 昙( 娑) + 拿一一v u ( , ；t d i v u ) ( 2 9 c ) s ：= i 【i ) + i ( 习+ i ( i ) + _ o xo z ( o zo zo zo z 一般来讲，卧曲和s z 是小量，对于粘性为常数的不可压流体，s ，= s ，= s ：= 0 。 式( 2 8 ) 是动量守恒方程，简称动量方程。 ( 3 ) 能量守恒方程 能量守恒方程定律包括有热交换的流动系统必须满足的基本定律。该定律可表述 为：微元体中能量的增加率等于进入微元体的净热流量加上与面力对微元体所做的功。 该定律实际是热力学第一定律。 流体的能量e 通常是内能i 、动能k = 妻 2 + 1 ，2 + w 2 ) 和势能p 三项之和，针对总 能量e 建立能量守恒方程，扣除动能的变化，从而得到关于内能i 的守恒方程。内能i 与温度t 之间存在一定关系，即i = c p t ，其中c p 是比热容。以温度t 为变量的能量守 恒方程为【1 9 】： _ o ( # r 1 3 t - d i v ( p u t ) ：咖( 土鲫刀) + 岛 ( 2 1 0 ) 叫 0 p 其中，邰是比热容，丁为温度，k 为流体的传热系数，研为流体的内热源及由于粘 性作用流体机械能转换为热能的部分，有时简称所为粘性耗散项。 大连理工大学硕士学位论文 ( 4 ) 控制方程的通用表达式 为了便于对各控制方程进行分析，并用同一程序对各控制方程进行求解，现在建立 各基本控制方程的通用形式。 比较四个基本的控制方程，可以看出，尽管这些方程中因变量各不相同，但它们均 反映了单位时间单位体积内物理量的守恒性质。用矽表示通用变量，则上述各控制方程 都可以表示成以下通用形式【1 9 1 ： 掣+ 旃1 ，(