（流体力学专业论文）三维机翼定常、非定常气动力NS方程计算.pdf

南京航空航天大学硕卜学位论文 摘要 本文研究了二维、三维定常n a v i e r s t o k e s 方程的有限体积多重网格解法 和非定常n a v i e r s t o k e s 方程的有限体积解法，并在此基础上编制了相应的计算 程序。采用简便的代数法生成c 型( 二维) 、c h ( 三维) 定常结构网格，用一 种快捷弹性变形技术生成非定常动态网格。在定常的计算中，采用5 步r u n g e k u t t a 时间推进，引进人工耗散，在人工耗散中应用了各向异性尺度因子，并应 用了当地时间步长、残值光顺、焓阻尼等加速技术。另外，采用了加速收敛最 有效手段一多重网格法，并用一种简单回插算子克服了粘性时多重网格的插值 困难，使计算稳定收敛。非定常计算中，由于瞬时性要求，必须采用最小时间 步长进行推进，这样将使推进相当缓慢丫本文同时采用当地和虚拟时间步，即 双时间推进，在每个真实时间步上采用虚拟时间进行定常迭代，从而不但使时 间迭代具有二阶精度，并使计算时间大大缩短，而且程序与定常相差不大，使 程序修改简洁、方便。紊流计算采用了b a l d w i n l o m a ) ( 代数紊流模型，文中对 紊流外边界的探测进行修正，改进了边界域的处理，加速了迭代过程。 ，二维、三维定常、非定常的计算结果表明本文工作和方法是有效的。 奉文的改进和有新意的工作在于： 1 在多重网格计算过程中，由粗网到密网的回插算子采用原值回传，克服了 线性插值对网格的依赖性，不仅节省了每个多重网格循环过程，而且提高 了效率，加速了整个迭代进程。 2 在非定常过程中，采用双时间推进，使时间推进具有二阶精度，提高了程 序高效性，而且加速了非定常计算速度。 3 在紊流模型的控制过程中，、雾日孑边界探测，远场无须紊流模型的介入， 这种方法节省了迭代时间。， 关键词：有限体积、非定常、双时间、多重网格、紊流 三维机翼定常、非定常气动力n s 方程计算 a b s t r a c t t h ef i n i t e v o 】u m em u 】t i g r i ds c h e m ef o rt h e “v oa n dt h i c ed i m e n s i o n a ls 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 n d 也ef i n i t e v 0 1 u m ee l e m e n tm e t h o df o ru n s t e a d vn a v i e r - s t o k e se q u a t i o n sa r es 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 og e n e r a t ec t y p e ( 2 d ) a n dc ht y p e ( 3 d ) s t e a d ys t n j c t i l r e dm e s h e s a n dt h er a p i dd e f o 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 丘i c i e n c y ，a n 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 e s t e p ， r e s i d u e s a v e m g i n ga n de n t h a l p yd a m p i n ga c c e l e r a t i o nt e c l l l l i q u e s a r eu s e di nt 1 1 e f i v e s t e p r u n g e - k u t t at i m es t e p p i n g t h em u l t i g r i ds c h e m e ，t h em o s tp o w e m la c c e l e r a t i o n t e c h n i q u e ，i ss u c c e s s 凡l l yu s e d i nt h ev i s c o u sn o wt om a k et h ec o i n p m a t i o n c o n v e 嚷e s t e a d i l y f o rt h er e q u i r e m e n to fi n s t a n t a n e i t y ，t h em i n i m a l t i m es t e pm u s tb eu s e di n t h e u n s t e a d y n o w c o m p u t a t i o n ， s oi tm a r c h e s s l o w l y a n i m p l i c i t r e a l t i m e d i s c r e t i s a t i o ni su s e d ，锄dt h ee q u a t i o n sa r ei n t e g r a t e di naf i c t i t i o u sp s e u d ot i m e t h i s a p p r o a c ha l l o w st h er e a l t i m es t e pt oh a v et w o o r d e ra c c i l r a c ya i l da c c e l e r a t e st h e s p e e do fc o n v e 唱e n c e a n dt h ep m g r a mi se d i t e dc o n v e n i e n t l y i l lo r d e rt os i m u l a t e t i l r b u l e n tf l o w s ，t l l eo r i g i i l a lv e r s i o n so f t h eb a l w i n l o m a xt i l f b u l e n c em o d e lh a sb e e n e m p l o y e di nt h ep r e s e n ts t u d ya n d t h ei t e r a t i o ni ss p e e d e d u pw i t h t h ei m p r o v e m e n t o f s e n s i n gt ot h eo u t e rt l l r b u l e n tb o u n d a r v t h es i m u l a t i o nr e s u l t so ft h et w o 锄dt t l r e ed i m e n s i o n a ls t e a d ya n d u n s t e a d yn o w s s h o wt h a tt h em e t h o d sd e v e l o p e di nt h et h e s i sa r ee 丘e c t i v e t h ei m p r o v e d 柚dn e wi d e 弱i nm i sp a p e ra r e 鹊f o n o w s ： l 、i i lt h e m u l t i 鲥ds c h 锄e ，t h ed e p e n d e n c eo fl i n e 盯i n t e r p o l a t i o no nm e s hi s o v e r c o m eb y u s 妯g t h eo r i g i n a lv a l u e b a c k 胁s f b r r e d 叩e r a t o l t h i sm e 也o ds a v e s 也et i m eo f t h ee a c hi t e r a t i o na n da d v a n c e s 也ee 丘i c i e n c vo f 也ei t e r a t i o n 2 、i i lt h eu n s t e a d yn o w sc o m p u t a t i o n ，m ed u a l - t i m es t e p p i n gm e t h o di su s e d t h e a c c u r a c yo f t l l es c h e m ei s 咐o o r d e ra 1 1 dm ec o n v e 疆e n c ei sa c c e l e r a t e d 3 、i i lt h eb a l w i n l o m a xt i l r b u l e n c em o d e l ，s e n s i n gt 0t h eb o u n d a r yi su s e d s ot h e t i l r b u l e n c em o d e li sn o tn e c e s s a r vt ou s e di nt h ef h rn o wf i e l da n dt h i ss a v e s 也e c o m p u t a t i o n t i m c k e yw o r d s ：f i n i t e v o l u m em e 也o d ，u n s t e a d yf i o w d u a l t i m es t 印p i n g ，m u i t i - g r i d s c h e m e t i l r b u l e n c em o d e l 2 南京航空航天大学硕士论文 1 1 问题的背景 第一章：绪论 随着现代航空航天技术的高速发展，实际飞行中的飞行器周围的流场变的 越来越复杂。复杂流动现象的主要研究手段是实验分析法和数值分析法二者的 结合。实验分析法的结果比较真实，是检验理论结果和计算结果的重要标准， 但存在着实验耗资昂贵、周期长、受模拟条件限制、不能显示和测量所有的流 动现象、不能同时满足几个相似准则、有测量误差等不足之处。而对于描述复 杂流动现象的数学模型方程，用解析方法求解是不可能的。因此，在六十年代 产生了用数值方法分析流动现象的学科一计算流体力学( c f d ) 。 自c f d 出现以来，随着其研究工具电子计算机的飞速发展和各种稳定、 精确、快速算法的出现，c f d 如今已形成一门独立的学科，显示出解决科学理 论和工程实际问题的巨大潜力，目前在航空、造船、气象、海洋、水力、液压 及石化等工程领域都有广泛的应用。如今它已和风洞实验、飞行实验一起成为 航空航天飞行器设计流动机理研究的三种相互支持、相互完善的手段。可以说， 其在工程设计、流体流动现象的研究中的作用已处于不可替代的地位。 计算流体力学以理论流体力学和计算数学为基础，是这两门学科的交叉学 科。主要研究把描述流体运动的连续介质数学模型离散成大型代数方程组，建 立可以在计算机上求解的算法。可以从流动现象出发，直接建立满足流动规律 的、适定的离散数值模型，而不必经由已有的流体力学偏微分方程组。但是由 于理论流体力学给出的数学模型十分成熟，仍可应用它们作为研究的基础。通 过时间、空间的离散，把连续的时间离散成间断有限的时间，把连续介质离散 成间断有限的空间模型，从而把偏微分方程转变成有限的代数方程。因此，数 值方法的实质就是离散化和代数化。离散化：把无限信息系统变成有限信息系 统；代数化：把偏微分方程变成代数方程。而离散的数值解一般可用两种形式 给出：网格点上的近似值，如差分法；单元中易于计算的表达式，如有限元法、 边界元法。 如前所述，计算流体力学的应用范围是很广泛的，并且有着很广阔的应用 前景，概括地说它的主要研究方向为： 揭示流动机理 处理真实外形 提高数值计算的效率和稳定性 计算结果的后处理问题 计算程序的通用性问题 经过多年的发展，计算空气动力学已经成长为计算流体力学( c f d ) 在航空 领域的一个主要部分。这些发展既有助于我们理解复杂流动的物理本质，同时 也是对风洞实验的一个补充，可以大大缩短飞机设计和研制的周期，降低整个 飞机的设计费用，特别是在气动弹性领域。 气动弹性在飞机的设计和开发过程中扮演着重要角色，对于要求更高的灵 活性的现代飞机来说，这点尤其重要。飞机的柔性部件之间的流动干涉现象， 会限制飞机的性能，从而带来危险。例如，大后掠翼飞机会产生由涡诱导的气 弹性抖动，以及在跨音速区域内由于激波的出现和移动所引起的各种非设想的 三维机翼定常、非定常气动力n s 方程计算 气弹性现象。 气弹性实验需要大量的资金，一个气弹性实验的费用是普通气动实验的许 多倍，而如果采用数值计算来补充风洞实验，那么整个飞机的设计费用就会大 大降低。 在这方面，准确求解跨音速非定常气动力是解决气弹和颤振问题的关键。 目前国内常采用的是跨音速小扰动位势方程的有限差分解法，该方法基于线化 理论，可计算的衰减频率和攻角有限；虽然说非定常全位势方程的求解进一步 提高了非定常计算的实用性，但是为了弥补它在强激波情况下的不足，还需要 进一步采用欧拉方程和n s 方程，以使结果更准确。本文正是针对与此，研制 了非定常n s 方程的解法，并编制了相应的程序。 1 2 本文的工作 本文针对前面所提的问题，用多重网格和有限体积法求解了二维、三维定 常和非定常n a v i e r s t o k e s 方程，在紊流的计算过程中采用了b a l d w i n l o m a x 代数模型，并研制了相应的计算程序，文中分别对多重网格，非定常过程，紊 流过程进行了有效的修正和改进，加快了整体计算效率。具体如下： 1 c 型( 二维) 、c - h 型( 三维) 粘性网格生成 采用简便的代数法生成c 型( 二维) c h ( 三维) 定常结构网格，用一种 快捷弹性变形技术生成非定常动态网格。 2 二维和三维定常和非定常的n a v i e r s t o k e s 方程的求解 在定常的计算过程中采用了j a i i l e s o n 的与时间相关的显式有限体积法和五 步r u n g e k u t t a 推进格式；在多步推进过程中，为了节省计算时间，采用只在 奇数步计算耗散，其它各步用前面的线性插值，并使用了当地时间步长和残值 光顺等加速收敛技术； 3 多重网格加速收敛技术 传统多重网格在回代过程中对网格展弦比的依赖性较大，尤其是粘性网格 的物面附近，常发生插值困难。文中改进了回代的插值过程，用原值回插取代 线性插值，取得了较好的效果，不但加速了整体迭代过程，而且在计算过程中 对网格展弦比的影响基本消除。 4 用双时间进行非定常时间推进 在非定常问题的求解中，如果采用当地真实时间求解n a v i e r s t o k e s 方程 将是非常缓慢和不经济的。因此，对时间格式进行了修正，采用双时间推进， 使计算能够快速和有效的进行。 5 紊流模型的加入 在紊流计算中，采用b a l d w i n l o i i l a x 代数模型。由于紊流在物面附近影响 较大，而在远场其作用却不大，因此为了减少不必要的计算，本文用涡量来判 断紊流的边界，节约了计算时间，加快了收敛速度。 在以上工作的基础上，研制了“多重网格定常、非定常二维、三维层流、 紊流n a v i e r s t o k e s 方程的计算程序”。研制出的程序计算准确、效率高。在翼 型和不同的三维机翼上已经验证了这一点。 南京航空航天大学硕士论文 2 1 引言 第二章：静态和动态网格生成 对于机翼及其周围流场的网格生成，本文采用了一种基于无限内插或多变 量内插理论的直接代数方法。该方法只需在边界上给出定义，并对局部网格质 量有很好的控制。动态网格生成是非定常气动力的一项关键技术。利用结构网 格的特点，本文采用一种快捷的弹性技术生成动态网格，为非定常气动力计算 打下了基础。 2 2 静态网格代数法生成 本文所采用的网格生成方法是一种直接代数生成法，称为无限内插法或多 变量内插法。 由g o r d o n 和h a l l 所提出的无限内插理论是种有关多变量内插的理论， 它的基本思想如下： 令f ( u ，v ，w ) = x ( u ，v ，w ) ，y ( u ，v ，w ) ，z ( u ，v ，w ) ，代表定义在区域 “l “p ，y l v ，w i w w ，内的三个变量u ，v ，w 的一个向量函数。在整个 区域内这个函数是未知的，只是在特定平面有： ，( ，y ，w ) = 吒( y ，w ) ：k 2 1 ，2 ，p 厂( “，v i ，w ) = 瓯( “，w ) ：k 2 1 ，2 ，q 厂( “，v ，w t ) = q ( “，v ) ； k 2 1 ，2 一，r ( 2 1 ) 为了在这些指定的平面之间插值，我们需要定义一系列单变量的融合函数 儿( w ) ；k 2 1 ，2 ，r 它们只需要满足下列条件： ( 砷) = 屯；屈( u ：) = 如；靠( ) = 这里： 6 h = 0 ；七fj h = l ；七= ， ( 2 2 ) ( 2 3 ) ( 2 4 ) 三维机翼定常、非定常气动力n s 方程计算 于是有如下的无限内插函数 工( ，w ) ：杰( “) d 。( v ，w ) = l ( 。，。，w ) ： ( 。，w ) + 圭反( v ) 钆( “，w ) 一 ( “，v 。，w ) t ，l ，( “，y ，w ) = 厶( “，v ，w ) + 以( w ) 【q ( “，v ) 一以 ，v ， ) 】 ( 2 5 ) i = l 现在函数f 就定义了由u ，v ，w 空间内区域“，“。，v isv v 。，w ls w w ， 向x ，y ，z 空间内任何形状区域的转换关系。理论上来说，该格式的每一步都包 含了在三个可能方向( u ，v 和w ) 的一个单变量插值，而且，可以将这三步格式 缩减成一个表达式。可以证明，如果函数f 在各空间平面间是连续的，那么在 三个方向的插值顺序是不会影响插值结果的。 对于一个典型的三维网格生成问题，有必要既减少所需的输入几何数据，又 能保持对插值过程的控制，至少在所关注的表面附近应当如此。如果有足够的中 间表面得到定义，那么，上面所描述的插值过程可以保证任何程度的控制。但 是，没有定义内部表面( f 只定义在参量区域的外表面) ，那么，对插值的控制 将是很糟糕的。为了改进在这种情况下的控制，定义了一种使用函数f 在表面 外向的导数的无限内插方案；定义f 的表面外向导数是为了引入对表面附近的 分布函数的直接控制。给定的数据可写为： 等m ，w ) 电( v ，w ) ；h ，2 n 2 0 ，l ，2 ，仇 斋弛，v ) = 矿b ，w ) ；k - 1 ，2 n 2 0 ，l ，2 ，a 斋m v 咿0 气_ l ，2 n = o ，l ，2 ，p t ( 2 6 ) 这可简化f 以及f 在区域“2 ，v i v v 2 ，w i w 蔓w 2 的外表面的外向 导数的定义。为了将f 插入参量区域的内部，我们定义了一组新的单变量融合 函数： 口：4 ( “) ；k 2 1 ，2n 2 0 ，1 ，2 ，既 4 南京航空航天大学硕士论文 y i ”( w ) ；k = l ，2n = o ，1 ，2 ， 它们只需满足下列条件： 昙气“) = 如k ； 熹聃咖如。； 号棚咿西“ 下面给出了无限内插格式的一般形式： 2 ( 州，w ) = 口乳咖( v ，w ) t = 1 月= o ( 虬v ，w ) = ( “，v ，w ) + 喜弘) ( v 删w ) - 斋脚洲 + ( v ) ( 圳) 一言 ( ) 】 t ；l = 0 v7 厂( ，n 计= ：( 以u 砷 ( 2 7 ) ( 2 8 ) + 喜粪( w ) 【c p ( “，v ) 一专，2 ( “ y m ) ( 2 9 ) 这时，函数f 定义了由u ，v ，w 空间内区域“i “2 ，v l y v 2 ，w i w w 2 到 x ，y ，z 空间内任意形状区域的一个变换。也可以证明，单变量插值在三个方向 的顺序不影响结果 对于c h 型网格，其生成过程可概述如下： 1 ) 定义机翼表面( 对应于参量空间的w = w l 表面) ： 2 ) 计算机翼表面外向导数； 3 ) 定义外边界( 对应于参量空间的w = w 2 表面) ； 4 ) 定义w 方向( 在外表面和机翼表面之间) 插值的融合方程； 5 ) 在w 方向插值( 无限内插的第一步) ； 6 ) 定义对称面( 对应于参量空间的v = v 1 表面) ： 7 ) 定义v 方向插值的融合方程； 8 ) 在v 方向插值。 本文在计算中二维，三维网格均由同一三维网格生成程序生成，只是在 三维机翼定常、非定常气动力n s 方程计算 二维的计算中网格点取三维的一个剖面。因为n s 方程的求解对网格内层的高 度有很严格的要求，因此我们在网格生成程序中加入了对这一高度的控制，并 且使其处于输入参数，我们不需对网格程序进行修改，只需控制输入程序就可 以得到满意的网格。在文后的附图中我们给出了三维网格的一些剖面。可以看 出，该方法生成的网格正交性好，网格贴体。满足n s 方程计算的要求。附图 ( 1 8 ) 中给出了我们生成的网格，从计算效果来看，该网格完全适应计算的要 求。 2 3 网格的几何特性 y 图2 1 ：用8 个角点定义的c a n e s i a i l 坐标系下的六面体 有限体积法最吸引人的特点是它适用于任何坐标系统。它不需要进行全场 坐标变换，而只需了解网格单元八个结点的坐标；也就是说不必为计算坐标变 换矩阵系数而构造局部曲线坐标。实际上，与其相当的物理量可由几何原理来 严格地得出。例如，对于全部所需的1 0 个矩阵系数；每个单元的体积q 和它的 三个表面积s ，e ，& 的三个分量( 图2 1 ) ，如果一个表面的四个角点共面， 那么它的面积为它的对角线乘积的一半s = 三k k ：如果不共面也可以计算 出，只是稍微复杂些。体积的计算方法如下： 每个六面体可由五个四面体组成，每个四面体的体积可写为 ( 2 1 0 ) 这里，互：。的下标是定义这个四面体的四个角点。六面体的体积就是这五 个四面体的总和。 即使网格变换是奇性的，如：单元的一条边退化为一点，一个面退化为一 条线或一个点，或者它的几个角点共面，这些几何处理仍然可以得出有意义和 准确的面积与体积，并且相应的数值解不需在程序中作特别的处理。 ；6 z z z z 儿儿儿儿 而砌崩靠 l 一6 i i 6五 南京航字航天大学硕士论文 7 2 图2 2 ：六面体的体积是5 个四面体的体积总和 q = 正2 3 6 + t 8 6 7 + 互3 4 8 + l 8 1 6 + 1 6 8 5 ( 2 1 1 ) 2 4 动态网格生成 为了避免方程复杂化，非定常计算时采用远场固定、物面随动的动态计 算网格。由于每一时间步重新生成计算网格过于费时，在主要用于气动弹性计 算的小幅运动情况下，普遍采用以上一时间步为基础解静平衡方程来得到新的 计算网格。但这一方法仍不经济。 考察绕圆柱的二维不可压位流，令圆柱半径为a ，来流速度圪，则流场流 涵数为： y ：圪fr 一竺陆毋 ( 2 1 2 ) l， 放大圆柱半径至a + e ，则过( r ，o ) 的原流线将产生径向位移r ，经推理可 得：这一结论直接应用于动态网格生成较为困难。经进一步的简化处理，三维 可压非定常流的动态网格点坐标( 以下标u 表示) 可用下式确定 血：。一。耸冀 i 。= i ，一g ，一i ，) g ( 2 1 3 ) 其中下标s 表示作为初值的静态网格点坐标，下标r 表示全流场静态网格点随 物面作刚性运动的瞬间坐标值。而g 值设为网格点序号的函数 其中下标b 表示对应的物面点，下标f 代表远场点。 对于二维情况只需去除一个方向。 ( 2 1 4 ) 2 5 小结 网格生成的质量直接决定计算结果的准确性，本文采用的代数方法生成 的定常网格具有高效、简洁而且易于控制等优点，高质量动态网格的生成是非 定常计算的基础，用上述方法生成的静态、动态网格通过后面的实例说明了其 有效性。 门 等n 从嚣剖 川u一 南京航空航天大学硕士论文 第三章二维定常n a v i e r s t o k e s 方程解法 3 1 引言 真实的流体是有粘性的，用e u l e r 方程无法模拟由于粘性而引起的二次 分离以及计算飞行器的摩擦阻力等问题；其次在跨音速流动计算中，e u l e r 方 程不能真实反映穿越激波的流体微团的熵增，因而计算的激波位置过于靠后， 致使物面压力分布及相关的升力计算不准：最后在大迎角绕流中，流场中存在 复杂蜗系的相互干扰、旋涡的破碎、小涡的耗散等现象。对于这类流动，如果 采用e u l e r 方程来计算，都无法反映流场的真实特征，因此需要求解控制真实 流体的n a v i e r _ s t o k e s 方程。 就实用外形的复杂流动来说，用n a v i e r _ s t o k e s 方程来计算收敛速度相当 慢，因此需要高效算法。在所有的加速收敛方法中，多重网格是最有效的加速 收敛措施。文中采用了多重网格加速收敛，取得了较高的收敛速度。 3 2 数值模拟 3 2 1 控制方程 考虑直角坐标下的n s 方程 业+ 堑+ 塑：f 塑+ 塑1 a f 氟砂l 苏砂j 其中： w = ，= p 册 p 施 俄 肛2 + p j d “v 辨e + 肾 g = m p 忧 2 + p 斛e + v p ( 3 1 ) ( 3 2 ) 三维机翼定常、非定常气动力n s 方程计算 这里 盯“= 2 删，一亏( 虬+ v ，) 盯，= 2 ，一詈( “，+ v ，) 盯叫= 盯”= ( “y + v ，) 吼：一七婴 戚 d “ “j 。_ 嬲 口v v j2 _ 缎 q 。：一七婴 。 却 ( 3 3 ) ( 3 4 ) ( 3 5 ) ( 3 6 ) “。= 婴 ( 3 7 ) 。 砂 v ，= 娑 ( 3 8 ) 。 却 女：j ，_ 坐y ：1 4 ( 3 9 ) y ln 、。 p = ( ，一1 ) p 【e 一丢 2 + v 2 ) 】 ( 3 1 0 ) p = 肚r ，r 由该式推出。 在以上的方程中，p 为压强，户为密度，u 、v 为速度的分量，t 为温度，为 粘性系数，由s u t h e r l a n d 定理确定，七为热传输因子，p ，为p r a n d t l 数。 将方程3 1 无量纲化，令p 。、p 。分别表示自由流的压力和密度，无量纲 变量定义为( 右上角加星号的表示无量纲变量) ( 石+ ，y + ) = ( 兰，上) ；p + = 旦：p = 旦；f + = 三旦生 cc p。p。 cp 。 ( “+ ，v + ，e ，丁) = g 一 盯p 。叻呦肌 g 一 掣 盯v+ h o 岍 r一叠=一凡e一儿一凡 ，一仨上仨 南京航空航天大学硕士论文 = 老訾。k 。 ：二冬； ( 3 1 1 ) y 一1 只 所以无量纲的自由流值为 p ：= ；p ：= 1 ；e = 击+ 圭蛾y ； “：= 万m 。c o s 口；v ：= 万m 。s i n 口； ( 3 1 2 ) 其中口为迎角，y 为比热比，通常其值为1 4 。所得的无量纲化n a v i e r s t o k e s 方程形式和( 3 1 ) 完全一样，只是将原方程中的变量用无量纲变量表示即可。 为了书写的方便，以后的表达中将略去无量纲变量右上角的星号。 3 2 2 有限体积格式的建立 有限体积即利用n s 方程的积分形式，直接用于每一个曲边格子，将空间 的积分利用格林公式化为格子边上的曲线积分，再利用格子的微小化，采用近 似积分法，变为直接通量的加减，这就是有限体积法。这种方法可以方便的应 用于任意形状的单元而不必引入复杂的形函数。下面详细介绍这种方法。 将n s 方程离散化，首先要将计算空间划分为结构网格( 图3 1 ) 的集合， 采用格心格式，方程离散为： 图3 1单位结构网格示意图 丢( ) + ( q 。( ) + q 。( ) 。+ 。( ) g ) = 。 ( 3 1 3 ) 。代表单元的面积，q c 为非粘性通量项，q ，为粘性通量，d 为耗散项， 下面分别介绍各自的求法和具体做法： 3 2 2 1 非粘性通量项q 的求法 三维机翼定常、非定常气动力n s 方程计算 q = g ( 傍一础) ， 在单元的积分上，由于结构网格单元基本 出 都是由矩形构成，而且其为微小单元，因此用近似的积分法，既用每条边上的 值与边长的积代表通量在该边上的积分。因为每条边都有相应的x 方向的长度 d x 和y 方向的长度d y ，因此有： q = j ( 纳一础) = 弦1 一础l + 肋2 一础2 击 + 肋3 一础3 + 傍4 一础4 ( 3 1 4 ) 其中蛳，出f 分别为单元各边的长度在y ，x 方向的分量，而厂，g 分 别为单元边上的中点通量，该值由相临单元格心值平均得到。 3 2 2 2 粘性通量项q ，的求法 q ，= 一水r 砂一鼢) ，求法和非粘性通量的求法一样，不过r ，s 中 由 包含速度的一阶导数，下面以娑， 优 罢的求法来说明一阶导数的算法： 洲 罢丢謦2 沙= 如嘞z 脚s 哪4 ， 考一；挚一 p 一知础z 恸s 毗。， 其中咖，捌分别为单元各边的长度在y ，x 方向的分量，而为单 元边上的中点速度，该值由相临单元格心值平均得到。 3 2 2 3 耗散项d 由于本文采用的有限体积格式相当于二阶中心差分，所以为了克服中心差分 所固有的奇偶不关联性，以及避免激波附近解的波动，需在当地引入带开关的 耗散项，即在压力梯度较大的区域只采用二阶耗散项，而在压力梯度较小的区域 则采用四阶耗散项。 d ( 矽) = d ，( ) + d ，( ) ( 3 1 6 ) 这里，d ，( 形) ，d ，( 缈) 为两个曲线坐标方向的对应耗散项，可写为守恒形 南京航空航天大学硕士论文 式 d ，( ) = d i ，一d h i j卜j 7 d y ( 矽) = d i dl + j一j 右边项有下面的形式，如 ( 3 1 7 ) ( 3 1 8 ) 其中，旯。为平均尺度因子，其求法为： h i 7 t 。，= 丢眈) u + 仇) ，+ 阮l 小+ “l ，】 ( s 。) 其中谱半径t ，以分别为： 丸= 砧i + c ；丑= v i + c 7 ( 3 2 0 ) “v 为逆变速度，c 为音速，占2 和s 4 是由流动的压力梯度来确定的。 设 那么： v f ，= s 墨，= k 壮m a 】【( v + l l + = j s 遵，= m & x 。，c 七( 哪一s 遵， 七( 2 和尼( 通常取为 ( 3 2 1 ) ( 3 2 2 ) o 棼 三维机翼定常、非定常气动力n s 方程计算 2 ) ：1 t ( 4 ) ：上( 3 2 3 ) 。 3 2 由上述公式可看出，激波以外的流动光滑区，v = 0 ( ) c 2 ) ，则2 很 小，s 4 o ：而在激波区，v = o ( 1 ) ，占4 = o ，则s 2 很大，以二阶耗散为主。 3 2 2 4 各项异性的人工粘性 中心差分格式方法中，人工粘性模型对方法的成功应用起着关键性的作 用。人工粘性拟制解在激波附近的振荡，又阻尼解在光滑区域内的高阶误差， 对解的线性稳定和收敛于定态是很重要的。虽然计算实践表明标量人工粘性模 式对无粘流场的计算是方便的，也比较可靠，但因其在各方向的尺度因子五是 个平均标量，在计算中往往使人工粘性值较大，会使无粘流场中激波的扑捉往 往占到3 4 个网格点：在粘性流场，特别在计算壁面附近的边界层时也会带来 较大的误差，这种不确定在网格不太精细时尤为突出，因为三维流动的计算量 巨大。网格不可能很精细，就更有必要改进人工粘性模型。此外，有些流动区 域( 翼型后缘附近) 虽为连续流动也必须将人工粘性值限制得很小才能有较好 的计算结果，况且标量人工粘性对于高超声速流动的计算结果也很不理想。这 些问题都不是简单地在标量人工粘性模型中减小常数就能解决的，必须改进模 型。 式( 3 一1 9 ) 和( 3 2 0 ) 中定义的尺度因子一般称为各向同性的尺度因子， 此尺度因子对无粘流场长宽比d ( 1 ) 的典型网格是可以的，随着网格长宽比的增 大就会引起太大的人工粘性。一般粘流近壁处网格的长宽比可达d ( 1 0 3 ) ，此尺 度因子会引起过大的人工粘性而极大的降低摩擦应力的计算准确度。改进模型 采用各个方向的平均值带替整体平均。即以： 。：广丢帆) 。+ 阮h 】； 。= 告帆l ，+ k l 。】 ( 3 2 4 ) 代替式( 3 1 9 ) ，这就是各向异性的尺度因子。在采用多重网格技术时，式( 3 2 4 ) 在流向的尺度因子会太大，使阻尼高频误差的粘性太小而大大降低收敛速 度，故再进一步改进流向尺度因子： 钆忆，：去阮1 ，( _ k 】 ( 3 2 5 ) 其中( 五，) u = 纯，( y ) ( 旯，) u ；识( ，) = l + ，己，o 口 1 ，一般取口= 1 2 。 y = 五。乃，而在y 方向，可以定义： ( 一l = “吉) ( 饥， ( 3 2 6 ) 南京航空航天大学硕士论文 本文在编制程序时应用了上述各向异性的尺度因子代替各向同性的尺度因 子，并用了多重网格技术，取得了很好的结果。 3 2 2 5 边界处理 在方程的求解中，由于计算域不可能无限大，因此必然在计算中会遇到两 个边界，其一为机翼表面，我们称之为物面边界；其二为远场计算边界，我们 称之为远场边界。 3 2 2 5 1 物面边界 b l 斗。 图3 1 粘性流体由于粘性的存在，物面上的速度值为零，而温度采用法向梯度为 零，即： “= v = 0 ：订v r = 0 ： 在实际计算中，让壁面温度和附近网格温度相同即可，即如图3 1 3 2 2 5 2 远场边界 边界条件不能随意给定，它要求数学上满足适定性，在物理上具有明显的 意义。远场边界的处理很重要，直接关系到方法的成败。本文采用j 锄e s o n 提 出的远场无反射边界条件，即扰动波不会反射回流场。由于高维问题的远场边 界条件的给定，目前尚无严格的理论依据，因而往往采用一维类推，即把每个 方向看成一维问题，直接使用一维理论。 图3 2 三维机翼定常、非定常气动力n s 方程计算 以一维为例，只考虑边界上的无粘非守恒气体动力学方程 这里，q = 望+ 爿塑：o aa x 爿= “ p o“ 0o 0 傀2 0o 0l 口 “0 c“ ( 3 2 7 ) c 为当地音速 则a 的特，仕值a 2 “，“，“+ c ，“一c 确定了b ，f j 空间特征线斜率，我们可以证明： “c 沿特征线a = “c 不变，称为融e m a n n 不变量。 y l 对一般亚音速问题，我们记彤2 + 鲁，r 一= 吒。+ 鲁，式中，n 表 示法向，m 表示来流值，e 表示由内点上的值外插得到，由此，我们可以求得 边界上的法向速度和音速c ，即 圪= 三( r + 何) ；c = 孕( r + - 州： ( 3 2 8 ) 根据法向速度的符号( 当 o 为出流边界) ，可以从来流 值或内点值得到边界上的点的熵s 和切向速度矿。 可分为如下情况： ( 1 ) 亚音速入流边界( m 。 l ， 1 ， l ， 0 ) ，有 p 2p 。 “= “f v = k p = 见( 3 3 2 ) 根据边界上的，c ，j ，巧值和关系式j 七，可求得的边界上的“，v ，p ，p 值。进 d p 一步可计算出边界上的各守恒变量的值。 3 2 2 5 3 边界处理的其他注意问题 在实际计算中，计算域的边界是有限的，因此在亚音速时，物体产生的 扰动会传播到离物体很远的地方，对于二维升力问题，如果远场边界取得不太 远，则需要考虑环量诱导速度场的修正，修正方法可按l j t h o m a s m d s a l a s 提出的基于小扰动方程的方法进行。 对位流计算，绕翼型的环量需按后缘k u t t a 条件确定，因此k u t t a 条件是 位流计算的唯一性条件。在n s 方程的计算中，只有满足k m t a 条件的解才是 n s 方程的解，在n s 方程的时间推进中，采用5 步m m g e k u t t a 推进，可以 满足k u t t a 条件。 3 2 3 紊流模型 三维机翼定常、非定常气动力n s 方程计算 工程上实用的计算湍流的方法是由不同湍流模型封闭的平均n s 方程 组，这里我们采用b a l d w i nl o m a x 代数湍流模型。 b a l d w i nl o m a ) 【代数模型在c e b i c i 代数模型的基础上进行了改进，避免 了求附面层外边界的必要性，将物面粘性层分成内外两层，分别用不同的代数 式确定紊流粘性系数以，即： ，f 0 ，) 。，y y 雎2 1 0 ，) ：，y y 。 y 为物面向外的法向距离， ( 3 3 3 ) y = m i n l y 。i ( ，) 。= ( ，) 。，d fy = y 。)( 3 3 4 ) 3 2 3 1 内层的紊流粘性系数 在内层中，应用p r a n d t l e 一、d r i e s t 公式 ( ，) 。= 2l 国i ，= 纠l e x p ( 一y + 爿+ ) 】 其中l l 为涡量大小 i 珊卜l “，一叱i y + ：型：巫 产m卢。 几和儿为壁面处气体密度与粘性系数，f 。是壁面摩擦应力。 3 2 3 2 外层的紊流粘性系数 ( 以) 一= 圈毛胪乙“( y ) ( 3 3 5 ) ( 3 3 6 ) ( 3 3 7 ) k 为c l a u s e r 常数，c 。为附加的常数。 f 。= m i n y 。，o 。，c 。y 。u f 。) ( 3 3 8 ) _ y 一为f ( y ) 取得时的y 值 南京航空航天大学硕士论文 ，( y ) ：y j j l c x p ( 叫+ 爿+ ) 】 ( 3 3 9 ) 在尾迹区中，上式中的幂项置零。( y ) 为k l e b a u o f f 间歇因子 讪) = 1 + s 文警，6 _ l b 4 u 。为速度型中最大总速度和最小总速度的差 【，。：f 了。一万了。 ( 3 4 1 ) 在固壁处【o 中的第二项取零，外层公式也可以用于尾迹中。 3 2 3 3 转捩的影响 如果某个固定的x 位置有0 ，) 。脚 c 卅。，则设这个速度型中 以= 0 。 3 2 3 4 常数的定义 上面公式中出现过的常数，定义如下： 彳+ = 2 6 ；c c p = 1 6 ；c m = o 3 ；c “2 o 2 5 ； 足= o 4 ；k = o 0 1 6 8 ；p r = o 7 2 ；p l = 0 9 ；c 。埘= 1 4 ( 3 4 2 ) 3 2 4 显式时间推进 在求解n s 方程的过程中，目前常用五步r k 方法推进，即 w ( 0 ) = w “ ( 1 ) ：w ( o ) 一口l 譬【q ( w ( 。) 一风】 俨) - w ( o ) 喝等【q ( ”) _ d l 】 w t ，) ；w ( m 一口，竽 q ( w ( z ) 一d ：】 三维机翼定常、非定常气动力n s 方程计算 其中 w ( 4 ) ：w ( 0 ) 一口。竽【q ( w ) d 3 n w ( 5 ) = w ( 一口5 竽 q ( w ( 4 ) d 4 九 w ”+ 1 = w ( 5 ) d 。= d l = d ( w o ) d 2 = d 3 = ( w 2 ) + ( 1 一户) d o d = 归( w 4 ) + ( 1 一，) d 2 p = 0 5 6 ，y 2 0 4 4 ( 3 4 3 ) 口i = 1 4 ，口2 = 1 6 ，口3 = 3 8 ，口4 = 1 2 ，口5 = l ( 3 4 4 ) 物理粘性仅在第一步上计算，并将其值冻结于以后计算步中，这并不太影 响整个计算格式的稳定性，却可节省大量的计算时间。而人工粘性只是在第一、 三、五步计算耗散项，其余各步冻结，每一步耗散前的系数作了调整。这样， 当真实粘性项占优时，可更好的保持数值稳定。 3 3 加速收敛措施 3 3 1 当地时间步长 在求解方程的过程中，每个网格允许的时间推进步长都不同，如果取统 一的时间步长，只能取他们的最小值，以确保流场中每个网格都满足稳定性的 c f l 条件，这样