基于大涡模拟方法的多层动网格技术识别平板气动参数

振动与冲击　第３Ｏ卷第４期　ＪＯＵＲＮＡＬ　ＯＦ　ＶＩＢＲＡＴＩＯＮ　ＡＮＤ　ＳＨ０ＣＫ　基于大涡模拟方法的多层动网格技术识别平板气动参数　刘祖军，葛耀君，杨泳昕　（同济大学桥梁工程系，上海２０００９２）　摘　要：从弱耦合的，ｆＪ度出发，对流体计算软件ｆｌｕｅｎｔ进行二次开发，利用其用户自定义函数（ＵＤＦ）描述结构的　运动状态并结合动网格技术实现流固耦合。在保证动网格运动速度符合空间守恒法则的条件下，针对固体模型在流场中　运动受网格尺寸限制且易造成网格变形过大导致计算失败的问题提出了多层动网格的解决方法。流体动力计算时考虑　湍流的作用，采用大涡模拟方法求解Ｎ—Ｓ方程。数值模拟了平板做单自由度强迫振动的断面绕流流场，通过最小二乘　法拟合气动力时程曲线获得气动导数。仿真结果与通过Ｔｈｅｏｄｏｒｓｅｎ理论导出的平板气动导数具有良好的一致性。　关键词：多层动网格；大涡模拟；气动参数；ＵＤＦ；流固耦合；空间守恒法则　中图分类号：Ｕ４４１　文献标识码：Ａ　Ｉｄｅｎｔｉｆｉｃａｔｉｏｎ　ｏｆ　ａｅｒｏｄｙｎａｍｉｃ　ｐａｒａｍｅｔｅｒｓ　ｏｆ　ａ　ｆｌａｔ　ｐｌａｔｅ　ｗｉｔｈ　ｍｕｌｔｉ－ｍｏｖｉｎｇ　ｇｒｉｄ　ｔｅｃｈｎｉｑｕｅ　ｂａｓｅｄ　ｏｎ　ｌａｒｇｅ　ｅｄｄｙ　ｓｉｍｕｌａｔｉｏｎ　ＬＩＵ　Ｚｕ－ｊｕｎ，ＧＥ　Ｙａｏ－ｊｕ￣，ＹＡＮＧ　一ｘｉｎ　（Ｔｏｎ￣ｉ　Ｕｎｉｖｅｒｓｉｔｙ，Ｓｈａｎｇｈａｉ　２０００９２，Ｃｈｉｎａ）　Ａｂｓｔｒａｃｔ：　Ｂａｓｅｄ　ｏｎ　Ｆｌｕｅｎｔ　ｓｏｆｔｗａｒｅ　ｒｅｄｅｖｅｌｏｐｍｅｎｔ，ｔｈｅ　ｍｏｔｉｏｎ　ｏｆ　ａ　ｓｔｒｕｃｔｕｒｅ　ｗａｓ　ｄｅｓｃｒｉｂｅｄ　ｗｉｔｈ　ＵＤＦ　ａｎｄ　ｍｏｖｉｎｇ　ｇｒｉｄｓ　ｗｅｒｅ　ａｌｓｏ　ｕｓｅｄ　ｔｏ　ｒｅａｌｉｚｅ　ｆｌｕｉｄ—ｓｏｌｉｄ　ｃｏｕｐｌｉｎｇ　ｆｒｏｍ　ｖｉｅｗ　ｏｆ　ａ　ｌｏｏｓｅｌｙ　ｃｏｕｐｌｅｄ　ｍｏｄｅ１．Ａ　ｎｅｗ　ｍｅｔｈｏｄ　ｏｆ　ｍｕｌｔｉ－ｍｏｖｉｎｇ　ｇｒｉｄ　ｔｅｃｈｎｉｑｕｅ　ｗａｓ　ｐｒｏｐｏｓｅｄ　ｔｏ　ｒｅｓｏｌｖｅ　ｐｒｏｂｌｅｍｓ　ｏｆ　ｓｔｒｕｃｔｕｒｅ　ｍｏｖｅｍｅｎｔ　ｒｅｓｔｒｉｃｔｅｄ　ｂｙ　ｆｌｕｉｄ　ｍｅｓｈ　ｓｉｚｅ　ａｎｄ　ｃａｌｃｕｌａｔｉｏｎ　ｆａｉｌｕｒｅ　ｃａｕｓｅｄ　ｂｙ　ｌａｒｇｅ　ｄｅｆｏｒｍａｔｉｏｎ　ｏｆ　ｇｒｉｄｓ　ｕｎｄｅｒ　ｔｈｅ　ｃｏｎｄｉｔｉｏｎ　ｔｈａｔ　ｔｈｅ　ｖｅｌｏｃｉｔｙ　ｏｆ　ｄｙｎａｍｉｃ　ｍｅｓｈｅｓ　ｗａｓ　ｉｎ　ｃｏｎｆｏｒｍｉｔｙ　ｗｉｔｈ　ｔｈｅ　ｓｐａｃｅ　ｃｏｎｓｅｒｖａｔｉｏｎ　ｌａｗ．Ｉｎ　ｈｙｄｒｏｄｙｎａｍｉｃ　ｃａｌｃｕｌａｔｉｏｎｓ，ｔｈｅ　ｍｅｔｈｏｄ　ｏｆ　ｌａｒｇｅ　ｅｄｄｙ　ｓｉｍｕｌａｔｉｏｎ　ｗａｓ　ｕｓｅｄ　ｔｏ　ｓｏｌｖｅ　Ｎ—Ｓ　ｅｑｕａｔｉｏｎ　ｉｎ　ｏｒｄｅｒ　ｔｏ　ｃｏｎｓｉｄｅｒ　ｔｈｅ　ｉｍｐａｃｔ　ｏｆ　ｔｕｒｂｕｌｅｎｃｅ．Ｎｕｍｅｒｉｃａｌ　ｓｉｍｕｌａｔｉｏｎ　ｏｆ　ｆｌｏｗ　ｆｉｅｌｄ　ａｒｏｕｎｄ　ａ　ｓｉｎｇｌｅ　ｄｅｇｒｅｅ　ｏｆ　ｒｅｅｄｏｍ　ｆｏｒｃｅｄ　ｖｉｂｒａｔｉｏｎ　ｐｌａｔｅ　ｗａｓ　ｇｉｆｖｅｎ．１ｔｓ　ａｅｒｏｄｙｎａｍｉｃ　ｐａｒａｍｅｔｅｉｓ　ｗｅｒｅ　ａｃｑｕｉｒｅｄ　ｂｙ　ｆｉｔｔｉｎｇ　ａｅｒｏｄｙｎａｍｉｃ　ｔｉｍｅ—ｈｉｓｔｏｒｙ　ｃｕｒｖｅｓ　ｗｉｔｈ　ｔｈｅ　ｌｅａｓｔ—ｓｑｕａｒｅ　ｍｅｔｈｏｄ．Ｔｈｅ　ｒｅｓｕｌｔｓ　ｗｅｒｅ　ｉｎ　ｂｅｔｔｅｒ　ａｇｒｅｅｍｅｎｔ　ｗｉｔｈ　Ｔｈｏｄｏｒｓｅｎ　ｔｈｅｏｒｅｔｉｃ　ｓｏｌｕｔｉｏｎｓ．　Ｋｅｙ　ｗｏｒｄｓ：ｔｈｅ　ｍｕｌｔｉ—ｍｏｖｉｎｇ　ｇｒｉｄ　ｔｅｃｈｎｉｑｕｅ；ｌａｒｇｅ　ｅｄｄｙ　ｓｉｍｕｌａｔｉｏｎ；ａｅｒｏｄｙｎａｍｉｃ　Ｐａｒａｍｅｔｅｒ；ＵＤＦ；ｆｌｕｉｄ—ｓｏｌｉｄ　ｃｏｕｐｌｉｎｇ：ｓｐａｃｅ　ｃｏｎｓｅｒｖａｔｉｏｎ　１ａｗ　目前桥梁风工程中的气动导数一般通过风洞节段　模型试验来获得，采用的方法主要有自由振动法和强　高和计算流体动力学的发展，采用数值模拟方法进行　气动导数识别在断面选型中应用越来越广泛。　迫振动法。前者通过给模型一个初始位移或脉冲激　数值模拟　采用动网格方法识别气动导数如果仅　控制结构表面一层的网格节点运动就要求网格划分得　足够细，结构振幅不易过大，否则容易出现网格节点间　的间距过小，导致网格形成很小的锐角或网格空间位　置发生翻转而产生错误无法计算，而振幅太小不能模　励，使其在气动力和结构阻尼下做自由衰减运动，后者　采用机械装置驱动模型作可控频率、振幅的强迫运动，　测得模型位移、加速度或表面压力等时程，通过各种参　数识别方法获得气动导数。自由振动法试验设备简　单，识别精度低，折减风速范围小，强迫振动法虽然可　以克服这些缺陷，参数识别精度高，能得到大范围折算　拟结构体的真实振动情况，为了兼顾振幅与网格密度，　本文采用多层动网格技术来解决这一问题。为考虑湍　风速下的气动导数，但是由于其试验设备要求高，在工　程中应用并不很广泛　。随着计算机硬件能力的提　流影响采用大涡模拟方式（ＬＥＳ）求解Ｎ—ｓ方程，通过对　流体计算软件ｆｌｕｅｎｔ的二次开发，利用其用户自定义函　数描述结构的运动状态，采用最／Ｊ￣－－乘法拟合气动力　基金项目：国家自然科学基金委重大研究计划（９０７１５０３９；９１０１５０１３）；国　家科技支撑计划（２００８ＢＡＧ０７Ｂ０２）联合资助　时程曲线来识别断面的气动导数。　收稿日期：２０１０—０３—０８修改稿收到日期：２０１０—０７—０５　第一作者刘祖军男，博士生，１９７８年生　第４期　刘祖军等：基于大涡模拟方法的多层动网格技术识别平板气动参数　１５７　１基于多层动网格的气动导数识别　１．１流体动力学计算　流体域内计算结构在一个振动周期内作用在其表　面的气动力是通过将结构表面的流体网格节点的简谐　振动位移等分的离散到／Ｖ个时间步上，将第ｉ个时间　步上变形后的节点坐标作为第ｉ＋１步的网格文件，每　个时间步上求解Ｎ—ｓ方程。　１：■　和Ｌｏｍｂａｒｄ明确提出了用来确定网格速度的法则也称　为空间法则，该法则直接从连续方程中推得：　，｜ｒ　ＪｐｄＶ＋Ｊｐ（Ｍ—Ｍ６）‘ｎｄＳ＝０　式中流体速度为常数，因此几何守恒法则表达为：　，１，　Ｕ　（３）　（４）　ｆ　ｄ　—Ｊ“６・ｎｄＳ＝０　Ｊ　因此只要是动网格问题都必须保证网格速度满足　几何守恒法则，否则可能在动网格计算过程中导致人　Ｐ・　ＯＵ＋Ｐ（ｕ・　）ｕ＝ｐ　—　ｐ＋　（，　△　＋　（　・ｕ）　（１）　－）　式中：Ｐ为流体密度，　为流体速度矢量，　为动力粘性　系数，Ｆ为体积力，Ｐ为压强　为Ｈａｍｉｌｔｏｎ算子，△为　Ｌａｐｌａｃｅ算子。　湍流　的存在使得求解上述Ｎ—ｓ方程变得十分　困难，通过大量学者的不懈努力已提出较多的湍流模　型，较为常见的有：雷诺时均方法（ＲＡＮＳ），直接数值模　拟（ＤＮＳ）以及介于二者之间的大涡模拟（ＬＥＳ）。由于　ＲＡＮＳ方法平均的结果是将时空变化的细节一概抹平，　丢失了包括在脉动运动中的大量信息，并且为了使　Ｎ—Ｓ方程封闭，提出了Ｒｅｙｎｏｌｄｓ应力模型和涡粘模型，　这些模型有一定的局限性，都存在着对经验数据过分　依赖和预报精度较差的特点。ＤＮＳ采用计算机直接数　值求解三维非定常Ｎ—ｓ方程，对湍流瞬时运动直接模　拟，可认为是一种精确的方法。主要用做湍流基础研　究，目前仅限于科学研究领域。　大涡模拟　一　是放弃了对全部尺度范围内的涡运　动进行直接模拟的梦想，改为只对尺度大的涡运动通　过数值方法直接求解Ｎ—Ｓ方程，对尺度小的涡运动不　直接求解而是通过建立模型模拟小涡运动对大涡的影　响。其具体实现过程为：首先通过滤波方法将小尺度　涡和大尺度涡分离开来，即将流动划分为低频可解和　高频不可解（亚格子部分）两部分。对不可解部分建立　亚格子尺度模型来模拟。目前常用的亚格子模型主要　可以分为三类：涡粘模型，相似性模型和混合模型。本　文ＣＦＤ计算时采用Ｓｍａｇｏｒｉｎｓｋｙ—Ｌｉｌｌｙ涡粘形式的亚格　子雷诺应力模型其表达式如下：　７＝ｔｔｉｕｊ—Ｕｉ　ｆ＝２ｖ　Ｓ　７一　一　占　Ｊ　Ｉ　１　２（Ｃ　△）　ｓ　（２５　５　）　一÷丁轰６　Ｊ　（２）　Ｓｍａｇｏｒｉｎｓｋｙ模型相当于混合长度形式的涡粘模　型，其中ｃ　△相当于混合长度，△为过滤尺度，Ｃ　为　Ｓｍａｇｏｒｉｎｓｋｙ常数，Ｌｉｌｌｙ于１９８７年通过计算建议Ｃ　一　０．１８。　１．２动网格空间守恒条件　采用动网格方法必须满足几何守恒条件，Ｔｈｏｍａｓ　为的质量源项。　本文提出的多层动网格方法是通过设定每层动网　格做刚体的平动或扭转运动，相邻动网格层之间通过　弹簧近似方法连接形成变形网格层，变形网格区，网格　只发生变形没有网格速度。因此只需证明发生刚体平　动的动网格层满足几何守恒条件即可。　如图１所示的控制体中，Ｗ，Ｓ，ｅ，凡分别为运动前网　格的四个边界，当网格发生竖向平动时，边界Ｓ平移到　Ｓ　的位置，ｎ平移到ｎ　位置，由于是刚体平移，故ｂｙ　＝　ｂｙ，，边界Ｗ和ｅ没有变化。　【　　１Ｊ　ｌ　　｛Ｊ　Ｉ　（＂＋１）　Ｐ　Ｓ　１　Ａｙ　，　图１　网格运动图示　Ｆｉｇ１．Ｔｈｅ　ｍｏｔｉｏｎ　ｏｆ　ｍｅｓｈ　在图示控制体上，对连续方程进行离散；　＿［　－　Ｉ＿＾　＿ｆ　　＋ｐ［（、　一　　）…　　一　（Ｍ—ｕ　）　］　・（Ａｙ）　＋Ｐ［（　一　）　一　（　一　）　（Ａｘ）”　＝０　（５）　式中：Ｕ，　代表笛卡尔坐标系下的水平速度分量和竖　向速度分量，ｎ和ｎ＋１对应两个速度时刻，ＡＶ为控制　体体积，假定流体运动速度不变，因此有关流体速度项　可以抵消，上式简化为：　＿［　垒　）＿　＿十　＿＿逝十ｐ［（　一　，　，　）．　（Ａｖ）　“＋Ｐ（　一　）（ｂｘ）”　＝０　（６）　按照前面定义网格做竖向刚体运动时：　Ｕｂ，ｅ－０’　＿０，ｚ　Ｊ　６ｌ　＝　代人上式则有：　［（△　）　一（ＡＶ）　）］＝０　（７）　故每层网格的运动都能保证几何守恒法则。　ｌ５８　振动与冲击　２０１１年第３０卷　１．３　多层动网格的实现　动情况，振幅与网格密度不能兼顾，因此需要通过多层　动网格技术来实现。　在工程实际问题中经常遇到流体与固体相互耦合　作用的情况，通常计算域的边界是运动的，从而使计算　域的形状随时间变化。动网格技术随时间变化对流场　多层动网格方法将结构的振动位移按照一定的比　例加到结构周围的Ⅳ层网格节点上，一定程度上加大　了结构的可移动范围，使流体网格在Ⅳ层范围内实现　模型进行更新。在某一时刻，边界发生一定位移后，边　界内部的网格应该如何确定，已经有很多研究成果。　运动。通过设定每层动网格层做刚体的平动或扭转运　动，相邻动网格层之间通过弹簧近似方法连接形成变　形网格层，由于采用多层动网格后相邻网格的平动速　度差值变小，因此变形区的网格不容易出现网格翻转　或过小尖角，易保证网格质量。多层动网格的具体做　较常用的有弹簧近似法，弹性体法，采用Ｌａｐｌａｃｅ方程　控制的方法，其中以弹簧近似方法应用最多。　弹簧近似方法原理简单　（图２），使用方便。在　实际应用中当边界位移较大时，容易出现网格节点问　的间距过小，从而导致网格形成很小的锐角（图３）或　网格点的空间位置发生翻转（图４），导致计算失败。　因此在ＣＦＤ计算时仅控制结构体表面一层网格节点运　动具有一定的局限性。如果要模拟出流体域中的附面　层不同大小的涡时就要将流体的网格划分得足够细，　一法是：按照各层间的比例关系给出各层节点的运动位　移（如图５），结构表面距离为Ｒ　的第ｉ层网格节点的　位移最大等于结构的振幅即Ｓ　＝Ａ　，结构表面距离为　的第ｎ层网格位移为０即Ｓ　＝０，位于这两层之间的　一　４　般结构表面的网格在０．Ｏ０　００１　ｍ的数量级，所以结　第　层网格节点位移进行线性插值　ｉ　＝　ｎ　，这　ｉ　构表面节点空间位置的变化对其周围流体单元的形态　影响较大。结构振幅较大，容易导致网格变形过大，而　产生错误无法计算；振幅太小，不能模拟结构的真实振　样就实现了运动结构周围流体网格的多层运动。这样既　提高了网格可动范围，又减小了相邻两层动网格之间的　相对位移，避免了网格位移过大造成的计算失败。　变形　。　＼　＼　＼　／　／　，　变形后网格位置／　ｌ＼　／　Ｉ　＼　／ｆ＼　＼＼／／　形前网ｔ／／／　格位置　／／／　图２动网格弹簧节点连接　Ｆｉｇ．２　Ｓｐｒｉｎｇ　ｎｏｄｅ　／／　／　、＼＼　＼Ｉ　／　＼１　变形前网格位置　图３　网格变形后出现过小的尖角　Ｆｉｇ．３　Ｐｒｏｄｕｃｅ　ｓｍａｌｌ　ｃｕｓｐ　图４　网格变形后出现翻转　ｇ４．Ｐｒｏｄｕｃｅ　ｆｌｉｐ　ａｆｔｅｒ　ｄｅｆｏｒｍａｔｉｏｎ　ｍｅｓｈ　ｏｆ　ｄｙｎａｍｉｃ　ｍｅｓｈ　ａｆｌｅｒ　ｄｅｆｏｒｍａｔｉｏｎ　ｍｅｓｈ　Ｍ＝　Ｕ２（２Ｂ　）・　［ＫＡ１　ｈ＋ＫＡ２　Ｂｏ　ｇ＋Ｋ２Ａ３　＋Ｋ２Ａ２　ｈ］（９）　理想平板是完全理想化的流线型结构，可以通过　图５多层动网格示意图　Ｆｉｇ　５．Ｔｈｅ　ｍｕｌｔｉ—ｌａｙｅｒ　ｄｙｎａｍｉｃ　ｍｅｓｈ　Ｔｈｅｏｄｏｒｓｏｎ建立的理论公式得到振动平板所受到的气　动力，从而导出气动导数理论值。平板振动时所受到　的非定常气动升力和升力矩可表示为：　Ｌ＝一２ｐ￣ｂＵ２・’　１．４气动参数的识别原理　气动力可以用气动导数来表示，而气动导数是系　统运动折减频率的函数，其函数关系决定于断面的气　动外形，一般很难有解析表达式，大多通过节段模型风　洞试验来识别。气动力用试验测到的８个气动导数表　示为：　：　｛ｃ（　）［　＋　ｈ＋［１＋ｃ（　）］　号］）（１ｏ）　Ｍ＝ｐ＇ｎｂ　Ｕ２・　｛ｃ（　）　ｈ＋［１一ｃ（ｋ）３　Ｔ６　（２Ｂ　）・　（１１）　式（８）一式（１１）中：Ｕ为来流速度，Ｂ＝２ｂ为平板宽　［删　ｈ＋ＫＨ；　＋　埘　＋　百ｈ］　（８）　度；Ｋ：ｗＢ／Ｕ为折算频率；　为振动频率；ｈ，ａ分别为　平板竖向位移和扭转角；　第４期　刘祖军等：基于大涡模拟方法的多层动网格技术识别平板气动参数　１５９　Ｃ（ｋ）：Ｆ（ｋ）＋ｉＧ（ｋ）为Ｔｈｅｏｄｏｒｓｏｎ函数；　平板模型如图７，宽高比２２．５，模型前后缘尖锐，　Ｈ　，Ａ　（ｉ：１～４）为气动导数。　根据Ｔｈｅｏｄｏｒｓｏｎ建立的上述自激力公式可以导出　小攻角下流线形特征明显，流场计算域外边界为矩形，　参考同济大学ＹＪ一１风洞尺寸设为２×６　ｔｎ。　平板断面的８个气动导数，并表达为折减频率函数：　Ａ１＂㈩＝　，Ａ２　）＝　［　＋　］　Ａ３　）＝　Ｆ一　ｋＧ＋　Ａ４　）一　图７平板模型图　Ｆｉｇ　７．Ｔｈｅ　ｍｏｄｅｌ　ｏｆｆｌａｔ　ｐｌａｔｅ　Ｈ　）＝一　，Ｈ２　）＝一　［　＋Ｆ＋　］　（　）＝一　［Ｆ一　］，　（　）＝詈［÷＋　Ｇ］　（１２）　强迫平板单独作竖向正弦振动或扭转正弦振动，　将结构表面压力值和摩阻力沿表面积分，可得到在绝　对坐标系里的气动升力和弯矩。然后由得到的气动力　和式（８）一（９）中对应的气动导数之系数，根据最小二　乘法，即可得到相应的气动导数。竖向振动时采用上　述方法可以提取　，磁，４　，　，４个气动导数，做扭　转运动可以确定另外的４个气动导数彤，埘，４　，Ａ　。　２平板的颤振导数识别　通过流体分析软件Ｆｌｕｅｎ　的用户自定义函数来　描述模型的运动状态，并利用其大涡模拟的湍流模式　进行数值计算（图６）　图６数值计算流程　３　０　２　０　】　０　０　８　—１　—０　—２　．３　：８　０　６．０　６．２　６．４　６．６　６．８　７．０　６．０　６．２　６．４　６．６　６．８　７．０　ｅＵ　ｔＵｆＢ　图８低折算风速下纯竖向运动产生的　气动升力系数和升力矩系数时程　Ｆｉｇ．８　Ｌｉｆｔ　ａｎｄ　ｍｏｍｅｎｔ　ｅｏｅｆｆｉｃｉｅｎｔｓ　Ｖｓ　ｔｉｍｅ　ｏｎ　ｈｅａｖｉｎｇ　ｖｉｂｒａｔｉｏｎ　ａｔ　ｌｏｗｅｒ　ｒｅｄｕｃｅ　ｗｉｎｄ　ｓｐｅｅｄ　平板竖向振动的振幅０．０８ｍ，扭转幅角３。。人口　为速度边界条件，人口风速为Ｕ＝５　ｎｒ／ｓ，出口为压力　边界条件，表压设为０，参考压力点在出口边界中心处；　上下边界设为对称边界，以加快计算速度，时间步长设　为０．００１　Ｓ。建立二维计算模型，网格线与物面平行，周　向网格线在靠近物面附近加密，离开物面间距逐渐加　大。交于物面的网格线在模型前后缘位置均加密．网　格数为３．０２万。　模型从零攻角的平衡位置开始运动．图８、图ｌ０是　低折算风速（　＝５）下的运动，图９、图１１是高折算风　速（Ｖｒ＝２０）下的运动．平板分别作纯竖向振动和纯扭　转振动，气动力数据在２个周期后开始采集。从图８一　图１１可知，在这两种振动情况下，升力系数和升力矩　系数均表现为对无量纲时间零均值的一次稳态谐波曲　线特性，升力、升力矩与位移之间的良好的线性关系表　明流场与平板一起进入谐振动状态，流场没有出现不　稳定情况。　在不同的折算风速下，分别计算平板做纯竖向振　动、纯扭转振动的气动力，当气动力的变化趋势变得平　稳时，开始采集气动力数据，并用最最小二乘法提取气　动导数。本文比较了计算得到的气动导数与通过Ｔｈｅ．　ｏｄｏｒｓｅｎ理论¨叫推导出的气动导数之间的差别，见图　ｌ２，从图中各颤振导数值的比较来看，当折减风速较低　时，各导数值吻合较好，随着折减风速增大，各导数值　偏差增大，其中　和Ａ　的误差相对大些，但均与通　过Ｔｈｅｏｄｏｒｓｅｎ理论导出的气动导数值有合理的一致　性，其中几个关键导数Ｈ　，埘，　，　及　吻合的　很好　Ｏ　Ｏ　Ｏ．２０　Ｏ　０．１５　０　Ｏ．１Ｏ　０　０．０５　一Ｏ　０．００　—．Ｏ　０．０５　—０＿｝０　—Ｏ　．０．１５　．Ｏ　．Ｏ．２０　８　０　９．０　１Ｏ．Ｏ　ｌ１．Ｏ　ｆ（　Ｂ　图９高折算风速纯竖向运动产生的　气动升力系数和升力矩系数时程　Ｆｉｇ．９　Ｌｉｆｔ　ａｎｄ　ｍｏｍｅｎｔ　ｃｏｅｆｉｆｃｉｅｎｔｓ　Ｖｓ　ｔｉｍｅ　ｏｎ　ｈｅａｖｉｎｇ　ｖｉｂｒａｔｉｏｎ　ａｔ　ｈｉｇｈｅｒ　ｒｅｄｕｃｅ　ｗｉｎｄ　ｓｐｅｅｄ　１６０　振动与冲击　２０１１年第３０卷　０　———０　０　０　０　０　０　０　０　０．３　０．２　０．０８　０．１　０．０　—０　１　０．２　０．３　一０．０２　—：　，．０．０８　图１０低折算风速纯扭转运动产生的　气动升力系数和升力矩系数时程　Ｆｉｇ．１０　Ｌｉｆｔ　ａｎｄ　ｍｏｍｅｎｔ　ｃｏｅｆｆｉｃｉｅｎｔｓ　Ｖｓ　ｔｉｍｅ　ｏｎ　ｐｉｔｃｈｉｎｇ　ｖｉｂｒａｔｉｏｎ　ａｔ］ｏｗｅｒ　ｒｅｄｕｃｅ　ｗｉｎｄ　ｓｐｅｅｄ　图１１　高折算风速纯扭转运动产生的　气动升力系数和升力矩系数时程　Ｆｉｇ．１ｌ　Ｌｉｔ　ａｎｄ　ｍｏｍｅｎｔｆ　ｃｏｅｆｉｃｉｆｅｎｔｓ　Ｖｓ　ｔｉｍｅ　ｏｎ　ｐｉｔｃｈｉｎｇ　ｖｉｂｒａｔｉｏｎ　ａｔ　ｈｉｇｈｅｒ　ｒｅｄｕｃｅ　ｗｉｎｄ　ｓｐｅｅｄ　Ｏ　一１　２　３　——５　２０　Ｏ．８　Ｏ．６　一　Ｏ．４　０．２　图１２平板气动导数　Ｆｉｇ．１２　Ｆｌｕｔｔｅｒ　ｄｅｒｉｖａｔｉｖｅｓ　ｏｆ　ｌａｔｆ　ｐｌａｔｅ　３　结　论　本文通过对流体计算软件ｆｌｕｅｎｔ的二次开发，利用　ｆｏｒ　ｆｌｏｗ　ｏｖｅｒ　ａ　ｓｅｌｆ－ｅｘｃｉｔｅｄ　ｖｉｂｒａｔｉｎｇ　Ｂｏｄｙ［Ｊ］．Ｔｓｉｎｇｈｕａ　Ｓｃｉｅｎｃｅ　ａｎｄ　Ｔｅｃｈｎｏｌｏｇｙ，１９９９，４（４）：１６８８—１６９１．　［４］Ｆａｓｅｌ　Ｈ　Ｆ，Ｓｅｉｄｅｌ　Ｊ，Ｗｅｒｎｚ　Ｓ．Ａ　Ｍｅｔｈｏｄｏｌｏｇｙ　ｆｏｒ　ｓｉｎｍｌａｔｉｏｎｓ　ｏｆ　ｃｏｍｐｌｅｘ　ｔｕｒｂｕｌｅｎｔ　ｆｌｏｗｓ　『Ｊ］．　Ｊｏｕｒｎａｌ　ｏｆ　Ｆ１ｕｉｄｓ　Ｅｎｇｉｎｅｅｒｉｎｇ，２００２，１２４（５）：９３３—９４２．　其提供的用户自定义函数模块，结合动网格技术实现　了平板的气动导数识别。由于固体模型在流场中运动　受到流体网格尺寸限制，当模型运动幅度较大时，常会　造成网格的变形过大，导致计算无法进行。针对这一　问题，本文采用了多层动网格方法较好地解决了这一　问题。能较好的兼顾模型振幅和网格密度。采用本文　［５］Ｋａｈｔｅｎｂａｃｈ　Ｊ　Ｈ，Ｓｃｈｕｍａｎｎ　Ｕ，Ｇｅｒｚ　Ｔ．Ｌａｒｇｅ　ｅｄｄｙ　ｓｉｍｕｌａｔｉｏｎ　ｏｆ　ｔｕｒｂｕｌｅｎｔ　ｄｉｆｆｕｓｉｏｎ　ｉｎ　ｓｔａｂｌｙ—ｓｔｒａｔｉｉｆｅｄ　ｆｌｏｗ　Ｊ］．Ｆｌｕｉｄｓ　Ｍｅｃｈ，１９９４，２８０：ｌ一４Ｏ．　［６］Ｋｏｎｄｏ　Ｋ，Ｍｕｒａｋａｍｉ　Ｓ，Ｍｏｃｈｉｄａ　Ａ．Ｇｅｎｅｒａｔｉｏｎ　ｏｆ　ｖｅｌｏｃｉｔｖ　ｌｕｃｔｆｕａｔｉｏｎｓ　ｆｏｒ　ｉｎｆｌｏｗ　ｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎ　ｏｆ　ＬＥＳ，『Ｊ］．Ｊｏｕｒｎａｌ　，方法所获得的计算结果与通过Ｔｈｅｏｄｏｒｓｅｎ理论导出的　平板气动导数具有良好的一致性。　参考文献　ｏｆ　Ｗｉｎｄ　Ｅｎｇｉｎｅｅｒｉｎｇ　ａｎｄ　Ｉｎｄｕｓｔｒｉａｌ　Ａｅｒｏｄｙｎａｍｉｃｓ１９９７，６７：　５１—６４．　［７］Ｔｅｔｓｕｒｏ　Ｔａｍｕｒａ，Ｙｏｓｈｉｙｕｋｉ　Ｏｎｏ．ＬＥＳ　ａｎａｌｙｓｉｓ　ｏｎ　ａｅｒｏｅｌａｓｔｉｃ　［１］Ｓｃａｎｌａｎ　Ｒ　Ｈ，Ｔｏｍｋｏ　Ｊ　Ｊ．Ａｉｒｆｏｉｌ　ａｎｄ　ｂｒｉｄｇｅ　ｆｌｕｔｔｅｒ　ｄｅｒｉｖａｔｉｖｅｓ　［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｅｎｇｉｎｅｅｒｉｎｇ　Ｍｅｃｈａｎｉｃｓ　Ｄｉｖｉｓｉｏｎ，１９７１，９７：　１７１７一ｌ７３７．　ｉｎｓｔａｂｉｌｉｔｙ　ｏｆ　ｐｒｉｓｍｓ　ｉｎ　ｔｕｒｂｕｌｅｎｔ　ｆｌｏｗ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｗｉｎｄ　Ｅｎｇｉｎｅｅｒｉｎｇ　ａｎｄ　ＩｎｄｕｓｔｒｉａｌＡｅｒｏｄｙｎａｍｉｃｓ，２００３，９１：５ｌ一６４．　［８］郭正，李晓斌．用非结构动网格方法模拟有相对运动的　［２］Ｙａｎｇ　Ｙｏｎｇｘｉｎ，Ｇｅ　Ｙａｏｊｕｎ，Ｘｉａｎｇ　Ｈａｉｆａｎ．Ｉｎｖｅｓｔｉｇａｔｉｏｎ　ｏｎ　ｌｆｕｔｔｅｒ　ｍｅｃｈａｎｉｓｍ　ｏｆ　ｌｏｎｇ—ｓｐａｎ　ｂｒｉｄｇｅｓ　ｗｉｔｈ　２ｄ．．３　ＤＯＦ　ｍｅｔｈｏｄ．ＷＩＮＤ　ＡＮＤ　ＳＴＲＵＣＴＵＲＥＳ，多体绕流［Ｊ］．空气动报，２００１，１９（３）：２８—３５．　［９］王福军．计算流体动力学分析——ｃＦＤ软件原理与应用　［Ｍ］．北京：清华大学出版社，２００４．　［１　Ｏ］Ｔｈｅｏｄｏｒｓｅｎ　Ｔ．Ｇｅｎｅｒａｌ　ｔｈｅｏｒｙ　ｏｆ　ａｅｒｏｄｙｎａｍｉｃ　ｉｎｓｔａｂｉｌｉｔｖ　ａｎｄ　ｔｈｅ　ｍｅｃｈａｎｉｓｍ　ｏｆ　ｌｆｕｔｔｅｒ［Ｒ］．ＮＡＣＡ　Ｒｅｐｏｒｔ，ｌ９３５．　２００７，１０（５）：　４２１—４３５．　［３］Ｃｈｅｎ　Ｚ，Ｈｏｎｇ　Ｌ，Ｘｉａｏｆｅｎｇ　Ｗ．Ｆｕ１１　ｆｕｎｃｔｉ０ｎ　ｎｕｉｎｅｒｉｃａｌ　ｍｅｔｈ０ｄ