（流体力学专业论文）机翼跨音速颤振的频域计算.pdf

南京航空航天大学硕十学位论文 摘要 本文研究了机翼跨音速颤振的频域计算方法。以给定机翼模态分布下机翼j 二各点 的模态值作为运动幅值，以三维非定常全位势方程和三维非定常e u l e r 方程为控制方 程，求解三维机翼简谐运动下的非定常气动力。剐全位势方程采用精确隐式近似因式 分解差分格式、牛顿内迭代算法。对g u l e r 方程采用有限体积法和双时间推进。所求 得的气动力作为已知值运用于颤振方程的求解。而颤振方程仍采用传统频域解法中的 v - g 法进行求解。对得到的一系列的阻尼、速度和频率进行了线性插值，从而得到颤 振速度和颤振频率。三维非定常气动力和颤振的计算结果是令人满意的。 关键词：非定常、全位势方程、欧拉方程、颤振、频域 塑望堕童鲨塑堡堕鉴茎生簦一一 a b s t r a c t c o m p u t a t i o no fw i n g f l u t t e rf o ru n s t e a d yt r a n s o n i cf l o w si nt h ef r e q u e n c yd o m a i ni s s t u d i e di nt h i st h e s i s t h ew i n gi sm o d e l e ds t r u c t u r a l l yb yv i b r a t i o n a lm o d e s ，w h i c ha r e u s e da sa m p l i t u d eo fm o t i o nb yaf i n i t ee l e m e n ta n a l y s i s c o n s e r v a t i v ef i f l l p o t e n t i a l e q u a t i o na n dt h r e ed i m e n s i o n a lu n s t e a d ye u l e re q u a t i o n a r eu s e dt oc o m p u t et h eu n s t e a d y t r a n s o n i cf l o w sa r o u n da i r f o i la n dw i n g su n d e rs i m p l eh a r m o n i cm o t i o n t h ef u l lp o t e n t i a l e q u a t i o ni ss o l v e db yt h et i m e a c c u r a t ea p p r o x i m a t ef a c t o r i z a t i o na l g o r i t h m a n di n t e r n a l n e w t o ni t e r a t i o n s ，a n dt h ee u l e re q u a t i o nb yt h ef i n i t e v o l u m er e s o l u t i o ns c h e m ea n d d u a l t i m es t e p p i n g t h ec a l c u l a t e da e r o d y n a m i cf o r c e sa r ei n p u tt os o l v ef l u t t e re q u a t i o n s w h i c ha r et r e a t e db yv - gm e t h o do ff r e q u e n c yd o m a i n t oas e r i e so fd a m p i n g s ，v e l o c i t i e s a n df r e q u e n c i e s ，t h el i n e a r i n t e r p o l a t i o n i s a p p l i e d a n dt h ef l u t t e r s p e e d a n df l u t t e r f r e q u e n c ya r eo b t m n e d t h ec a l c u l a t e d r e s u k so ft h r e ed i m e n s i o n a lu n s t e a d yf l o wm a d f l u t t e rh a v e p r o v e dt h a tt h e m e t h o d sp r e s e n t e di nt h et h e s i sa r es a t i s f a c t o r y k e yw o r d s ：u n s t e a d yf l o w , f u l l p o t e n t i a le q u a t i o n ，e u l e re q u a t i o n ，f l u t t e r , f r e q u e n c y 堕塞：堕至堕盔叁堂堡：! 堂堡笙苎一一一 1 1 问题的背景 第一章绪论 飞机颤振是飞机设计过程中经常遇到的一个技术难点。飞机飞行中如果发生颤振 就会引起飞机结构的破坏，危及b 行员和乘客的生命。所以在飞机设计时，在飞行包线 范围内不允许发生颤振。颤振研究涉及到空气动力学和结构动力学问题，涉及面广， 技术难度大。 颤振实验需要大量的资会，一个颤振实验的费用是普通气动实验的许多倍，而如 果采用数值计算来补充风洞实验，那么整个飞机的设计费用就会大大降低。 对颤振方程的数值解法分为时域和频域两种解法。时域法”是在给定初始扰动下 计算机翼随时间的运动变化情况以判断收敛和发散。具体的说，已知某一个时刻f 的 运动，通过非定常气动力方程和结构运动方程的耦合求解求得下个时刻h 出的运 动。通常，出必须足够的小。而厂义坐标的时i n l l ：, 应可以判断出飞机是否发生了颤 振。与时域法不同的另一方法就是在频域巾处理问题。与时域法相比，当前这种方法 更受欢迎，得到了更广泛的应用。其最重要的原因或许在于简谐运动的气动力理论远 比任意非定常运动的发展得更彻底这。事实。 非定常气动力的计算是气动弹性分析计算中的最重要部分。而准确求解跨音速非 定常气动力口o 又是目前解决气弹和颤振问题【4 5 1 的关键。对线性问题，早在7 0 年代 就发展了亚音速偶极子格网法 6 一”求解颤振中的气动力。其基本思想为：1 将翼面沿 展向和弦向分成若干梯形网格，并认为每个网格上的压力差是常数。2 网格的侧缘必 须平行于来流，而且前后网格的侧缘必须在同一条顺来流方向的直线上，不能互相错 丌。翼面前后缘如有弯折，则弯折点必须位于网格角点。如果有操纵面，则操纵面边 缘必须同网格边缘一致。对网格的这些要求都是为了保证前面提到的假定每个网格上 的压力是常数。3 只在每个网格中切面的3 4 弦长点( 称为控制点) 满足积分方程。 这样，积分方程就成为以每个网格的压力差为未知量的线性方程组。4 认为压力偶极 子只分布在网格1 4 弦长线上，并且认为，网格上的压力作用点位于网格中切面的1 4 弦长点，成为压力点。在以上工作的基础上再求解气动力。从此可以看出，偶极子格 网法在亚音速时，求解气动力准确、经济，计算时间短。但由于此方法基于线形理沦， 在颤振经常发生的跨音速阶段，往往采取极为粗糙的凹坑余量扣除。但各种飞机的跨 机冀跨音速勤振的频域计算 音速颤振形态很不千h 同，有的甚至不出现凹坑的情况。这样的余量的扣除已不再适应 现代飞机设计的需要。由此町知，偶极子格网法的准确性不够。 存国外，9 0 年代前后开始发展以e u l e r p 川和n s 方程为丰管方程的跨音速非定 常气动力1 1 1 , 1 2 , 13 1 的求解方法，以及气动弹性分析方法，并应用于三维机翼的颤振特性的 研究。国内主要采用风洞实验研究，但关丁非线性气动力计算也丌始了尝试，并正获 得不断进展。但发展基于非线性。e 动方程的跨音速颤振计算方法仍有许多工作要做。 跨音速非定常气动力i t 算方法分为小扰动方程、全位势方程 j 4 1 、欧拉方程 1 5 , 1 6 , 17j 及n s 方程m 解法。由他们得到的非定常气动力和结构运动方程1 1 9 2 0 , 2 1 1 耦合求解i ：”， 在时域内计算广义坐标的响应以判断飞机是否发生颤振，此方法即为时域法 2 3 川。最 近几年来，国内外已经发表了不少有关这方面的研究报告和学术论文1 2 5 , 2 6 。由上述的 各种解法求得的非定常气动力也町以用于频域法的颤振计算。这种计算颤振的方法目 前国内还很少有人做过，本文正是针对与此，研制了求解颤振的非定常气动力的频域 解法，并编制了相应的程序。 1 2 本文的工作 本文针对前面所提的问题，研究了机翼跨音速颤振的频域计算方法。以三维非定 常全位势方程和三维非定常e u l e r 方程为控制方程，求解三维机翼简谐运动下的非定 常气动力，在此基础上运用了频域法求解颤振。对全位势方程采用精确隐式近似因式 分解差分格式、牛顿内迭代算法。对e u l e r 方程采用有限体积解法和双时问推进。由 于所求气动力是用于频域法的颤振计算，所以对两个方程的非定常气动力的求解过程 进行了相应的修正。以给定机翼模态分布f 机翼上各点的模态值作为运动幅值求解简 谐运动下的非定常气动力，并把得到的非定常气动力作为已知运用于颤振计算。对颤 振方程运用频域解法中的v - g 法进行求解。对得到的一系列的阻尼、速度和频率进行 了线性插值，从而得到了颤振速度和颤振频率。对非定常气动力和颤振方程的求解过 程进行了有效的结合，对此编制了程序。具体工作有如下三部分：非定常全位势方程 的求解、非定常e u l e r 方程的求解、运用频域法求解弹性运动方程。在以上工作的基 础上，研制了“机翼跨音速颤振频域解法的计算程序”。在翼型和不同的三维机翼匕 得到了f 确的计算结果，说明发展的方法和程序的合理性，可以为国内飞机型号设计 部门提供新的颤振计算分析手段。 南京航空航天人学硕十学位论文 2 1 引言 第二章三维非定常全位势方程的求解 颤振计算极大地依赖于非定常。t 动力的计算精度。跨音速时，颤振边界出现跨 音速凹坑，颤振速度有较大下降，往往限制了飞行包线。因此，准确地确定跨音速 颤振速度在飞机设计中有很重要的作用。但线性理论不能预示凹坑现象，必须用 非线性理论来研究跨音速非定常气动力。 基于跨音速小扰动位势方程的有限差分解法最早受人们重视，目前已成熟。用欧 拉方程甚至n s 方程来计算跨音速非定常气动力的研究工作也已经开展。全位势方程 方法由于精度和适用范围l c d , 扰动方法前进一步，计算机时又比欧拉方程方法小一个 数量级，在工程实践中仍有实用价值。 本文采用c h 型有限差分网格”，用时问精确隐式分解筹分格式 2 8 , 2 9 1 、牛顿内 迭代算法求解非定常全位势方程。边界条件及非定常尾迹条件均使用隐式嵌入处理技 术以提高计算效率。鉴于粘性对激波位置有较大影响，解准二维附面层积分方程得到 附面层位移厚度，采用粘位流迭代的方法计及粘性影响。 2 2 控制方程及差分解法 贴体坐标下非定常全位势守恒型方程可写为 ( p ，) 。+ ( p u j ) f + ( p v j ) 。+ ( p w j ) f = o( 2 1 ) 式中p 为密度，u 、矿和是逆变速度分量，1 为时间( - - t ) ，为变换j a c o b i a n 。而 。= t 一孚蛾( z + , + q2 - 1 ) r q 2 5 ( u 一 。) + ( 矿一t 1 ，) $ 。+ ( 一 ，) 十( u 2 f + a | | + + a 1 2 十1 + a l3 q b c 矿= ”j + 。】2 中 + a 2 2 巾n + 口2 3 咖( 机冀跨音速颤振的频域计算 w = r - i - 口1 3 中 + a 2 3 4 ) | 1 十巾c 其中妫速度位，a ，至口。为度规系数。 方程( 2 1 ) 可表示成 f ( m ) = 0 其中十为当前时问层上待求速度位值。用牛顿迭代法解之 砌，+ 州剐 亿：， 式中十为当前时间层上的已知值或假设值。迭代收敛时蛳：十一十。将趋于零。由于 篝只对收敛进程重要，其精度对最后收敛结果影响不大，其中的塞项只保罔在时 间项内。采用近似因式分解，从方程( 2 2 ) i ：1 。uj ， 毛于蛳的线性方程 上 上_ l f = 尺慨)r 2 3 1 其中 名2 1 + 记一去名妇。名 厂 三。= - + t v o - 吉帅“a 。 _ 1 + t 聊一言吐p “乜，吼 而p “为产生人工粘性的逆风密度m 】，r 为方程残值。 ( 2 3 ) 式分三步求解 a 三f 妒= r b - 口。= ” c f = 矽。( 2 4 ) 每步均为一三对角方程组。求解一遍后所得m 新值为 = + + 2 3 网格生成 4 南京航宅航天人学硕十学位论文 2 3 1 代数生成方法 采用无限内插法实现二维网格生成 令： i , l i “u 口，1 4 2 w w d 代表在区域上定义的二变量u ，v 的函数，这个隔数只有部分是已知的： f ( u 女，w ) = a k ( w ) ；k = 1 ，2 p ( “，w 女) = b 。 ) ；k = l ，2 r 于是： 厂( “，w ) = ( w ) ，( “，o ) + ( w ) - o 。- 誓- f v ( “，o ) + 畸“- 厂( “，1 ) 函数f 定义了从u ，w 到x ，z 的变换。 2 3 2 解椭圆型方程的网格优化 设物理空间为直角坐标系( x ，z ) ，计算空间为( 亭，0 。在代数方法生成的网格避础 上按剖面采用线松弛解法解下列l a p a c e 方程进行嘲格优化 。+ ：= 0 叮。+ ：0 ( 2 6 ) 2 3 3 物面网格正交化及压缩处理 以给定的物面上第一层网格高度沿曲线进行压缩，经物面正交化处理得到准三维 网格。 2 3 4 动态网格 5 机宾跨音速颤振的频域计算 为了避免方程复杂化，非定常计算州采用远场周定、物面随动的动态计算网格。 由于每一时阳j 步重新生成汁算劂格过于费时，在主要用于气动弹性计算的小幅运动情 况f ，普遍采用以上时间步为基础解静平衡方程来得到新的计算网格。但这方法 仍不经济。 考察绕圆柱的_ 维4 i 可压位流，令圆柱半径为a ，柬流速度圪，则流场流函数 为： y = 亿卜昝n 口 b ， 放大圆柱半径至a 十s ，则过( r ，0 ) 的原流线将产生径向位移ar ，经推导可得： a r = g - g - 蟛 r + a ( 2 8 ) 这一结论直接应用于动态网格生成较为困难。经进一步的简化处理，三维可压非定常 流的动态网格点坐标( 以下标u 表示) 可用下式确定： 置= i ，( j ，一i 。) g( 2 9 ) 其中下标s 表示作为初值的静态网格点坐标，下标r 表示全流场静态网格点随物面作 刚性运动的瞬间坐标值。而g 值设为网格点序号的函数 其中下标b 表示对应的物而点，下标f 代表远场点。 2 4 边界条件 一a x 2 ，( 鬻 2 ，( 剖2 仁，。， 生成贴体网格后，l 方向选为翼剖面表面曲线方向，t 1 为展向方向， 为由内向外 的径向。这样，尽管方程( 2 4 a ) 、( 2 4 b ) 的未知量物理意义无明确含义，但所涉及的边 界条件仅为远场边界条件和对称边界条件，不难加以实现。 方向的方程( 2 4 c ) 可分为 一培口分求解，即应用物面边界条件的物面区及应用尾迹条件的自由尾迹区。 2 4 1 隐式物件面边界条件 6 南京航空航犬人学硕十学位论文 翼面上的法向速度为零，即w - o 。令i 、j 、k 分别代表、t 、 o ) 。它有两个自由度，既机翼的丘f 平移h ( 向下为正) 及转轴e 的转 动( 转角为口，抬头为正) ，它们又分别称为浮沉运动及俯仰运动。 于是机翼在e 点后r 处的某点的位移可表示为： z = + r 口 ( 5 1 ) 用拉格朗日方程导出颤振运动方程。对于一个可用离散的广义坐标 q ，( f = 1 2 ，h ) 来表示其状态的系统，有： 南京航空航大人学硕十学f 口论文 i d i 8 t ) 一挈+ 婺：q ( h 厶 )( 5 2 ) 万瓦一百+ 瓦2 q ( 7 2 1 厶”)喊到 其巾t 为系统的动能，0 为位能，q 为与吼相对应的广义力。( j 。2 ) 式即为拉格朗同 方程。 一元机翼的动能为： 丁= 圩( 籼融 ：_ ，lm h 2 + s 。矗+ 三，。砬2 ( 5 3 ) 1 d 一 ”， 其中： m 2j 丽出 ( j 3 a ) s 。j 丽出= m x 。b ( 5 3 b ) 乞= f b r 2 融= 硎名2 b 2( 5 3 c ) 图5 1二元棚翼图 以上各式中，丽为单位面积的机翼质量，卅为单位展长的机翼质量。& 为单位 展长机翼对转轴e 的质量静矩，x 。为重心到转轴e 的距离的无量纲量( 重心在后为 ! ! 篓堕童鲨墅塑塑塑塑生篓l 正) 。，。为单位展氏机翼对转轴l 的质量惯矩，0 为机翼对转轴e 的回转半径的无量 纲量。 弹能u 为： u ：当k 。矗：+ 去世。2 ( 5 4 ) ” “ 应用拉格朗f 方程，得 m h + s 。a 七kh h = q h s 。h + ，。舀+ k 。d = 级 ( 5 5 ) 其中，g 为与矗相应的广义力，即振动引起的气动力l ( 向下为正) 。绒为与d 相应 的广义力，即振动引起的对f 点的气动力矩m ，( 抬头为正) 。上式即为二元机翼的 颤振运动方程。 5 3 三元机翼的颤振 5 3 1 概述 三元机翼比二元机翼复杂的多。三元机翼是一个复杂的连续的弹性体，它的自由 度数有无限多个。因此，对于真实机翼，要完全精确的研究它是不可能的。因此，要 研究三元机翼的颤振，首先要把真实的、复杂的机翼简化成个能反映其主要弹性及 惯性特性的比较简单的物理模型以及相应的数学模型。 通常，u j 以用有限元的方法把连续的结构离散成若干部分，而用结构的有限个 离散点的广义位移来近似地描述一个有无限多自由度的变形。例如，在研究梁的弯曲 时，可以把梁分成n 段，把每段的质量集中在各段的中点上。这样，这个梁便成为一 具有n 个集中质量的无自重的弹性梁。粱的变形便可以用这n 个集中质量的位移来描 述。这便简化成一n 个自由度的系统。对于机翼，则可把机翼离散成若干个有限元， 而用若干结点的位移来描述机翼的变形。但是，为了较精确地描述结构的变形，通常 这样简化的系统的自由度数仍然是很高的。 5 3 2 颤振运动方程 2 4 南京航空航灭人学硕士学位论文 个n 个自由度的弹性、惯性系统，具有n 个自振频率珊，和白振型1 ，： ( i = 1 , 2 ，n ) 。若对应每一自振型f ，选取一广义坐标q ，则利用自振型的正交 性，w j 以把系统的自由振动的运动方程，表为n 个独立的单自由度系统的运动方程。 这就是况，这样选取广义坐标，便把一n 个自由度互相耦合的系统变为n 个单自由度 系统| e | 勺迭加。而振动时，各单自由度系统对振动的响应大小取决于自振频率与振动频 率的关系。它们愈靠近，则响应愈大。反之，则响应愈小。而刈实际工作有兴趣的颤 振频率都与比较低的几个自振频率比较靠近。因此，著对应自振型选取广义坐标，我 们便可以只选取与若t 二低阶自振型相应的广义坐标，并以它们来描述系统的变形。 这样机翼的变形可表示为f 1 个自振型的迭加，即 上 z ( x ，y ，：，) = i ( x ，y ，z ) q 心) ( 5 6 ) ，= i 显然，月 n 。如前所述，实际的振型都用若十节点( 例如n 个节点) 的位移来表示， 即t ，：通常表示为一( n ，1 ) 阶的向量。于是，e 式可写成矩阵形式如： z i = 巩( ，) ( 5 7 ) 其中， z ， 都是( n ，1 ) 阶向量。 因为 是自振型，所以它满足自由振动运动方程。即 0 ) 2 阻u = k 协)( f - 1 ，2 ，胛) 其中阻 为质量矩阵，k 】为刚度矩阵，它们都是( n ，n ) 阶对称矩阵。 自振型还具有正交性，即 l ：且肛蟛 溺 及 u 辩船j 2 ( 5 8 ) ( 5 9 ) 曷 ( 51 0 ) 其中，m ，k ，分别是与第i 个自振型相应的广义质量和广义刚度。 当机翼的变形可用自振型的迭加来表示时，机翼的动能为： 2 5 机疆跨音速蓟振的频域计算 丁= 窆k ( x ，y ) 姒f ) 2 而( z ，y ) 出砂 ( 5 其中，而为单位面积机翼的质量。积分区域s 表示整个机翼面积。上式的矩阵形式 为： ，：丢窆主j 阻】矗，口 1 ，l = 去m ，口? ( 5 1 2 ) 上式变换中，引用了自振型的正交性。 弹能为： ，2 圭喜，( l y h ( 。露( 墨儿g 朋喜乃( g 劫钆( ，) 由咖出方 ( 5 13 ) 其中，露( x ，y ，q ，1 1 ) 为刚度影向系数。其矩阵形式为： u ：i 1 乞n n j 医 ) = 去m ，2 9 ，2 ( 5 1 4 ) 将( 5 1 2 ) ，( 5 1 4 ) 代入拉格朗几方程得颤振运动方程： m ，口，+ m ，c o ，| q ，= q ( f _ 1 2 ，月)( 5 1 5 ) 由上式可见，对与惯性力和弹性力项，各自由度问是不耦合的。但q 却与各广义坐 标有关，故上式需要联立求解。 机翼上的广义气动力可用矩阵的形式表示为川： i 1 2 b 】白 其中 q l ( 聊，n ) = k 必】 c f 扣”为第m 个自振型的力作用点的法向位移，c 岁为第n 个自振型的c 。分布，s 为 南京航空航天人学硕十学位论文 l 酬格的面积。为了书写的方便，以卜的q | 伞部去掉卜标1 。 5 4 颤振方程的量纲分析 f 5 1 6 ) 在求解颤振方程之前，先对其进行虽纲分析m j ，这早用的是与v g 解法对应的 颤振运动方程，如f ： t 珊2 m + ( 1 + 喀) k r 1p v 2 爿 g = 0 ( 5 1 7 ) 其中国是颤振频率，1 秒： m 是广义质量，认为振动形态是无量纲量，m 具有质量的量纲，千克： k 是广义刚度，认为振动形态是无量纲量，k 具有刚度的量纲，牛顿米； p 是结构阻尼系数，无量纲： p 是大气密度，千克米3 ： y 是飞行速度，米秒； b 是参考长度，米： a 是，5 义空气动力系数，a = 厂s ( ，厂为振动形态是无量纲的，a s 为网 格的面积，米2 ，c p 一_ a j ，x 具有长度的量纲，所以a 的量纲为米： a x q 是广义坐标，米。 引用符号k 表示模型的。谚数。同实物的谚数x 。之比 型的运动方程可以用实物参数写为： k q q 一。2 尼。，。2 m 。+ ( 1 + i k 。g 。) 女e k 。 一妄女p k ，2 。k 。p 。圪2 b 。a 。 = 0 k 一= k t k “k 。 由此可圃，如粜 即= _ 。i x 。，则模 ( 5 1 8 ) r 5 1 9 1 2 7 计”n险p 1 2 f | g_-i k+ 1，j ny一 扩西 ，iul m rl 勰 为式形阵e*的 口壬 t方阵颤则 机翼跨音速颤振的频域计算 k k = 。2 k 。= 女f k = k p k ，2 k ，k k ( ， ( 5 2 0 ) 则模趔和实物具有相同的颤扳方程。 山( 5 2 0 ) 式必须有。= 1 。同时，如果计算的模型和实物具有相同的m ，p ， 减缩频率b e o 矿，相似的空气动力外形和振动形态，还要求计算得到的速度即为真实 的速度，有 k = 女。2 k ，。= k h k ，k m k ( ，( 5 2 1 ) 在实际计算过程中，若把算例的所有具有长度量纲的参数以及模态( g 口为振动形态) 都缩小工倍，则有 = 1 l k ，= 1 l k “= 1 三2 k ( 。= = 1 l 由( 5 2 1 ) 可得到 k k = k2 k 。，= 1 l 3 考虑到减缩频率相同，k b 。k 。= l ，月a ，= l ，则有。= l ，故 k 。= 三 女。= 1 l 5 在计算中如果用已知的模态在计算气动力，有时会出现发散的情况 要把算例的模态缩小。若仅把算例的模态缩小上倍，则有 轧= 1 三 由( 52 1 ) 可得到 k k = k f k m = 1 l 考虑到减缩频率相同，k b k 。k ，= 1 ，且 。= 1 ，k 。= 1 则有。= 1 ，故 ( 5 2 2 ) ( 5 2 3 ) ( 5 2 4 ) 因此，需 f 5 2 5 1 南京航卒航天人学硕十学位论文 k 。，= 1 五 ( 5 2 6 ) 由此可知，在实际计算过程中，若把算例的所有具有长度量纲的参数以及模念( 即为 振动形态) ，缩小l 倍，为了确保计算的速度与真实的速度相同，则应把m 缩小l 5 ， k 缩小l 3 ，g o 扩大三倍；若仅把算例的模态缩小三倍，则只需把懈口翩自小r 倍。 还有另一种量纲情况：认为振动形态厂有长度的量纲，米，则广义坐标g 无量纲， j 、义质量m 的量纲为，千克。米z ，( 1 p 一d f _ ，所以( ，。无量纲，那么4 ：厂。a s 。u p ， 4 的量纲为米3 ，其余参数的量纲同上。 5 5 求解颤振方程 5 5 1 概述 不论何种方法求解颤振方程，除形成m 及k 矩阵外，均需先讨算某马赫数下的若 1 二个减缩频率k 值的广义气动力。为此应按第二章和第三章所述方法进行气动力计 算，求得广义气动力。在求气动力前，还要进行插值，以将结构节点振型转化为气动 点( 各网格力作用点) 处振型( 参见第四章) 。 在以上工作的基础上，再用v g 法求解颤振运动方程。 5 5 2 广义质量及广义刚度的求法 广义质量可取自计算或实验结果。若取自计算结果，则因自振型的正交性，m 矩 阵为一对角矩阵。若取白共振实验结果，则般m 是满矩阵，即m 。0 ( f ，) 。 但不论广义质量矩阵是否对角矩阵，均设广义刚度矩阵为对角矩阵，且 e 乏兰f j；i=jij ； z ， 世。= 0( ) 其中，g o ，为第i 阶自振频率。 计算广义质量矩阵时，m ，代表了第f 个自振型的广义质量且被定义为： 机翼跨音速颤振的频域汁算 m ，= f f m ( w ) 2 d s ( 5 2 8 ) 其中呱x ，y ) 为在x ，y 点处当地单位面积质量，为第i 个自振型。 广义质量耦合项具有正交性，则可被表述为： m ，= j m ( ，) ，乃嬲= 0 ( ，- ，) ( 5 2 9 ) ( 5 2 8 ) 式可用下面( a 1 ) 和( a 2 ) 定义的无量纲量表示为： 6 0 f 砌( 和甜够 ( 5 r30 ) i = 三( a 1 ) 歹2 孚 ( a 2 ) 其中，0 表示弦向角度的变化值，6 为当地半翼弦长度，b 。为半翼根弦长度，为 半展长。运用积分系数和，。( 合适的积分系数冬使得歹：只，合适的积分系数， 使得口= 配。c ，s 分别代表弦向，展向的坐标) 则( 5 3 0 ) 式呵写成如下的求和形式： m = 7 莓莩厶鲁。m ( 配。) ：2 ( 眈，只) s i n 色 ( 5 3 1 ) 写成矩阵相乘的形式如下： m i = ，一生b ob ( 吼，y , , ) i cs i n 艮血2 ( 吼，歹。) ( 5 3 2 ) 列矩阵协2 ，可从相应的方阵 的列元素得到。惯性项吖。被扛：( 眈，只) 代替， 交叉相乘为 k f ) 时。 当弦向积分系数，( ，和展向积分系数，。的代数齐次多项式已知，机翼质量分布 m ( o ，歹，) 和自振型，已知时，就可计算出广义质量矩阵 m ，从而进一步求得广义刚 度矩阵k 。 5 5 3 弹性方程的频域解法i ， 3 0 南京航空航天人学硕t 学位论文 5 5 3 1 v - g 法( 参数法) 这是目前颤振计算中最常用的方法。 用v g 法求解颤振运动方程时，设系统作简谐振动，则 抽 = q 。扣( 5 3 3 ) 其中国为振动频率，且引入遐想的结构阻尼系数g 。实验证明，振动时结构阻尼 力的大小与振动位移量成正比，因而也和弹性恢复力r 构：k 4 、成正比，而与频率尢关。 此外，它应和速度反向，即和位移有9 0 度的相位羞。则( 5 16 ) 式可写作： m 口) + k o + 辔) g ) = 三p y2 q g ) 将( 5 3 3 ) 带入上式，并做一些变换得： 咖m ) 叫阻卜譬【q 】 g b 或 丑 g = ( 阻】+ 豢蚶阻? i g ( 5 3 4 ) 其中，i 为减缩频率，：堕；b o 为参考长度。 v 彳= 篙1 g = 岳1g c 吲； 以= = 一f 1 一f f l 1 + ! 。 则： 2 青一等 f 2 霄，一青 解( 5 3 4 ) 式，得 a u2a 胁十i 2 j 为振型数，i 为减缩频率数。在求解( 5 3 4 ) 式前，必须先给定p 并没定组k 值。 对每一个以，可求得： 驴一鬻 c s s 鄹 塑塑堕童鲨墅堡塑垫蔓盐笪一 曲“= a 蜘- 2 g ( j 3 6 ) = ( 53 7 ) 由( 5 3 5 ) 式给出的g 值的物理意义是：假象这个机翼的结构阻尼系数等于这个g 值时，则在上述p ，v 条件n 机翼作简谐振动。若机翼的真实结构阻尼小于此值， 表示要再给机翼加一些阻尼，它才能作简谐振动，故机翼实际作发散振动。若g 为负 值，表示要加上一负阻尼即激振力，机翼才作简谐振动。实际上，机翼末受此激振力， 故机翼的振动是衰减的。 对于给定结构，计算颤振速度时，要给定p 值，并设组k ( 或1 k ) 值，求得 相应的v ，g ，0 9 ，通常绘制成g 矿及g o v 曲线，而以g 等于实际结构阻尼值时所对 应的矿值作为颤振速度。 实际上，结构阻尼系数比较难以测定，而设g = o 可以得。保守解。故通常取g = 0 所对应的v 值作为颤振速度，记坐。，相应的可确定0 9 昧：。也可以认为结构阻尼 不为0 ，例如认为9 2 0 0 2 ，这时，可在v g 图上求得相应的_ 瑚。及忙。 在v g 图上求得多个颤振点时，其中颤振速度最小的即为该m 数及空气密度p 情况下的颤振速度。 程序也可自行确定颤振点，其步骤是： 1 在求得的n k k 2 旧_ 解中，依次找出g ，。 o 的点。 2 按线形关系预测g = o 的= 七+ 值。 3 按第二和第三章所述的方法计算k 时的广义气动力，求得新的颤振解。对 新求得的g ：，若 g 。 占 则认为已找到满足要求的颤振点。其中动一小值，例如可取s = 1 0 1 0 一。并 求得相应的础及v + ，它们即是满足精度要求的颤振解。 4 若不满足g ： s 时， 3 2 南京航空航天人学硕十学位论文 若g ： o ，且f g ：j o ，且g ： g h 。，则由，g ：，预测g = o 的新的颤振点。 以下按2 - 4 的步骤重复，直至找到满足要求的解。 5 若迭代次数超过规定的最大迭代次数或不满足4 中的要求时，则迭代结束 取最小的g 值的解为近似的颤振解。 5 5 3 2 p - k 法 其中 在p - k 法中，设 幻) = q 。) g “( 5 3 8 ) p 3o j ( r z j( 5 3 9 ) ，为瞬时衰减率系数。 以( 5 3 8 ) 代入( 5 1 6 ) 式，得 p 2 眦 f g ) + 世) - p v 2 2 q r + i q ，】 g )( 5 4 0 ) 若记 i q ， g w 旦 q ， 办 则( 5 4 0 ) 式可写作 ( p 2 卜p ( p b v 2 ) q ， k + (