（流体力学专业论文）球锥动态稳定性导数数值研究.pdf

国防科学技术大学研究生院学位论文 摘要 本文以非定常n a v i e r s t o k e s 方程作为球锥体外形绕流的控制方程，并采用 n n d 格式数值求解。在非定常计算中应用四阶r u n g e k u t t a 方法和变系数残 值光顺技术来提高计算的时间精度和效率。在e t k i n 理论前提下，本文研究了 超声速流动下动态稳定性参数的计算方法，并计算给出了球锥外形的俯仰阻尼 导数和滚转阻尼导数。 ( 在验证求解阻尼导数的模型及上述求解方法的基础上，本文还对以下动导 数求解的影响因素作j ，初步的研究： ( 1 ) 动网格的生成方法。 ( 2 ) 非定常计算时，壁面压力边界条件的提法。 ( 3 ) 振心位置对俯仰阻尼特性的影响。 ( 4 ) 振幅和减缩频率对阻尼导数的影响。) 本文计算与有关文献实验结果的比较表明，文中所发展的非定常流场的计 算方法及动导数计算结果是可信的，并可以推广到求解其它各类型的阻尼导数 i 。 关键词：非定常流动动导数俯仰阻尼导数滚转阻尼导数数值模拟 第l 页 国防科学技术大学研究生院学位论文 a b s tr a c t t h e u n s t e a d y n a v i e r s t o k e se q u a t i o n sa r et a k e na st h eg o v e r n i n ge q u a t i o n so f t h ef l o wo v e rt h es p h e r e - c o n eb o d i e s ，w h i c ha r es o l v e db yt h en n ds c h e m ei nt h i s p a p e r i nt h e u n s t e a d yc o m p u t a t i o n ，t h ef o u r s t a g er u n g e k u t t as c h e m e a n di m p l i c i t v a r i e dc o e f f i c i e n tr e s i d u a ls m o o t h i n ga r ea p p l i e dt oi m p r o v et h et i m ep r e c i s i o na n d i n c r e a s et h ec o m p u t a t i o n a le f f i c i e n c ya c c o r d i n gt ot h ee t k i n st h e o r y ，t h ed a m p i n g i n - p i t c ha n d t h ed a m p i n g - i n r o l ld e r i v a t i v e sa l en u m e r i c a l l ys t u d i e ds o m ef a c t o r s t h a ti n f l u e n c et h ec a l c u l a t i o nr e s u l t so f t h e s ed a m p i n gd e r i v a t i v e sa r es t u d i e di nt h e p a p e r a sf o l l o w s ： ( 1 ) t h em e t h o d s o f c r e a t et h em o v i n gm e s h ( 2 ) t h eb o u n d a r yc o n d i t i o no fp r e s s u r eo n t h ew a l l ，a tu n s t e a d yc o n d i t i o n ( 3 ) t h ei n f l u e n c eo f t h er e d u c e df r e q u e n c ya n d t h ea m p l i t u d eo f o s c i l l a t i o n ( 4 ) t h ei n f l u e n c eo f t h ec e n t e ro f g r a v i t y t h er e s u l t sb e i n gi d e n t i c a lw e l lw i t ht h ee x p e r i m e n ti n d i c a t et h a tt h em e t h o d s d e v e l o p e di nt h i sp a p e r t os o l v et h ed a m p i n g - - i n - p i t c hd e r i v a t i v e sa n dt h ed a m p i n g - - i n r o l ld e r i v a t i v e sa r es u c c e s s f u lt h er e s u l t sp r e s e n t e di nt h ep a p e ra r ec r e d i b l ea n d t h em e t h o d sc a r lb ee x t e n d e dt oc a l c u l a t eo t h e rd a m p i n gd e r i v a t i v e s k e yw o r d s ：u n s t e a d yf l o w d a m p i n g i n - p i t c h n u m e r j c a ls i m u l a t i o n d y n a m i cs t a b ji t yd e r jv a t iv e s d a m p in g in - r o i l 第1 1 页 国防科学技术大学研究生院学位论文 符号说明 当地库朗数 受格式稳定性限制的临界库朗数 轴向力系数 阻力系数 滚转力距系数 滚转阻尼导数 绕质心的俯仰力距系数：c 。：_ 生 兰p 0 s l 。f 俯仰阻尼导数 法向力系数 压力系数：c ，= 华 ：p * p 。 升力系数 单位体积总能 笛卡尔系下的无粘通量 贴体系下的无粘通量 笛卡尔系下的粘性通量 贴体系下的粘性通量 j a c o b i 行列式值 热传导系数 参考长度 来流马赫数 压力 普朗特数 雷诺数 参考面积 温度 周期 迭代推进时间 笛卡尔坐标系下的求解矢量 第1 i i 页 口 v v g | ( u g g 峨 九凡乐乐 一上 丫 一 o 、 、 o 哑巴o g 勺巳 q e c e正反k乐，。k帆p竹k s i ， 国防科学技术大学研究生院学位论文 下标 w 上标 贴体系下的求解矢量 笛卡尔系下的速度分量 笛卡尔坐标系坐标 攻角 侧滑角 比热比 光顺参数 贴体坐标系坐标。亭：流向；卵：周向；f ：径向 俯仰角 粘性系数 密度 滚转角 来流参数 壁面参数 对时间求导 有量纲量 第页 w ： f 、 、 、 v y 叩 、 、 、 一u x 口 y 占善口 p y 国防科学技术大学研究生院学位论文 第一章引言 随着目前我国载人飞船、小型化弹头以及高机动飞行器的研制及发展，使 得如何进一步准确确定动导数的问题更为重要。大攻角和带侧滑情况下的机动 飞行出现的非线性流动现象会对飞行器的动稳定性产生严重影响，要预测在此 条件下的飞行器的动态特性和设计其飞行控制系统，以及对飞机进行失速尾 旋分析，就必须获得其动态稳定性参数。航天再入体在再入过程中，飞行速度 跨越了高超、超声速及亚跨声速区域，飞行器动态特性呈强烈的非线性变化， 并且由于烧蚀和侵蚀的影响，使得气动外形不对称而产生不对称气动力，从而 可能使再入体发生滚转共振、滚速过零和滚转异常，电可能由于配平攻角过大 而产生过大的横向过载，带来结构强度方面的问题。因此，动导数数据不仅与 再入安全直接相关，而且也是控制系统设计必不可少的气动数据之一。 在决定动态稳定性的诸多参数中，俯仰阻尼导数是决定飞行器攻角振幅衰 减的重要参数，而滚转阻尼导数在飞机、导弹、航天再入器的研制中所起的作 用亦十分重要。如现代飞行器在做大攻角机动飞行时所遇到的机翼摇滚现象就 与滚转阻尼特性的变化密切相关，而小型化弹头运动过程中预测误差的大小及 滚转共振、滚速过零等现象的研究均与滚转阻尼导数相关。 国外早在六七十年代对载人飞船的研制过程中，就对飞行器动态特性的问 题高度重视，如美国在“阿波罗”载人飞船的研制过程中，单是动稳定性风洞 实验就用j - 9 套模型在1 4 座风洞中进行了亚、跨、超和高超声速实验共约7 0 0 多小时。在“双子星座”飞船的研制中，仅俯仰阻尼导数和偏航阻尼导数就在 各种攻角、侧滑角、马赫数下进行了大量实验。从这些可以看出动导数研究对 型号设计的重要性。 目前，国内外在实验和计算方面都取得了很大进展，但同时各自还存在着 不少问题尚待解决。 在实验方面，确定动导数的主要方法有强迫振动法、有限自由振动法、自 由翻滚法和模型自由飞等方法。强迫振动法技术设施较简单，但振幅较小，振 心位置难以与实际重心相一致，并存在尾支架干扰；有限自由振动法振幅较大， 但实验只允许在配平稳定区，并存在支架干扰；自由翻滚法的攻角范围不受限 制，振心位置与实际重心位置一致，但有横杆支承干扰：模型自由飞实验没有 支架干扰，但实验的振动周期较少，纵向和横侧向运动的相互干扰不易分清等。 另外，在滚转阻尼的实验中，由于其量级往往很小( 如短粗体) ，对实验装置 和测试手段要求很高( 如机械阻尼、支架干扰等) ，因此实验测试难度较大。 第1 页 国防科学技术大学研究生院学位论文 此外，一些新兴的方法，如磁悬浮风洞技术仍处在起步阶段。由于以上各种原 因，实验研究动导数还有较大误差，据有关文献分析，实验确定动导数其分散 度在2 0 左右。 9 0 年代以前，国内外在计算方面主要集中在一些经验和半经验方法上。如 修正牛顿理论、内伏牛顿流理论、牛顿一玻尔兹曼理论、修正激波一膨胀波理 论等。由于这些理论未计及气流分离、再附和尾流等效应，对于复杂外形大攻 角流动难以适应。 数值求解动导数的文献在9 0 年代以前十分少见，9 0 年代以后，有h u i 在 摄动理论基础上发展的谐振摄动法，该理论在小攻角范围内对细长体外形获 得了成功，但对短粗体外形或攻角较大的流动，由于小扰动理论的局限性，结 果难以令人满意。在直接采用非定常e u l e r 方程或n a v i e r s t o k e s 方程求解扰动 流场方面，比较典型的文献有任玉新、刘秋生、沈孟育求解e u l e r 方程计算 小攻角球锥外形的俯仰阻尼导数。刘伟、张鲁民利用“双时间法”求解薄 层近似的n a v i e r s t o k e s 方程模拟有粘扰动流场，给出了“类联盟号”飞船钝体 外形的俯仰阻尼导数。此外，w e i n a c h tp 等主要通过模拟锥运动来确定俯仰阻 尼导数，但攻角一般都不大i ”1 1j 。 在采用数值方法确定动导数的种类上，目前所见到的主要以俯仰特性为主， 其它阻尼导数尚不常见。典型的有刘伟【8 】在高超声速下给出了短粗体偏航阻尼 导数的计算方法。而滚转阻尼导数计算文献十分少见，这一方面是由于动导数 涉及非定常流动，理论上还不十分成熟；另一方面则是由于滚转阻尼导数量级 往往很小，并涉及表面摩阻等技术问题，计算难度很大。 本文从e t k i n 理论【9 】出发，即：当飞行器的状态变量在时刻t 以前的无穷长 时间内( 一o o ，墨t ) 是时间的解析函数时，气动力是状态变量及其对时间的 各阶导数的函数。在数值模拟中采用n a v i e r s t o k e s 方程描述非定常流场，应用 n n d 格式【”】离散方程，并采用四阶r u n g e k u t t a 方法1 及变系数残值光顺技术【1 2 1 提高时间精度及计算效率。在充分验证方法的基础上计算了不同攻角下钝锥俯 仰阻尼导数及滚转阻尼导数，与有关实验 1 3 儿“1 进行了比较，结果令人满意。 本文还研究了动网格生成技术、非定常壁面压力边界条件提法、振动中心及振 幅对动导数的影响等。 此外，本文还结合文献 1 5 】的理论分析，验证了在目前网格条件下，n a v i e r s t o k e s 方程与薄层近似的t l n s 方程在气动力的计算方面基本无差别。据此， 在本文的计算中均采用t l n s 方程。 第2 页 国防科学技术大学研究生院学位论文 第二章控制方程 在本文计算中，采用了贴体坐标系下的n a v i e r s t o k e s 方程、薄层近似的 n a v i e r s t o k e s 方程( t l n s ) 。它们都是由笛卡尔坐标系下的n a v i e r s t o k e s 方程 变换，简化得来，下面对这三个控制方程进行介绍。 2 1 笛卡尔坐标系下的控制方程 坐标系选取如图1 ，取为右手系，与一般飞行动力学中机体坐标系略有不 同，即x 轴指向与来流方向一致。定义如下无量纲量 g ，y ，z ) = g ，歹，享) p = 汤。 t = f | f 。 t = r e o z 0 ，u w ) = g ，铲，访) 吃 p = 芦瓦吃： u = p | h 。 式中：带”者为有量纲量，p 。，瓦，疋，瓦为目由来流参数。三为物 体特征长度，在本文中球头、球双锥、球锥f 均取为头部半径。对于忽略体积 力的常比热的气体，三维非定常n s 方程可写成下列无量纲形式 一o u + 丝+ 堡+ 箜：坠+ 丝+ 堑 ( 2 1 ) 动a ) ： 砂 a z瓠 砂 e z 式中 u = 1 e 。面1 p p u p v p w e e = p u 伽。+ p 珊v p u w ，+ p ) u 0 f 搿 r 廿 r 世 ，罢 ，f = l p u v p v + p 口w ( e ，+ p ) v r e g = , o w p u w p v w o w 。+ p ( e ，+ p ) w 0 r f t 母 r ，| 】 考 第3 页 国防科学技术大学研究生院学位论文 g 1 - 上 r e 印告+ j 1 p , 2 _ _ ! ：2 - f i v 2 ) 七2 ( y - 1 儿坐 p r r 。： f 坐+ 竺+ 竺1 + 2 坐 k 瑚【瓦+ 丽+ ij + 2 p 瓦 k ：z f 竺+ 尘+ 一o w1 2 竺 以l 瓦+ 万+ 西j 十2 p 可 l ：a f 塑+ 堡+ 型1 + 2 丝 k “l 瓦+ 面+ 西j + 2 p 瓦 2 川、号+ 刮 叫【瓦+ 瓦j 铲程+ 割 k 2 l 石+ 夏 fo v功1 k 2 l 瓦+ 万j 、 方程中的无量纲参数 r e ：盈型 t 。 p r ：星空 k 。 五满足s t o k e s 假设 “由s u t h e r l a n d 公式 m 。：垒 口。 y = 14 由完全气体假设可得无量纲的气体状态方程 音速方程可由状态方程推出 c = 1 1 0 4 k | t 。 p ：l p 丁 1 抄，。 口= 4 b 第4 页 1il塑品 o 。k l 栅 妒 rv+ 荔 r” ，1，。l 、i， c c o + 一十 = 一， f，l 功：p 一 一一 姒 国防科学技术大学研究生院学位论文 2 2 贴体坐标系下的控制方程 为了适用于一般物形和计算，引入岵体坐标系g ，仇f ，f ) r = t j 喜= 善；：，j ，z ；t ( 22 ) 1 叩= 叩g ，y ，z ；t ) p 卅 k = f g ；y ，z ；t ) 则n s 方程变换为 掣+ 豢+ 要+ 等：鲁+ 誓+ 拿 a z 8 8 研a ga a q0、。 式中：驴为求解矢量，e 一、f 一、g 一为无粘通量，e 一、瓦、巨为粘性通量。 订：。p l e = j 1 蟾t u + 正十；+ 曲= j f = j 一1 0 ，u + 叩z + 叩，+ ，7 二g ) = t ， 弓= t ，一1 9 ，u + f 五+ f ，+ f f ) = l ， 丘= ，。晤乒，+ 善只+ 善，g ，) f = ，。白乒、+ 7 ，e + 7 g ) 己，j = j 。、乜e ，+ f ，+ o 、 删 舀+ 善。p 印+ # 。p p w u 十吉：p ( e ，+ p 如一鲁p 口i ) u v + q 。p p 1 + 7 7 、p , o w v + 7 7 ：p ( e ，+ p p r b p p w 瓣七。p 印w + 。p 谛+ f ：p 暇+ p 弦一4 , p 第5 页 国防科学技术大学研究生院学位论文 五= 鼻+ 善。”+ 毒。v + 善：， 、1 = 7 7 f + 7 7 u + 7 7 r v + 7 7 ：w 谛= f ，+ f ，“+ f 、1 ，+ f = 各雅可比度量系数为 i ；= j b ，z ；一z 。y ； 告、= ，( ：。x i x 。z ； l 善：= ，0 。y 一o ；= j b ；：。一z ：y ， f ，= ，b k z 。 【f ：= l ，g ；一坎b h = t ，坼。一o y ) 7 7 、= l ，( ：x 。一x ：。) i ，7 ：= ，k 儿一儿_ ) f 喜，= 一b ，喜。+ y ，舌、+ ：。掌：) 玑= 一0 。7 7 。+ y 。町、+ ：，7 7 ：) i f ，= 一0 1 f 。+ y 。f 、+ ：。f ：) r 只是，的函数，与r 、y 、：无关：r ，= r 、= r = = 0 2 3 薄层近似方程( t l n s ) 在来流雷诺数很高的情况下，粘性主要在靠近物面的边界层内起作用。由 i 于边界层厚度占o c r e2 ，当r e 1 0 5 时，边界层很薄。本文计算算例均为大雷 诺数情况，并且由于计算机速度和内存的限制，三维运算时流向、周向所取的 网格数较少，只在物面法向进行了指数加密，因此物面法向的粘性作用被计入， 而流向和周向的粘性作用影响不大。在巨中忽略昙、善导数后形式如下： d go r 其中 j ：，一卫 r e o m 一：+ m2 x | 3 m ，+ 聊。f 3 ”l ”c + m 2 ：3 ；7 1 1 ，”；+ m ：曾，“+ f 。、，+ f ：w ) 3 一，1 = 幺2 + f 。：+ f ：2 ”7 ：= f ， + f ，+ f ：w ( 24 ) 第6 页 端 矧一 国防科学技术大学研究生院学位论文 m ，2 2q _ v 2 + 1 4 , 2 x + 矿而1 简化后的n s 方程可称为t l n s 方程，可写成 面葩舔丽舔 一+ + + = 一 6 z a 8 qa a f 25 ) 在大雷诺数且没有大的分离存在的情况下，t l n s 方程可以得出很好结果。 第7 页 国防科学技术大学研究生院学位论文 第三章数值方法 对于三维非线性的n a v i e r s t o k e s 方程，本文采用差分万法求解。对定常问 题，采用隆式n n d 格式计算；对非定常问题? 本文采用四步r u n g e k u t t a 方 法构造的显式n n d 格式，并引入变系数残值光顺技术进行求解。 31 差分方法 3 】差：兮格式 采用有限差分求解控制万程 对于一维气动万程组 丝+ 笪：o 。ia x 半离散化的n n d 格式为 差分离散采用张涵信院士提出的n n d 格式。 ( 甜= 一如i 矗 其中 f + ! 。 o 。 上述空间离散的n n d 格式也可写成如下形式，见文献 1 0 。 t + ；2 ；+ f + z ) 一i 1i 爿 ，+ ：，+ ； + 三爿：；m i nm 。a ( 呜u ，+ ；u 一1 ：a _ - , m i nm o d ：曩c ， ( 31 ) ( 32 ) 第8 页 ” q 0 0 0 0 国防科学技术大学研究生院学位论文 a = a + + a - = 荔为f 的j a c o b i a n 矩阵，矩阵。+ ，月分别为爿的正负特 征值分裂矩阵。 m + ；2 彳：；一爿j ；4 ：；为爿! 在，和沁二2 两点上值的某种平均，如算术 平均、r o c 平均等。 对于一般曲线坐标系下的三维n s 方程： 警+ 善+ 筹+ 毒= 专- 筹+ 等 o ， 8 t e 8 q0 0 e na 、j 差分离散形式可写成 。= 佤昏瞩+ f i ；g , - c 5 , e - c s , f - t 万) ( 3 s ) 万为中心差分算子，万即为上述n n d 算子。 当采用f t 显式离散时间项时 ( ，。i 卜u 。? t + r h s ，( 39 ) ，的选取：为简便，不妨取毒= a t = f = l 。记j a c o b i a n 矩阵 4 = 嘉b = 堡c g uc = 嘉 ( 31 0 ) 文7e u 。 其特征相似主对角矩阵a 。、a 。、a 。写成如下统一形式 人。= 西昭似，拦，硝，拦，霉)( 31 1 ) 特征值( = 1 , 2 ，3 4 ，5 ) 为 前= 葛= 鹚= 女，+ k x u + 女一，+ k w = 0 麓= 0 + 口女：2 + 女，2 + ： 雉= 0 一口以2 + 女! + 女= 2( 31 2 ) f 喜对a = l 玎对b i f对c 口为声速。在不考虑粘性影响的下，格式稳定时的时间步长为 ，t 删忐 第9 页 国防科学技术大学研究生院学位论文 l 力i 。= m a x l 笛i ，l 破l f = m i n u e t ，垃，a t f ) f 3 1 3 ) 对于一维流动，n n d 格式的稳定库朗数c n = 妄。 312 隐式n n d 如果采用显式格式求解n s 方程，花费时间非常大，因此本文在定常流计 算中应用经对角化处理，计算效率较高的隐式n n d 格式。该格式在空间上除 个别点外，均具有二阶精度，离散形式如下 f ，地一嘉“。嘉+ 毒嘉 = r h s 血b 4 ， 式中p o u s 即为上述n n d 格式的差分离散，引入因式分解并对左端标量化处理 得到 愿( + f 毒天；k ( ，+ a t g 。天。) r i l 心( ，+ 出毒天f k l 驴”1 = 心s - f 可”1 ：可“+ a u ”1 ( 31 5 ) 兀。= a a g 瞬，霉，君，冠，君) 为雅可比矩阵之特征值在i ：i + l 点上的某种平 均。这里占，、最、曼采用一阶迎风处理。再对上式进行a d l 分裂后得到a 型 的隐式n n d 格式 f ，+ f 占，天。) 盯+ ：r ：1 r t t s a t ( ，+ a t , ；f i 。) 可“= 巧1 r 可 ( ，+ f 毒天) a u 。+ = 月；1 r 。驴” 驴”1 = 趣可+ 驴”1 = 驴”+ 盯 ( 31 6 ) a d i 型的隐式n n d 格式可直接进行标量求解，避免了矩阵求逆。数值实验表 明，格式c n 数可取高达4 0 ，而每一步的计算量只比显式格式多3 0 ，因此 极大的提高了计算效率。 3 2 四步r u n g e k u t t a 方法 3 2 】四步r u n g e k u t t a 方法形式 第1 0 页 里堕塾兰垫查奎兰丝壅竺堕兰堡丝圣 上述隐式n n d 格式由于对角化以及线化处理，在时间上只有一阶粘度。 对非定常流动，为了提高时间精度，本立采用如下四步r u n g e k u t t a 方法 e ) _ 玩” 扩”= u j ”一a l a t r t l s 。o 可，扪= 可0 、一口，a t r i q s ( 1 玩。j = 可，。一嵋删帮 可，“= 盯一a a t r h s o ) 猡”1 = 可4 c1 7 1 ( 口，口：，a ，口。) = ( j 1 ，吉，- ) 该格式具有二阶时问精度1 ，r h s ( ”) 为第n 步推进时的n n d 格式离散项。， 为时间步长，在非定常流场计算中，a t 全场相同。对本文计算的谐振流场， a t 由周期中的推进步数决定 2 丌 k 百一 7 3 1 8 ) “ c l c k k a ， n 。为同一周期内的推进步数，为减缩频率。 3 22 四步r u n g e k u t t a 方法的稳定陧 格式的稳定性限制了非定常计算中a t 的选取 方法的稳定性。 对一维模型方程 塑+ 望：o f 是“的函数，且，= 口及臼：篓，口为其特征值。 半离散的n n d 格式可写成 = 去哎0 一。甜1 m i n m o j k a f ；蟛 + 去m i n m o a 【够：；_ a f * i 7 l l i i 一一j 下面分析四步r u n g k k u t t a ( 31 9 ) + 蛐n 2 , 5 x 一m l n 2 工 由m i n m o d ( x , y ) 的定义，上式写成 ( 黔一去哆j 1 一去t 够乏+ 去t ：，：； + 去t 够萼一击t 。够：； 这里0 k ，1 ，= 1 ，2 ，3 ，4 ( 3 2 0 ) 式也可写为 式中 瞎) = 一去哆：；一去哆：! 由( 32 1 ) 可知 土口s 三 22 ( 32 1 ) 式也可写为 瞎) 一以”量恤：； 舀= 去2 化_ 归4缸l ”厂 肛去( ，如咄归 由于口y 和口j “中必有个为0 ，所以只需分析特征值为正时的情况。 当口? o 时，a j “= 口? ，口i ”= 0 ，故夕= 0 ，( 3 2 3 ) 式成为 ( 32 0 ) ( 32 1 ) f 32 2 1 ( 3 2 3 ) 1 2第页 、，l、i ”、) 一，- 目，一2 鲈。 v ， n 一2 h 一： 酽， v ， r 州l r 叫l 啪 啪 川j1叫 l 4 ，阿，k 一 一 ，- 3 忙 ，2一2 一 一 一一rl i l i l 口 国防科学技术大学研究生院学位论文 ( 豺一“鼍 令： 尸0 j 。) = 抛“= & 0 7 一”：) 应用傅立叶分析，令? = a e “p 州“) 有 6 f ? 圳， ( 32 4 ) ( 3 2 5 ) c o s ( 纽) ) f + i as i n ( f l a x ) a t( 32 6 ) 9 叩 此处口= = 生称为波数。 五 对于一维标量方程，四步r u n g e - k u n a 方法如下 甜( o ) = “一 ”( t ) ：”( o ) 一! 炉“( o ) 1 4 、， ) ：扩) 一；f pc ，( 】 ) j - 1 委护0 ) ”( 4 ) = 。j a t p ( 3 ) 1 h n + 1 = ( 4 ) 由傅立叶放大系数g ：4 n _ + l ，并注意到( 32 6 ) 式，可得到 ” g = g o ) = ，一：+ 了z 2 一詈+ 丢 式中 ：= , 4 1 - c 0 4 0 ) ) + ，口s i n ( o r ) 稳定性要求 同理，对口? 综合以上两式 口，：坐f x i o = 8 缸 一l ( k ：- k i 归 口f - 万，】 法分析，可以得到最大稳定库朗数为 第1 3 页 ( 3 2 7 ) ( 32 8 ) 柳 珊 p 0 淼 用 = 到 采 m 得卜上以壮r 可陆他恤吼 o m 4纛 、k 一 况 国防科学技术大学研究生院学位论文 3 3 变系数残值光顺技术 显式四步r u n g e k u t t a 方法的时间步长受稳定库朗数的限制，因此要提高 非定常流场的计算效率，应对格式进行一定的处理，本文采用j a m e s o n 提出后 经j o r g e n s o n 和c m m a 改进的隐式变系数残值光顺方法。在r u n g e k u t t a 方 法每推进一步后对残值进行光顺处理 ( 1 一嚷x 1 一蛾。x 1 一哦) 盯= 可 ( 33 1 ) 可为每推进一步后的残值，p - 为光顺后的残值。s 为光顺参数，由当地库 朗数决定 比舭b a x 陪h 。s z ， c f l 。t 。tj 。a t c f l * 为受原格式稳定性限制的库朗数，c f l 。为当地计算库朗数。和常系数 光顺相比，在远离物面处，由于网格间距大，当地库朗数c f l 。，。远远小于 c f l * ，因此s 为零，这样选取光顺参数避免了远离物面处的物理信息的损失。 对( 33 1 ) 式进行三对角分解得 ( 1 一s 占。l 盯6 = 驴 ( 1 一峨。) 可“= 驴6 ( 1 一哝必可“= 可” 可+ = 驴6 “ u ”1 = u ”+ u + ( 33 3 ) 数值实验表明，由于残值光顺可分解为三对角标量追赶完成，计算量只增加了 1 5 ，而c f l 。却提高了1 0 倍以上，因此，残值光顺大大提高了非定常流场 的计算效率，与一般隐式格式相当。此外，残值光顺方法还保留了显式格式的 简明性，有利于复杂边界条件处理。经j a m e s o n “1 证明，隐式残值光顺在理论 上是恒稳的。 j o r g e n s o n 证明残值光顺对流场的影响为o ( 心3 ) 量级m 】。因此，只要s 不 是太大，隐式残值光顺对流场的计算精度没有影响。并且隐式残值光顺是在空 间上进行，因此它对时间精度没有任何影响。 第1 4 页 国防科学技术大学研究生院学位论文 第四章计算网格 本文在计算中采用贴体坐标系，坐标系选取如图1 。由于计算外形为球锥 体，采用代数方法生成的网格，其光滑性、均匀性、壁面正交性都能得到较好 的满足。另外，在非定常流场计算中，代数方法能大大节省动网格生成时间。 4 1 定常网格生成 如图1 代表物理区域，图2 代表计算区域。在图2 中，k ：1 代表物体壁 面，k = n 代表外边界。记b ，y ，z ) 和g 。，y 。：知) 两个点分别为同一条f 线上 两个面上的对应点坐标，i 为g ，y ，：1 ) 和g 。，y ，知) 两点连线上的单位向量， 第k 层坐标面上相应点的笛卡尔坐标可表示为 x 。，y 。，气) = g ，y ，& ) + l ( h ，y 。，知) 一x 。，y ，2 。】- 口g 户 ( 41 ) i 表示取模，a ( k ) 为壁面加密函数 州小h 管 ”阶( 甜 ：， 式中b = i * f 坐标线上网格点的分布由参数控制，为大于1 的量。当 口越接近于1 时，网格点在壁面附近加密越厉害。 4 2 动网格生成 本文所进行的非定常计算是球锥体绕质心作微幅简谐振动，非定常计算在 惯性系下进行，以球锥振动前定常状态为基准运动，并且以此时的机体坐标系 为惯性坐标系。因此，必须在不同的时刻重新生成网格。下面以俯仰振动为例 来阐述一下本文所采用的三种动网格生成方法。 4 2 1 动网格一( n e t 。) 外边界网格固定不变，物面边界随物体做旋转振动，流场内部点由内外边 界插值得到。假定初始物面坐标为b 、，y w o , 2 w o ) ，质心坐标b 。，y g , z 。) ，相对质 心平移后的物面坐标为b 0 ，y k ，z 0 ) b 0 ，y ，z 0j = b 、- x t t ，y 、一y 。，z 、一2 9 ) ( 4 3 ) 绕质心旋转口后物面坐标为b ：0 ，y ：o ，z ：0 ) 第1 5 页 国防科学技术大学研究生院学位论文 转动后的物面坐标在原坐标系中为0 。f f ) ，。( f ) ，：。( ，) ) 内层网格则由式( 41 ) 插值得到。 ( 44 ) ( 45 ) 4 22 动网格二( n e t 2 ) 整个物理平面的网格由一次旋转变换得到，在这种情况下，只需将上节中 物面坐标 ，( f ) y 。( f ) ：。o ) ) 替换为流场所有点坐标k ( ，) 儿( f ) 气0 ) ) 即可。 4 2 3 动网格三( n e t 3 ) 在求解区域内，以适当的等f 网格线为界，内部网格( 包括界线) 由一次 旋转变换得到，外部则由分界线与固定的外边界插值得出。这样得到的动网格 外边界仍然保持不动。插值时要尽量保证网格在分界线两侧连续变化。在旋转 转动幅度较小时，这样得到的动网格仍然比较均匀，畸变不大。 第1 6 页 口 乜。龇 一。一 、jill， 0 毛 ，。11。1。 + 、illli， ”- z y 工 ，。1。，l l i n 川“ o o ok 儿o ，。：，l 国防科学技术大学研究生院学位论文 第五章边界条件与初值条件 由于本文计算涉及到定常、非定常计算，它们的初边值条件是不一样的 并且，非定常计算中俯仰振动和滚转振动也有区别，因此，下面分别说明。 5 1 定常部分 51 1 初值条件 外边界给定为来流值，即 肚元 p 讽= 去 壁面提无滑移条件和等温壁条件 旷= 0 t = l 内层点由它与壁面的距离插值决定 = c k x 。+ ( 1 一c k i x 。 式中x = 0v ”p7 1 ) 7 c 七=i k ，y 。，气) 一g y ，：。) l g 。，y 。，：。) 一g 。，y 。，：。) p = p 。= 1( 51 ) 下标n 代表外边界w 代表物体壁面 k 代表求解区域内部。 51 2 边值条件 ( 52 ) ( 53 ) 如图4 ，外边界按流动划分为入口边界及出口边界，入口处流速为超声速， 特征值均大于零，特征线指向解区，因此给出所有来流条件 0 。v 。w 。p 。p n ) 7 = 0 。1 ，。w 。p 。户。) 7( 5 4 ) 出口处仍为超声速，特z 、 区x 外，因此不需要给解析边界条件，只 需提数值边界条件，本文采用简单的一阶线性外推 - o s ， 对上式进行一阶差分得 x x = x h( 5 6 、 第1 7 页 国防科学技术大学研究生院学位论文 = 恤1 wpp ) 。( 57 ) 对于采用一般坐标系变换而出现的奇性轴( ，= 1 ) ，j a c o b i a n 变换行列式 出现奇性，采用轴平均处理。 在俯仰振动流场计算中，左右流场对称，只需求解半边流场引入对称边 界条件处理对称面，以，一面为例 “o ，。十1 ，k ) - - ”o ，。一1 、七) v ( f ，一+ 1 ，k ) = 一v o ，i ，。一1 ，) w ( f ，一l ，七) = w o ，。一l ，k ) p ( f ，。、+ 1 ，_ i ) = p o ，j 一1 ，k ) p ( j ，j 一+ l ，女) = p ( i ，。一1 ，女) ( 58 ) ，= l 面处理同上。 在滚转振动流场计算中，左右流场不对称，需求解全空间流场。在周向 j = 1 与j = k 为同一子午面，提周期性边界条件 ”o ，。，j i ) = u ( i ，1 ，k ) 、，o ，i ，k ) = ，( ，1 ，k ) w o ，一，k ) = w ( f ，1 ，斤) p ( f ，。，k ) = p ( f ，1 ，k ) p o ，j 。，k ) = p ( f ，l ，女)( 59 ) 壁面条件：旷= 0 ，t = l 。 此外，由边界层近似对壁面压力提数值边界条件：! 翌：0 。 5 2 非定常部分 初值条件，以定常流场计算出的收敛解作为非定常流场的初场。 边界条件：壁面提无滑移条件和等温壁条件，外边界提法同定常情况。 非定常壁面压力数值边界条件的提法很大程度上依赖与数值经验，本文对 文献中常见的三种典型提法在数值计算中进行比较分析。 1 法向动量方程形式如下 一亲= 户曙+ “罢+ v 考+ w 老) 一 等+ 等+ 等惟 锄i l 西良却出l 缸却玉。 d。-lp ( 詈4 - ”塞+ v 堡+ w a y 喜 一( 誓+ 誓+ 誓m i ”+ v + l i 二+ 卫+ 二i 胁 la苏瑟jl 缸却” 第1 8 页 国防科学技术大学研究生院学位论文 + ip 詈+ ”瓦0 w 州万a w 十w 警 一 等+ 等+ 警汁：十ip l 百栅瓦+ 1 万十w 西j l 言+ 亏+ 蔷肝 2壁面加速度修正公式为 宴：一厕亓 式中厅：f 害，孚，掣1 为壁面加速度 a2 l 西西西j 为壁皿加速发 再= 0 ，押，心) 为壁面单位法向矢量。 3 n g n i t g 时，呈：o 嘶 ( 51 0 ) ( 51 1 ) 第1 9 页 国防科学技术大学研究生院学位论文 第六章动导数求解方法 飞行器动态稳定性的研究对于现代飞行器的控制系统设计具有十分重要的 意义，然而动态稳定性的研究十分困难，目前理论和实验都不成熟。在理论上， 主要是由于时间历史效应，飞行器的运动与其所受到的空气动力相互耦合，因 此在严格意义上，要研究飞行器的运动响应及动态响应必须联立求解飞行器 的运动方程以及绕飞行器流体的运动方程。这种方法虽然精确，但难于实现， 目前的计算条件尚难满足。 为了避免耦合运算，须引入描述空气动力与飞行器运动之间关系的数学模 型。在现有理论中，最为完整的属于t o b a k 的非线性指示泛函理论，但不幸的 是t o b a k 的指示泛函在应用上存在种种困难。而应用最为广泛的理论为e t k i n 模型。 e t k i n 假定可表示为：当飞行器的运动状态变量在时刻f 以前的无穷长时间 内( 一m ， ，) 是时间的解析函数时，气动力是状态变量及其对时间各阶导数 的函数。e t k i n 理论考虑了时间历史效应对动导数的影响。对于质心以定常速 度运动的飞行器来说，其状态变量包括：攻角口，侧滑角，绕x 、几z 轴的 旋转角速率p 、q 、r 。因此如果用上表示做非定常运动的飞行器所受作用力或 力矩，l 可表示为 ， l ( t ) = l 虹，应，p ，p ，芦，一，g ，口，辱，i ，- )( 61 ) 一般来说，非定常气动力( 力矩) 还应与时间，相关，但当非定常运动具有周 期解时，在上述定义中不显含时间变量f 。 6 1 俯仰阻尼导数求解方法 给定强迫振动形式 口q ) = 0 = a 。一口。c o s k t( 62 ) 式中a 为瞬时攻角，口为俯仰角，a 。为平均攻角，口。为振动幅度，k 为减缩 频率。k 0 时，= 口( o ) = a 。一a 。，即为基准状态的攻角。 对于飞行器的纵向运动，如果其质心速度不变，则在体轴系中，确立运动 的独立变量只有迎角口和俯仰角速度q ，以俯仰力矩系数为例，由e t k i n 假定， 俯仰力矩系数c 。，可写成 c 。= c 。( 口o l 应( f l 应0 l ，g o l 乎( f l 学( f l )( 63 ) 在满足假定 第2 0 页 国防科学技术大学研究生院学位论文 ( 1 ) 基准飞行状态为对称定直飞行且扰动幅度很小。 ( ”计算中只考虑一阶动导数，忽略高阶动导数。 ( 3 ) 气动力达到周期解。 泰勒展开并略去高阶小量，有 c ，( f ) = c 。+ c 一+ c 。a o ( 64 ) 式中 f c m o ：= ( c c m 。- + k f 2 ( c ) 一。+ ：c ( f 。) + + k c 。( 。c ) = = ? 2 + c ( es ) 【f = ( c 。+ f ) 一七：( f + c ) + 。 当t 不很大，且忽略高阶动导数的影啊时，上式可简化成 c ，。2c + c ( 66 ) c 。即是本文所求的绕质心振动时的俯仰阻尼导数。 当非定常振动过程达到谐振解时，c 。的解可通过对式( 64 ) 的数值辨识求 出，本文采用下面的积分公式 一去r c 觯) s l n 融 ( 6 7 ) f 。为积分起始时刻，t 为振动周期：瓦= t 2 7 9 。利用该式可以避开静态导数 c 。的计算，当非定常计算二至三个周期以后，流场就可以达到比较好的谐振 性。经文献【4 】分析，利用积分式辨识出的动导数与最小二乘法非常接近，但最 小二乘法须确定静态导数c 。 6 , 2 滚转阻尼导数求解方法 当滚转角妒以小振幅振动时 缈= y o y l c o s k t 当，= 0 时，y = y 。一p 。= 0 ，也可将上式写作 y = ( 1 一c o s m ) 是小量，攻角a ，侧滑角卢随滚转角y 的变化而变化。 设来流速度矢量为旷，初始攻角为口 f ”= v c o s l 7 0 1 ，= 0 1 w = vs i n 口o ( 68 ) ( 69 ) ( 61 0 ) 第2 1 页 国防科学技术大学研究生院学位论文 如图3 ，绕x 轴转动y ( f ) 角后，瞬时攻角口( f ) ，侧滑角o ) 可推出为 口= a r c t g t g a oc o s 矿( f ) 】 = a r c s i n s i n a os i n 妒o ) ( 61 1 ) 由泰勒展开可得 a = a o + o 妙2 ) = s i n “o y + o 3 )( 61 2 ) 略去高阶小量有 口= t 2 , o = 常数 声= s i n 口o ( 61 3 ) 因此，对于微幅滚转振动，状态变量只有