（流体力学专业论文）三维带操纵面机翼Euler方程非定常气动力计算.pdf

南京航空航天大学硕士学位论文 摘要 本文研究了三维带操纵面机翼定常和非定常e u l e r 方程的有限体积解法，并 编制了相应的f o r t r a n 程序。采用简便的代数法生成c 型二维网格和c h 准 三维的定常结构网格，用一种快捷弹性变形技术生成非定常动态网格。在定常的 计算中，采用5 步r u n g e k 删a 时间推进，引进人工耗散，应用了当地时间步长、 焓阻尼等加速收敛技术。在非定常计算中，本文同时采用当地和虚拟时间步，即 双时间推进，在每个真实时间步上采用虚拟时间进行定常迭代，从而不但使时间 迭代具有二阶精度，并使计算时间大大缩短，这样程序与定常程序在算法思想上 一脉相承之求解操纵面的不连续面上的插值时，基于d u k o w i c z 的理论和r a m s h a w 的方法发展了一种应用范围更广阔的插值算子算法。三维定常、非定常的计算结 果表明本文的方法是有效的。j 、一 三维机翼定常、非定常气动力n s 方程计算 a b s t r a c t t 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 ef o rt h et h r e ed i m e n s i o n a l s t e a d y a n d u n s t e a d ye u l e re q u a t i o n si s s t u d i e di nt h i st h e s i s t h ea l g e b r a i cm e t h o di su s e dt o g e n e r a t ect y p e2 d a n dc ht y p e3 d s t e a d ys t r u c t u r e dm e s h e s t h er a p i dd e f o r m i n g t e c h n i q u e sa r ed e v e l o p e dt og e n e r a t eu n s t e a d ym o v i n gm e s h e s i no r d e rt oi n c r e a s e t h e c o m p u t a t i o ne f f i c i e n c y , a r t i f i c i a ld i s s i p a t i o n v a r i a b l et i m es t e pa n de n t h a l p y d a m p i n ga c c e l e r a t i o nt e c h n i q u e sa r eu s e di nt h ef i v e - s t e pr u n g e k u t t at i m es t e p p i n g f o rt h er e q u i r e m e n to f i n s t a n t a n e i t y ，a l li m p l i c i tr e a l t i m ed i s c r e t i s a t i o ni su s e d ，a n d t h ee q u a t i o n sa r ei n t e g r a t e di na f i c t i t i o u s l yp s e u d ot i m e t h i sa p p r o a c ha l l o w st h e r e a l t i m es t e pt oh a v et w o - o r d e ra c c u r a c ya n da c c e l e r a t e st 1 1 e s p e e do fc o n v e r g e n c e b a s e do nt h et h e o r yo fm r j o h nk d u k o w i c za n dt h ea l g o r i t h mo fm r j o h nd r a m s h a w , w eh a v ed e v e l o p e da na l g o r i t h mo ft h es o c a l l e d i n t e r p o l a t i n go p e r a t o r f o ri n t e r p o l a t i o na c r o s st h ed i s c o n t i n u o u ss u r f a c eo f t h e c o n t r o ls u r f a c e t h er e s u l t so f t h r e ed i m e n s i o n a ls t e a d ya n du n s t e a d yf l o w ss h o wt h a tt h em e t h o d s d e v e l o p e di nt h e t h e s i sa r ee f f e c t i v e k e y w o r d s ：c o n t r o ls u r f a c e ，e u l e r e q u a t i o n s ，f i n i t ev o l u m em e t h o d ，u n s t e a d yf l o w 2 南京航空航天大学硕士论文 1 1 问题的背景 第一章：绪论 随着现代航空航天技术的高速发展，实际飞行中的飞行器周围的流场变得越来越 复杂。复杂流动现象的主要研究手段是实验分析法和数值分析法二者。实验分析法的 结果比较真实，是检验理论结果和计算结果的重要标准，但存在着耗资昂贵、周期长、 模拟条件受限、不能显示和测量所有的流动现象、不能同时满足几个相似准则、有测 量误差等不足之处。而对于描述复杂流动现象的数学模型方程，用解析方法求解是不 可能的。因此，在六十年代产生了用数值方法模拟和分析流动现象的学科计算流 体力学( c f d ) 。 自c f d 出现以来，随着其研究工具电子计算机的飞速发展和各种稳定、精 确、快速算法的出现，c f d 如今已形成- - i 3 独立的学科，显示出解决科学理论和工 程实际问题的巨大作用，目前在航空、造船、气象、海洋、水力、液压及石化等工程 领域都有广泛的应用。如今它已和风洞实验、飞行试验一起成为航空航天飞行器设计 流动机理研究的三种相互支持、相互完善的手段。可以说，其在工程设计、流体流动 现象的研究中的作用已处于不可替代的地位。 、 计算流体力学以理论流体力学和计算数学为基础，是这两门学科的交叉学科。主 要研究把描述流体运动的连续介质数学模型方程组离散成大型代数方程组，建立可以 在计算机上求解的算法。另外也可以从流动现象出发，直接建立满足流动规律的、适 定的离散数值模型，而不必经由已有的流体力学偏微分方程组来得到。但是由于理论 流体力学给出的数学模型方程组十分成熟，所以国内外主要应用它们作为研究的基 础。通过时间、空间的离散，把连续的时间离散成间断有限的时间，把连续介质离散 成间断有限的空间模型，从而把偏微分方程转变成有限的代数方程。因此，数值模拟 的实质就是离散化和代数化。离散化：把无限信息系统变成有限信息系统；代数化： 把偏微分方程变成代数方程。而离散的数值解一般可用两种形式给出：网格点上的近 似值，如差分法；单元中易于计算的表达式，如有限元法、边界元法。 如前所述，计算流体力学的应用范围是很广泛的，并且有着很广阔的应用前景， 概括地说它的主要研究方向有：揭示流动机理、处理真实外形、提高数值计算的效率 和稳定性、计算结果的后处理问题、计算程序的通用性问题。计算流体力学( c f d ) 应用到航空领域具体为计算空气动力学，它既有助于我们理解复杂流动的物理本质， 同时也是对风洞实验的一个补充，可以大大缩短飞机设计和研制的周期，降低整个飞 机的设计费用。 而非定常空气动力计算是动气动弹性计算中的重要组成部分。从飞机的发展史 三维带操纵面机翼e u l e r 方程非定常气动力计算 看，从固定翼飞机问世的第一天起，就遇到了气动弹性问题。随着飞机飞行速度的提 高，以气动弹性要求为设计条件的部件越来越多；而且，气动弹性要求不仅只影响到 构件的尺寸，还影响到结构布局，以至影响到空气动力郁局的选择。这样，从飞机的 总体设计阶段开始，就要求进行详细的气动弹性分析。 跟气动弹性一起成长的非定常空气动力计算的发展，一方面取决于飞机研制的需 要，同时还受到计算机能力的限制，特别是跨音速飞机的非定常气动力计算的发展。 二十世纪五十年代，跨音速飞机就已经成为现实。但是，受计算机能力的限制，只能 从亚音速和超音速两边向跨音速范围“逼近”，再加上大量的跨音速风洞实验验证； 随着计算机的发展，到了七十年代，跨音速非定常空气动力计算成为可能。对于要求 更高灵活性的现代飞机来说，气动弹性的设计要求显得尤其重要，相应的对非定常气 动力计算的方法和效果也提出了更高的要求。飞机的柔性部件之间的流动干涉现 象，会降低飞机的性能，甚至带来危险。例如，大后掠翼飞机会产生由涡诱导的气弹 性抖动，以及在跨音速区域内由于激波的出现和移动所引起的各种非设想的气弹性现 象。 操纵面气动载荷的精确预测对改善飞机的机动性能、起飞和着陆性能和稳定性， 优化选择助力器的功率和操纵方式，实现良好的操纵性有很大帮助。随着颤振主动控 制技术发展，机翼操纵面非定常气动力的求解成为关键。尽管跨音速时这一问题较为 困难，但由于颤振边界往往出现跨音速凹坑，发展较精确的跨音速机翼操纵面非 定常气动力计算方法更显迫切。 对于非定常计算，目前国内常采用的是跨音速小扰动位势方程的有限差分解法， 该方法基于线化理论，可计算的衰减频率和攻角有限；虽然说非定常全位势方程的求 解进一步提高了非定常计算的实用性，但是为了弥补它在强激波情况下的不足，还需 要进一步采用e u l e r 方程和n s 方程，以使结果更准确。粘性有时对结果影响很大， 基于n s 方程的跨音速非定常气动力计算方法理应最为精确。但对三维问题，该种 方法往往仅能在巨型机上实现，有时还不一定很理想| 2 l 。本文正是针对于此，研制了 非定常e u l e r 方程的解法，并编制了相应的程序。 1 2 本文的工作 本文针对带操纵面机翼，用有限体积法和双时间推进求解了三维非定常e u l e r 方 程，在剪刀口处基于j o h nk d u k o w i c z 的理论【3 】和j o h nd r a m s h a w 对他的理论的发 展”1 提出了种新的插值算子算法，并研制了相应的计算程序。文中对非定常过程进 行了有效的修正和改进，提高了整体计算效率。其具体工作有如下三个部分：c h 型 2 南京航空航犬大学硕士论文 ( 三维) 弹性网格生成、用双时间进行非定常推进、插值算子的计算和应用。在以上 工作的基础上，编制了“三维带操纵面机翼非定常e u l e r 方程的计算程序”。此程序 计算准确，效率高，所得结果已经验证了这一点。 三维带操纵面机翼e u l e r 方程非定常气动力计算 2 1 引言 第二章：静态和动态网格生成 对于机翼及其周围流场的网格生成，本文采用了一种代数方法及椭圆型方程线 松弛迭代优化方法直接生成准三维单区c h 贴体网格。动态网格生成是非定常气动 力的一项关键技术。利用结构网格的特点，本文采用种快捷的弹性技术生成动态网 格，为操纵面偏转及其非定常气动力计算打下了基础。为了便于操纵面的运动，本文 将二维网格生成的思想和方法利用插值直接推广到了三维，故下面仅仅介绍二维网格 生成的思想和方法。 2 2 静态网格代数法生成 所谓代数生成方法，也就是利用已知的边界值采用中间插值方式来产生网格系的 方法。给出一个包含某些待定系数的特殊曲线坐标函数f ( 孝，叩) ( 即插值函数) ，根据 拟合边界上或某些中间点上给定的笛卡尔坐标值，以及可能还有导数值( 根据各自方 法的特点而异) 的要求来确定这些待定系数。这样，在曲线坐标的常值处即可用此捅 值函数来确定笛卡尔坐标值，从而确定整个网格系统。因此，插值函数是代数方法的 核心。一维插值函数的一般形式可写成： 其中函数仉是满足 r ( 善) = ( 导) l ( 2 1 ) 月i l j 纸( 孚) = 屯 ( 2 2 ) 的多项式，称为型函数( “b l e n d i n g ”f u n c t i o n ) 。若采用l a n g r a n g e 多项式，则 州争彝畿洲圳 s ， 采用不同形式的插值函数，即可形成不同的代数方法。多维时需要多方向的插值，这 里介绍“超限插值”( t r a n s f i n i t ei n t e r p o l a t i o n ) 的思想”1 ，即二维时在每一个方 妻室堕窒堕盔盔兰堡兰笙苎 向单独采用多项式插值： r ( 告，叩) = 吼( 专) r ( 善。，叩) “ 1 ( 2 4 ) r ( 毒，叩) ：兰九( 导) r ( 弧，) 两式相加，则右端项为： s ( 善，叮) ：兰钆( 萼) r ( 六，叩) + 艺丸，( 号) r ( 孝，)s ( 善，叮) = 钆( 号) r ( 六，叩) + 丸，( 号) r ( 孝，) 在孝：o 边界处s ( o ，叩) - - - - r ( o ，7 7 ) + 兰丸( 导) r ( o ，) 。 这表明简单地在各方向单独采用多项式插值不能在边界处拟合给定值。为此超限 插值函数应取为： r ( 善，叩) ：兰( 毒) r ( ，叩) + 芝九( 号) r ( 孝，) 一兰兰吼( 毒) 丸，( 号) r ( ，) ( 2 5 a )r ( 善，叩) = ( 等) r ( ，叩) + 九( 二 ) r ( 孝，) 一吼( 号) 丸，( 号) r ( ，) ( 2 l i m i i o 月1 im z l 4o 这样可以在各边界处都拟合给定值。右端最后一项称为插值的张量一乘积形式。上式 也可改写成： r ( 孝，叩) ：兰丸( 号) r ( 掌，。) + n 纯( 事) r ( ，叩) 一m 。了r ) r ( 茧，巩) 、(25b)r 5 b ) r ( 孝，叩) = 丸( 号) r ( 掌，。) + 纯( 号) r ( ，叩) 一：乙。了) r ( 茧，巩) 、( 2 m l l u h 。j 1 m l l o 具体计算时： f l ( 孝，叩) ：兰庐。( 导) r ( 掌，叩卅) ：叩方向的插值 f 2 ( 孝，叩) = 仇( 号) r ( ，叩) 一f l ( ，叩) 】：孝方向的插值 ( 2 6 ) n - l 1 r ( 手，叩) = f l ( 孝，玎) + f 2 ( 孝，刁) 考虑到只给定两个边界面的函数较难控制所需要的网格的分布，且在壁面附近网格线 的正交性对流场的计算结果有较大的影响，因此可采用广义的超限插值1 6 j ，即不仅给 出函数本身，还给出外法线方向的导数，从而有效地控制在壁面附近插值函数的特点。 2 3 椭圆型方程生成方法 由求解椭圆型方程来实现网格生成的方法称为椭圆型微分方程生成法。这种方法 三丝整堡丛亘盟墨兰! ! ! ! 查堡j ! 塞堂墨垫查生墨 的最早代表为著名的t t m 方法。该方法取x ，y 平面中一对l a p l a c e 方程： 孝麒+ 毒= o ( 2 7 ) 叩。+ 玎w = 0 的解作为贴体坐标系，该解应满足边界条件 孝= 善t ( x ，y ) 【x ，y 】厂： ( 2 8 a ) r = ，7 l 孝5 参( x ，y ) 【x ，y 】r ： ( 2 8 b ) 1 = | | 2 其中厂。和r ：分别为物面边界和外边界，玑和r ：为任意给定常数，卣和彘分别为沿广 和r ，的任意选定的单调函数。 和( 2 8 ) 改写成： 缎嚣一：竺向：”：。 ；! 笔乏：y ；托 。：。， a y 嚣一2 砂如+ 口_ = 0 ；”，。1 边界条件为： 舻乏般j ( 川r l ( 2l0 a ) y = ( 孝，叩，) 1 扣自搀仉： ( 善，叩2 ) r2 + ( 2 1 0 b ) j ，2 9 2 悖，r 2j 其中厂，+ ，厂2 分别和r ，r2 对应，i ， 以及g ，g ：分别由物面边界及外边界确 定。将方程用中心差分格式离散，再用松弛法( 或a d i 方法) 逐次迭代即可最终求得 x ，y 。 此方法的优点是：1 所得网格线是光滑的；2 可以处理复杂的边界；3 若在方程 中置右端项( 源项) ，使之成为泊桑方程，即可控制网格线分布的疏密及倾斜；其缺 点是：1 较难给定右端项：2 较难实现内部点的控制：3 计算时间长。 为了控制内部点的分布，控制网格线分布的疏密，可以在l a p l a c e 方程右端置以 p ，q 项。方程写成一般形式后为： v 2 f 。= p ，i = i ，2 ( 2 1 1 ) 南京航空航天大学硕士论文 在实际计算中还可以采用更简单的形式，即考虑各个方向独自的控制，则 v 2 手= g l p 。，i = l ，2 ( 2 1 2 ) 经过变换运算，它们分别变成： 9 2 2 k + g i l 0 _ 一2 9 1 2 0 口+ ( d e t ( i ，) ) 2 ( j d k + q r q ) = 0 ( 2 1 3 a ) 或 9 2 2 ( + p k ) + g l l ( 0 _ + q 0 ) 一2 9 1 2 0 _ = 0 ( 2 1 3 b ) 其中r 为矢量，即x 和y 。 解此泊桑方程不会遇到任何麻烦，困难在于如何根据控制点的分布给出p ，q 项。 通常希望物面附近的网格线与物面正交( 9 0 0 ) 并控制物面处第一条网格线与物面之 间的距离。目前所有的源项控制方法都可归纳成如下计算步骤： 1 在边界处假定源项的初值( 最简单是令p = q = o ) 2 改变在边界处源项值( 即“源项控制”) ，使生成的网格尽量接近所期望的 网格： 3 将物面与外边界处的源项值内插至流场内部网格点处： 4 迭代求解该源项值时的泊桑方程。 重复步骤2 - 4 直到获得所期望的网格。步骤4 称“内迭代”，而步骤2 4 称“外 迭代”。注意只在外迭代过程中改变源项值，而在内迭代时保持源项值不变。目前常 用的源项控制方法大致可以分成两类，一类根据正交性和网格间距的要求直接导出右 端项p ，q 的表达式；另一类在迭代过程中根据源项变化的情况，采用“人工”控制而 实现所期望的网格。本文中采用的h i l g e n s t o c k 方法属于后一种。下面简单介绍此 方法。 令物面网格间距为 d = i r ( 孝，叩+ 1 ) 一r ( 善，叩) i ( 2 1 4 ) 要求的网格间距值为d ，若d ，一d = 0 ，表示已达要求，则源项q 值无需变化。若 d ，一d 0 ，表明间距太小，需要增大q 值；反之，d ，一d 0 ，则应减小p ：若口，一口 0 ：则应增大p 。故亦可简单地令 p i 。) = p 【2 1 + ( 口，一口) ( 2 1 8 a ) 或取 p ( 2 “) = p c n 一盯a r c t a n ( a ，一口) ( 2 1 8 b ) 在外边界可做类似的控制，只需将阻尼系数改一符号。 w h i t e 分别用t a m ( 堡二旦) 4 和t a i l l l ( 1 2 善) 4 函数来代替方程( ) 和中a r c t a n 函 口，“ 数，其中取0 8 。计算表明这样可提高收敛速度4 0 ，大大改进了原有的方法【8 | 。 2 4 操纵面的偏转和动态网格生成 为了避免方程复杂化，偏转操纵面来进行非定常计算时采用远场固定、物面随动 的动态计算网格。由于每一时间步重新生成计算网格过于费时，在主要用于气动弹性 计算的小幅运动情况下，普遍采用以上一时间步为基础解静平衡方程来得到新的计算 网格。但这一方法仍不经济。 考察绕圆柱的二维不可压位流，令圆柱半径为a ，来流速度比，则流场流函数 眦愕户 汜- 。， 放大圆柱半径至a + e ，则过( r ，e ) 的原流线将产生径向位移r ，经推理可得： 这一结论直接应用于动态网格生成较为困难。经进一步的简化处理，三维可压非定常 流的动态网格点坐标( 以下标u 表示) 可用下式确定 血一蜡a r 。+ 。 南京航空航天大学硕士论文 x 一。= i ，一( i ，一王。) g ( 2 2 0 ) 其中下标s 表示作为初值的静态网格点坐标，下标r 表示全流场静态网格点随物面作 刚性运动的瞬间坐标值。而g 值设为网格点序号的函数，其中下标b 表示对应的物面 点，下标f 代表远场点。 一a x f 等 2 i 篇 2 ) ( 篇 2 z ， 对于二维情况只需去除一个方向。 2 5 小结 网格生成的质量直接决定计算结果的准确性，本文采用的代数方法生成的定常网 格具有高效、简洁而且易于控制等优点，高质量动态网格的生成是非定常计算的基础， 用上述方法生成的静态、动态网格通过后面的实例说明了其有效性。最主要的是，为 了便于操纵面的控制，本文采用了具有二维特性的准三维网格。 三丝萱堡塑耍垫墨! ! 坚立堡斐塞芏墨垫垄生兰 第三章三维非定常e u l e r 方程的求解 通常，非定常计算是在定常计算的基础上进行的。由于非定常计算的时间步 = = 在 整个流场中应取最小时间推进，本文采用双时间推进，即在真实方程中人为加入虚拟 时i 白j 步长，可以大大加快收敛和计算速度。具体介绍如下： 3 1 控制方程 考虑直角坐标系下的非定常e u l e r 方程 a _ w + o f + 塑+ 塑：0 a t a ) c 匆跣 其中： j = p u p u u + p 6 b 伯+ p ) u + x , p w = p 毋 p w o e p v p v u g = ip v v + p p v w 0 + 力v + y , p p w d “ q = f p 矿v 硼w + p b + p ) w + z t p ( 3 1 ) ( 3 2 ) 这里： p = ( r - 1 ) p i e i ( u + y 2 + w 2 ) 】 ( 3 3 ) 在以上的等式中，p 为压强，p 为密度，e ；肛，u 、v 、w 为直角牛标系下速度的三个 分量，u = u x ，；v = v y ，；矽。= w z ，；x t , y ，z ，为网格点在x ，y z 方向上的运动 速度。 将方程3 1 无量纲化，令p 。、几分别表示自由流的压强和密度，无量纲变量定 义为( 右上角加星号的表示无量纲变量) ： 南京航空航天大学硕士论文 c x + ，z 卜( 苎c ，詈，p 2 尝；p 。去2 尝 c cp，p。 c p 。 ( 3 4 ) 所以无量纲的自由流值为 p ：1 ：疋：1 ：p ：士+ 丢m ：y ； y lz u 2 = , y m 。c o s 口；v ：= o ；w ：= 歹m 。s i n a ； ( 3 5 ) 其中口为迎角，y 为比热比，通常其值为1 4 。所得的无量纲化e u l e r 方程形式和( 3 1 ) 完全一样，只是将原方程中的变量用无量纲变量表示即可。为了书写的方便，以后的 表达式中将略去无量纲变量右上角的星号。 3 2 方程的离散及其有限体积格式的建立 有限体积法即积分形式的有限差分方法，直接对每一个微元应用方程的积分形式 利用散度定理将空间的体积积分化为六面体六个面上的面积分，再利用格子的微小化 采用近似积分法，变为面上直接通量的加减。在网格生成部分，己知采用结构网格， 有限体积法最吸引人的特点是它适用于任何坐标系，它不需要进行全场坐标变换，而 只需了解网格单元八个结点的坐标：也就是说不必为计算坐标变换矩阵系数而构造局 部曲线坐标。实际上，与其相当的物理量可由几何原理来严格地得出。 那么将e u l e r 方程离散化，首先要将计算空间划分为六面体( 图3 1 ) 的集合。有 限体积有格心格式和格点格式两种形式，前者把守恒变量定义在网格单元中心，后者 把守恒变量定义在网格点上。两种格式没有本质区别，只是变量定义不同。本文采用 格心格式( 彬。视作网格单元( i ，j ，k ) 格心处的守恒量) 的有限体积法，将e u l e r 方 程的积分形式： 鲁。，p + 。厉丽 c s e ， 用于一个有限体积( i ，j ，k ) ： 一儿一几e一儿几 。一仨上仨 三丝堂塑塑亘塑墨! ! ! ! ! 垄堡j ! 塞堂皇垫垄生簦 8 i + l ，j ，k 图3 1 六面体单元示意图 得到方程的离散形式为： 丢( 。) + g ( ) 舭= o ( 3 7 ) 此处，h y k 代表单元( i ，j ，k ) 的体积，q 为通量项。体积的计算方法如下： 每个六面体可由五个四面体组成，每个四面体的体积可写为： 。= 吉匡 ( 3 8 ) 这里，正。的下标是指这介四面体的四个顶点在六面体中的编号为1 , 2 ，3 和6 ( 如图 3 1 ) 。六面体的体积就是这五个四面体的总和： h = 正2 3 6 + 五8 6 7 + 正3 4 8 + 五8 1 6 + 五6 ” ( 3 9 ) 即使网格变换是奇性的，如：单元的一条边退化为一点，一个面退化为一条线或一个 点，或者它的几个顶点共面，这些几何处理仍然可以得出有意义和准确的面积与体积， 并且相应的数值解不需特别的处理。 3 3 双时间推进 非定常的计算需要时间的最小推进，这是因为每一步推进都有相应的网格变化， 而在每一次网格相对固定时解的推进要有同步性，不允许采用当地时间使解的推进不 同步。而统一采用同一步最小时间推进在三维e u l e r 方程将会非常漫长和不经济。因 此我们对格式进行改进，首先对方程( 3 7 ) 进行隐式离散。隐式离散的方程为： 鲁( 略1 ) 删。) = 。 ( 3 1 0 ) 1 2 南京航空航天大学硕士论文 对时问的离散采用向后k 阶精 a - i ：嗡。) = 嗡1 一 本文采用二阶时间精度，方程( 3 1 0 ) 变为： 堡塑兰堡鉴盟+ 尺( 略，) ：o 2 a f 、”7 将方程( 3 1 3 ) 的左边记为：r ( 略。) ，则方程( 3 1 3 ) 变为： ( 3 1 1 ) ( 3 1 2 ) ( 3 1 3 ) r + ( 雠“) = 0 ( 3 1 4 ) 在此方程的基础上迭加虚拟时间f ，得： 啄1 百d w j t + 胄( 嘴1 ) = 。 ( 3 1 5 ) 方程( 3 1 5 ) 在每步实际时间t 固定后相当于虚拟时间r 上的定常e u l e r 方程，只不过 在残值胄+ 中加入了真实时间项。方程( 3 ，1 5 ) 收敛的结果是便月变为零，也就是略2 在虚拟时间步上趋于稳定，即每一次虚拟时间的推进是为了当地真实时间步的稳定。 而在方程( 3 1 5 ) 的推进中，就相当于求解定常e u l e r 方程，因此定常e u l e r 方程的解 法称为非定常问题求解的一个相当关键的组成部分。那么我们就来介绍定常方程的求 解了。 3 4 定常方程的求解 3 4 1 通量项q 的求法 在结构网格中，由于单元基本上是一个长方体，而且是微元，因此用近似的积分 法，即用每个面上的值与该面面积的乘积代表该面上的通量对于如图3 1 的单元 ( i ，j ，k ) ，有如下表达式： g ( 矿) 社= 岛t r 直残和h积 ? 体 k 的 j 一阳 “ 。州 层 。一址 问 = 时 d 一西 + n表 代 十 n ： 标有上分 ： 中差 里 式度 这 三维带操纵丽机翼e u l e r 方程非定常气动力计算 = 户蜃i ，+ ；小一fj i ，一；。 + 户j l + ；，。一户i l 一；。 + f 蜃l + + ；一f - i “，。一； 3 ，1 6 其中：f = 貂，向量乱弓是单元c t ，“，和c 卜，* ，的交接面在t 方向上的法向 3 4 2 耗散项 在本文采用的有限体积格式中，相邻格子交接面处的数值通量采用平均计算，这 种格式对于矩形长方体网格相当于二阶中- t l , 差分。而中一t l , 格式的一个重要特点就是 对线性和非线性问题均有奇偶不关联性，这种特性使得数值格式不稳定。理论上讲， 物理粘性项是能够提供使数值格式稳定及捕捉到间断区所需的耗散性质的，但是即使 在n s 方程的计算中，那也需要流场内的网格很密，计算量也会相应增大。所以实际 应用中，无论是n - s 方程还是e u l e r 方程的计算中，仍有引入人工耗散项的必要，以 克服中心差分近似所固有的奇偶不关联性和抑制激波附近解的振荡，来保证此格式的 稳定及更有效地捕捉到间断区。j a m e s o n 等人为中心格式设计了一个二阶差分和四阶 差分项相融合的自适应耗散项，即在压力梯度较大( 激波附近) 的区域只采用二阶耗 散项，而在压力梯度较小( 光滑区) 的区域则采用四阶耗散项。 公式( 3 。7 ) 可改写为： o ( 矗形) + q ( w ) - d ( w ) = 0 ( 3 1 7 ) 口f 其中q ( w ) 是通量项，d ( w ) 是耗散项。有： d ( w ) = d ，( ) + d ，( 矿) + d ：( ) ( 3 1 8 ) 这里，d ，( r e ) ，d y ( 矽) 和d ：( ) 为对应于三个坐标方向的耗散项，可写为守恒形式： d ，( 矽) _ dl ，一d1 ( 3 1 9 a ) + j 0 4一j 0 4 d ，( 矿) = d1 一d l ( 3 1 9 b ) o j + 矿”j 4 d ：( 矿) = dl dl ( 3 1 9 c ) j 4 + jo j j 右边项有下面的形式，例如： 南京航空航天大学硕士论文 其中 4 小。3 肚噬拈正肚一小t 甜t c s 2 0 ) 占。彬 = 彬扎小一彬 。 k 彬 女= 彬+ 2 卅一3 彬“+ 3 彬女一彬吐； a 哇。= 圭c 畿+ 瓦h s 1 , j , k ， c sz t ， 这里，f + 是c f l 取1 时的时间步，h 是单元体积。系数占2 和占4 用来适应当地流动 的压力梯度，取： 占量。= 七啦m a x v u ，女) l + ：，# 占哇( 4 。= m a x 卜q 勘) 占哇。p l t ) v 。，。用来探测当地压力梯度 i pj + l 山一2 p u ，+ p j 1 ，j 女l l p 山，女+ 2 p u ，t + p 山l ( 3 2 2 ) ( 3 2 3 ) ( 3 2 4 ) 耗散系数女( 2 和4 通常取为： 一= o 2 5 o 5 ，= 百1 一面1 ( 3 2 5 ) 由上述公式可看出，激波以外的流动光滑区，v = 0 ( 缸2 ) ，则占2 很小，且小于，使 s 4 0 为一阶，d 为三阶；而在激波区，v = 0 ( 1 ) ，占2 ( 女“，则占4 = 0 ，四阶耗散自 动关闭，以二阶耗散为主。同时，耗散项边界的处理也需要注意。 3 4 3 显式时间推进 在求解e u l e r 方程的过程中，目前常用五步r k 时间推进，即 w 0 ) - ( w ) ” 三维带操纵面机翼e u l e r 方程非定常气动力计算 其中 吒“川t 务f ( w ( 0 ) ) w ( 2 ) - w ( ”叫z 务( w ( ) w 【3 ) = w ( ”吣，等( 2 ) 4 】= “叫t 等( w f 3 】) w ( w ( “叫，等趴w ( 4 ) ( w a ，+ 1 ) “= w d o = d l = d ( w o ) ( 3 2 6 ) d 2 = d 3 = p d ( w 2 ) + ( 1 一p ) d o ( 3 2 7 ) d 4 = r d ( w 4 ) + ( 1 - r ) d 2 2 0 5 6 ，y = 0 4 4 口l = 1 4 ，a 2 = 1 6 ，口3 = 3 8 ，口4 = 1 2 ，口5 = l 式中上标m 为虚拟时间推：i 茳m a r 。人工粘性只是在第一、三、五步计及，其余各 步冻结，每一步耗散前的系数作了调整。 3 4 4 边界处理 在e u l e r 方程的求解中，由于计算区域不可能无限大，因此在计算中必然会遇到 两个边界，其一为翼型表面，称为物面边界：其二为远场计算边界，称为远场边界。 3 4 4 1 远场边界 边界条件要求数学上满足适定性，物理上具有明显的意义。本文提到的远场边界 实际上是人为的有限边界，但是流场中的扰动会传到很远的地方。因而远场边界的处 理很重要，直接关系到方法的成败。而它本身的情况又比较复杂，不能直接给定具体 童室堕窒堕墨查堂堡主适茎 一一 的流场值，需要与流场内的值共同来确定。在本计算模型中，远场边界相对于来流保 持固定不动，采用无反射边界条件，即扰动波不会反射回流场。高维问题的远场边界 条件的给定，目前尚无严格的理论依据，因而往往采用一维类推，即把每个方向均看 成一维问题，直接使用一维理论。 为此，沿流场边界，作一维特性分析。现用n 表示法向，下标o 。表示自出流值， 下标e 表示从邻近边界的单元外插得到的边界上的数值。又令虿，i 和c 为速度向量、 单位法向向量和音速。 只考虑边界上的无粘非守恒气体动力学方程： 望+ 彳望：0 ( 3 2 8 ) 研ax 其中： q = 爿= u c o“ 00 0 p c 2 00 0 l p “0 c“ ( 3 2 9 ) 解出矩阵a 的特征值为旯= “，“，“+ c ，“一c ，这四个特征值确定了g ，f ) 空间特征线的 斜率，记： 丑= “一c ，如= “+ c ( 3 3 0 ) 它们对应的两簇特征线为： 宰：“一c ；c - 特征线 讲 拿；“十c ；c + 特征线 ( 3 3 1 ) d f 由数学知识可知，“c 沿特征线a = “c 不变，称“c 为r i e m a n n 不变量， y ly l 用r 表示，记： 1 r 一= “一二c ：r + ：z f + _ = o c( 3 3 2 ) y l，一l 这里，r 一沿入流特征线c 一是常数，可以用来流条件计算得到；r + 沿出流特征线c + 是常数，可以用流场内部外向插值计算。对一般亚音速问题，记 月+ ：圪。+ 垒等，r 一：+ 三等，由此，我们可以求得边界上的法向速度k 和音速。， y ly l 即： 三维带操纵面机翼e u l e r 方程1 f 定常气动力汁算 = 扣谢) c = 孕( r + - r - ) ； ( 33 3 ) 当v o 0 为出流边界，可以从来流值或内点值得到边界上的点的 熵s 和切向速度矿。根据r i e m a n n 不变量，按边界附近信息传播的性质把远场边界分 成以下四种情形： ( 1 ) 亚音速入流边界( m 。 1 ，k 1 ， 1 ，圪 0 ) ，无入流特征线，不需要在边界上规定 条件，边界上的值由流场内的值外插确定，有： p = p ! “= u 。 v = v p p = p 。 ( 3 3 9 ) d 根据边界上的吒，c ，s ，k 值和关系式s 一七，可求得的边界上的“，v ，p ，p 值。进一步可 口。 计算出边界上的各守恒变量的值。应用上述无反射边界条件，可以将远场取物体的5 到8 倍，即可获得准确结果。 3 4 4 2 物面边界 在非定常情况下，物面边界条件须在考虑到网格移动速度的前提下，满足物面无 穿透条件，即：让物面速度等于网格速度。如下： 三丝堂堡丛堕堡墨曼! ! ! ! 查型j ! 塞堂墨塾塑盐翌一 p ( ) 俐盯= m i rm i n 0 l p n ? l 芦”。p 盯 ( 3 4 0 ) p 胛#f ( x ，2 ，+ y f n y + 门：) j 3 4 4 3 后缘k u t t a 条件 对位流计算，绕翼型的环量需按后缘k u t t a 条件确定，因此k u t t a 条件是位流计 算的唯一性条件。在e u l e r 方程的计算中，只有满足k u t t a 条件的解才是e u l e r 方程的 解，e u l e r 方程在5 步r u n g e k u t t a 时间推进中，能自动满足k u t t a 条件，无需专- t 3 考 虑。 3 4 5 加速收敛措施 3 4 5 1 当地时间步长 在求解方程的过程中，每个网格单元在稳定性条件下允许的时间推进步长都是不 同的，如果取统一的时间步长，只能取他们的最小值，以确保流场中每个网格单元都 满足稳定性条件。这样在网格最密的地方，允许的a r 很小，使得在网格较粗的地方r 也不能取得大，因而增加了推进步数，降低了收敛速度。为提高收敛速度，采用变步 长方法，使每个网格的f 取当地稳定性条件或c f l 条件要求的允许值，这对最终定 常解是没有影响的。流场内每个网格单元都满足c f l 条件，大大加速了达到最终稳定 态的速度，节约了机时。 按当地c f l 条件可以确定当地时间步长为： 批c f z f 1 - - - 击a t + 壶a tj 慨t z ， ，_fj 其中： 峨2 可可雨再丽 n q e = 日s | ( 3 4 3 ) ( 3 4 4 ) 南京航空航天大学硕士论文 孑= 纠，雪，是单元在孝c 即i ，方向的面积矢量，s ( 2 】，s 为其三个分量a 本文采用 单元在善方向上的两个微元的面积矢量的平均来计算a 以 l ( i ，女) ，2 ( i ，k + 1 ) ，6 ( i ，j + l ，七+ 1 ) ，5 ( i ，+ 1 ，女) 四个点所在面为例，面积矢量j 为： 扛( 6 1 x 5 2 ) = - 吉区 歹 i y 6 1z 6 1 y 5 2z 5 2 ( 3 4 5 ) x ( f ，+ l ，k + 1 ) 一x ( f ，j ，k ) llx ( f ，+ 1 ，k ) 一x ( i ，j ，k + 1 ) i 其中：一6 1 = lj ，( f ，+ l ，k + 1 ) 一y ( i ，j ，| i ) | ，5 j = jy ( f ，+ l ，七) - y ( i ，j ，k + 1 ) l ，j 的三个分量分 lz ( f ，+ l ，k + 1 ) 一z ( i ，j ，k ) llz ( f ，j + l ，k ) 一z ( i ，j ，+ 1 ) j 别为其在f ，j ，k 方向的值，如下： 4 1 j = 一i 1 ( y 6 l 十z 5 2 一y 5 2 + z 6 1 ) ( 3 4 6 a ) 卜：z 于是有 a ( 2 1 ：i 1 ( x 。l * z 5 2 一x ，：+ z 。) ，一= t j ，fz ( 3 4 6 b ) 一生。= 一圭( x 6 1 * y 5 2 - - y 6 1 * x $ 2 ) ( 3 _ 4 6 c ) s j - ) = 去( 爿( 1 。+ 爿? 。) z 卜j j h j j s j 2 ) = 圭( 一遗。+ 爿遗。， s 尸2 i 1 ( 4 _ ：! ：；。+ 4 遗。) ( 3 4 7 a ) ( 3 4 7 b ) ( 3 4 7 c ) ，1 ，a t f 求法类同a t 。 由于虚拟时间f 要受到真实时间t 的制约，因此在对当地虚拟时间的推进中采用 如下修正： 三丝壹塑坐亘垫墨! ! 堕查堡! ! 塞萱墨垫垄生墨 一 屹。叫 - a , 百2 a t ( 3a s ) 其中f 的求法为上述当地时间步长f 的求法。 3 4 5 2 焓阻尼 定常流中，总焓h 沿流线是常数，如果流线均起源于远前方的均匀自由流，这 里总焓h 。为常数，则定常流中总焓将处处为常数，并等于。焓阻尼就是耿当时 当地的总焓与定常值h 。之差，并假设每个守恒变量矿随时涮的变化率里o t 与差 ( h h 。) 成正比，即 _ o w + f l w ( h 日。) ：0 ( 3 4 9 ) o t 口为用户选取的常数，利用该式改进每个r u n g e k u t t a 推进中最后一步的结果。用、 表示r - k 推进中最后一步的结果，则 w 一w + f l a t w ”1 ( h h 。) = 0 ( 3 5 0 ) 解得 “： 鲨 ( 3 5 1 ) 1 + f l a t ( h 一风) 此式仅用于连续方程和动量方程，对于能量方程我们采用： e ”“一p + r ( p ”+ 1 + p - p h 。) = 0 ( 3 5 2 ) 解得 矿1 = 志 ；一肚( ；一瓶) 】 ( 3 5 3 ) 实际计算中口一般取0 1 ，也只对能量方程进行焓阻尼。焓阻尼的使用，对加速收敛 有一定的效果。 南京航空航天大学硕士论文 第四章通量插值 偏转操纵面后，新网格和原来的网格必