（流体力学专业论文）低马赫数、低雷诺数翼型非定常流数值模拟.pdf

西北丁业大学砸士论文 摘要 本文以小型昆虫翼面运动的模拟为出发点，作为研究的第一步，将翼面运动的 研究简化成翼型的沉浮和俯仰两种运动的耦合，然后利用n s 方程去求解这种扑动 翼型的流场。由于在流场求解的过程当中，随着翼型的扑动，每一步都需要生成新 的网格，为了提高网格生成速度，我们采用效率较高的代数插值方法。另外我们 求解的是极低马赫数下的流场( m 0 0 5 ) ，极低的马赫数会导致可压流动方程十分 僵硬而难以收敛，对此，我们对n s 方程进行了矩阵预处理。在方程求解方面，采 用了有限体积方法，采用r u n g e - k u t t a 方法对方程进行双时间推进。为了加快解的 收敛速度，采用了当地时间步长以及隐式残值光顺等技术。本文对在低马赫数下 绕n a c a 0 0 1 2 翼型和n a c a 4 4 1 2 翼型的定常流场进行了求解，结果和实验基本吻 合，并对非定常运动情况进行了初步模拟研究，得出了一些有意义的结果。 2 两- l z 业大学顺 论文 a b s t r a c t t h i sp a p e rs t a r t sf r o mt h em o t i o no fi n s e c t s w es i m p l i f yt h em o t i o no f i n s e c t si n t o ak i n do fc o m b i n a t i o no f p l u n g i n ga n dp i t c h i n g ，a n dt h e nu s en se q u a t i o nt os o l v et h e f l o wa r o u n dt h ea i r f o i l i ne v e r yt i m e m a r c h i n gs t e p ，w em t l s tr e g e n e r a t et h em e s h i n o r d e rt o i m p r o v et h ee f f i c i e n c y o ft h em e s hg e n e r a t i o n ，w ee m p l o yt h et r a n s f i n i t e i n t e r p o l a t i o nm e t h o d ，i nl o ws p e e df l o w ( m 1 0 ) ；( 2 ) 诱导频率小；( 3 ) 雷诺数高( 雷诺数大干1 0 0 0 ，0 0 0 ) ；相反，昆虫的扑翼中，一般来讲，展弦比小于5 ， 诱导频率大于o 2 ，雷诺数小于3 0 ，0 0 0 。而且，在翼型剖面结构和翼表面状况之间 的差别也是很明显的。在飞行状态下，扑翼和旋转翼间亦存在重大的差别：( 1 ) 拍 动翼飞行周期中，其前冲速度是正弦导数而非恒定。( 2 ) 做顺桨运动以得到升力， 它主要是由在动力行程和回复行程的末端而不是旋转翼剖面的端，作用在机翼上 端的吸力。( 3 ) 顺浆铰链有大的角旋转，分别在半周期，内转，仰转的反向终端。 lu葛。弓辱8u鲁ueg 8 u 互兰山 一d 11c岩ug苎u o 世o i 色巴言z 西北工业大学碗士论文 因此，旋转翼气动力的定量结果不能直接应用。对机翼旋转在雷诺数为2 ，4 0 0 _ 卜的非定常气动性能效应的研究，证实了：在一个很高诱导频率( k - 2 5 ) 的中心涡 流线的影响下，在旋转阶段产生了一个高的升力系数( c 1 = 2 7 ) 。如前所述，薄翼 在头部总有一个分离泡。当分离泡部分从翼型脱落时，发生完全失速。涡进一步引 起尾缘附近反向旋转涡的脱落。 在拍动飞行中，在同样的速度和迎角下，其翼产生超过定常运动的升力。特别 强调的是：扑翼高升力的产生是流动分离或分离泡或尾缘及其附近的泡的效应，尤 其是由仰转过程的俯仰旋转而引起，这被作为昆虫飞行的一个新的机理。对昆虫飞 行来讲，困其机翼不是流线型的而是在头部为波状的或角状的，任何昆虫总是依赖 部分分开的机翼飞行。特别提到的是：在悬停飞行中，因为低的机翼速度而需高的 升力系数。此时，昆虫有时采用气动失速来产生大的平均升力。 若是在拍动的动力行程或下行程的终端附近，以高的前冲速度飞行，则能够产 生大的升力系数。在许多文献中，如前所述，强调绕叶片元件有一大的环量r ，它 主要由绕头部的快速流动引起。然而，若平移速度u ( 在气动中心或靠近c 4 处) 很小，例如内转中，则气动力会很小，因为升力l = p u f 。 对低雷诺数，低展弦比机翼，与其它机翼相比，昆虫的阻力系数是如此之大， 以至于最大升阻比发生在大升力系数处。在小昆虫的扑翼飞行中，在机翼转动的过 程中产生较大的升力，因此在动态滞止引起的最大升力系数增高的下行程过程中， 以高速产生大升力。 3 表面粗糙度的影响 国内外初步研究表明，一般情况下，表面粗糙度降低，升力线斜率、最大升力 系数、升阻比都变差。 4 高升力机理研究 在下拍过程中，空气绕头部形成漩涡，并形成一个强烈的前缘涡。对于弦线雷 诺数( 小于5 ，0 0 0 ) 尖头薄翼翼型，如果它在大迎角下运动，就会出现层流前缘涡。 前缘涡的环量增强了附着涡因此增加了升力。对天蛾运动的观察发现，在前缘涡里 西北t 业犬学顺七论丈 有一很强的沿着展向流动，当这个流动与前缘漩涡运动结合时，这就形成了螺旋形 的前缘涡。展向流动使涡向翼梢移动，在那里它与翼尖涡联合，从而防止了前缘涡 的脱落。因此，展向流动稳定了前缘涡，对整个下拍过程而言，延长了动态滞止的 有利效应，并产生足够的升力以支撑重量。对于昆虫来讲，沿着展向存在着拍动速 度梯度，从而引起了动压梯度，这种动压梯度就形成了展向的流动。 前缘涡周期性的增长和从机翼脱落( 在二维流动中) ，限制了平均升力系数。 展向流动及所引起的螺旋形的前缘涡，在层流雷诺数下是明显的现象，但当过渡到 湍流时则消失。这并不奇怪，因为即使是在层流的情况下，前缘涡也最终经历螺旋 形涡中心的分离和沿展向升力的减少。高雷诺数下湍流的引入可能破坏涡的完整， 并将之缩减到一般的二维非定常流动。 5 扑翼的一个完整的拍动过程 在下拍的多数过程中，有一强烈的圆锥螺旋前缘涡，中心有一很强的展向流动。 这导致上表面负压力范围的稳定增加，从而升力增加。 在下拍的后半部分：翼梢出现另一个前缘涡，中心也有螺旋流动，但是这个 涡朝机翼下表面移动。 在接下来的仰转中：这两个涡逐渐改变其形状。随着机翼到达机翼拍动的底部， 它们与一翼梢涡( 在仰转开始时出现) 结合，形成一钩状涡。 在仰转的后半部分：涡强烈变形，迅速瓦解，超出两个前缘涡在展向中点的结 合处，从尾缘脱落。 在上拍刚开始：不出现前缘涡上下表面的流动几乎光滑。 上拍的后半部分：出现前缘涡。 内转过程中：上拍时出现的前缘涡还在，但不在上表面。下表面下有一尾涡， 是内转上半部分机翼转动的结果。另外，上表面有一剪切层涡。 当机翼开始下拍的移动运动时，内转过程中观察到的涡都脱落因此机翼必须 产生一新的前缘涡以产生升力。实际上，也确实产生了前缘涡。 8 0 的升力产生在下拍过程和上拍过程的后半部分。 r 西北工业大学l i i - 论立 1 3 影响翼面涡产生和脱离的几个因素 正如f ； 所述，在低雷诺数的扑翼研究当中，涡是产生升力和推力的主要原因。 以前的研究主要是运用实验或者是数值分析的方法分析涡的产生形态，比如：日本 的早稻罔大学和法国的博伊特大学合作研究的“环绕大幅度振荡机翼的涡形态研 究”以及纽约大学的时h 翼飞行中的涡脱离形态以及频率的选择”等，他们都着重研究 了环绕扑翼的涡的形成和脱离。而且，他们一致认为：涡产生和分离的形态是和以 下几个因素密切相关的，也就是说升力和推力的产生形态和他们相关，那就是频率、 振幅、雷诺数、俯仰转动轴的位置、翼型的截面形状等 1 3 1 频率的影响 在较低的频率下，翼型横截面和俯仰转动轴的位置不会对涡的形成产生太大 的影响，首先会在前缘上表面形成一个较大的涡，这个涡在上表面的发展形成的初 始的非定常的尾迹，随后这个涡开始脱离，这使得上表面的流线呈反”s ”形，接下 来较小的后缘涡开始脱离。这个发展模式和不发生拍动的机翼几乎是相同的。一般 前缘涡是反时针的，而后缘涡是顺时针的。如果俯仰转动轴靠近前缘的话，那么前 缘涡的分离会更加平坦，而后缘涡的分离会更加稳定。在中等的频率上，后缘会随 着振荡交替的产生顺时针和反时针的涡，而且前缘的下转会加速前缘涡的发展，相 反，前缘的上转会抑制前缘涡的发展。在这种频率下，横截面和转动轴位置对涡的 影响更加明显了，在该频率下，前缘涡和后缘涡的发展是互不干扰的，各自独立地 发展，这被称为平行脱离模式。如果频率在继续增加的话，涡的发展模式和中等频 率相同，只是前缘和后缘的涡特征更加明显了，而且在这种情况下，转动轴和横截 面形状对涡的影响就更加明显了。如果转动轴靠近后缘的话，周期性产生的反时针 前缘涡就会一个个的在上表面重叠起来，最后形成一个单一的直立涡这种被称为 涡重叠模式。如果转动轴靠近前缘的话，周期性的前缘涡不会进入重叠区域，只是 9 西北| 二业大学坝l 。论史 有规律的向下游移动，和共同旋转的后缘涡( 反方向的后缘涡) 组台到一起，然后 有规律的向下游移动。最后以一排交替涡( 前、后缘涡) 的形式向下游脱离。这种 被称为同步脱离模式。 1 3 2 振幅的影响 现在谈谈振幅，总的来说振幅会带来两个明显的效果：第一，前缘涡和后 缘涡都会加强；第二，下游流尾迹会加宽。有时随着频率的增加，振幅的变化还会 带来一些其它的作用。在低频率的情况下，随着振幅的增加，前缘涡会变大，并且 在分离的过程当中变得更加稳定。如果我们增加活动范围( 初始迎角加上振荡迎 角) ，我们会观察到同样的现象。通过一系列的实验表明，在低频率情况下，对涡 发展影响最大的不是振幅，而是活动范围。轴位置的影响是很小的。 在中等频率下，频率和振幅起着相互抵消的作用。在转动轴靠近前缘的时候， 涡的发展模式主要由平衡位置的的值决定，而不是振幅。当转动轴靠近后缘的时候， 涡的发展模式会随着振幅的增加而逐渐变化。 在高频率的情况下，尾迹的模式在很大程度上取决于振幅。在较大的振幅下， 涡重叠模式和同步脱离模式交替出现；在较小的振幅下，它们是平行脱离模式。 1 3 3 雷诺数的影响 当雷诺数逐步增加的时候，初始尾迹的发展会被加速，从而减小了尾迹从初始 状态到稳定状态的过渡期。特别是当雷诺数达到1 ，0 0 0 0 时尾迹里面会出现湍流 影响，基本上来说，尾迹还是均匀的。同时，雷诺数对涡的影响是和其它参数相联 系的，也可以说，其它参数在某种程度上决定着雷诺数的作用。而且，和其它参数 相比较，雷诺数的影响是比较小的。 当振幅和振动范围以及频率都比较小的时候雷诺数的影响是可以明显看见的。 雷诺数的增加会导致前缘分离流随着振动迅速的转换。 o 西北工业 学颁j 论文 当频率增大时，雷诺数的一个明显的影响就是会促进尾涡的发展以及下游尾迹 的波动。在这个雷诺数下，涡的模式是介于平行脱离模式和同步脱离模式之间的， 但在较小的平衡位置下，涡的脱离模式趋向于同步脱离模式，这可能是因为雷诺数 的增加促进了后缘涡的发展的缘故。 最后，当频率进一步增加时，雷诺数的作用取决于振幅，当振幅稍大的时候， 涡的发展模式是同步脱离模式和重叠模式的交替出现。同步脱离模式对雷诺数的变 化不是十分的敏感。相反，重叠模式对雷诺数的变化十分敏感。特别是在雷诺数为 l ，0 0 0 的时候，无论是顺时针的还是反时针的前缘涡，它们很少是紊流，它们的发 展也十分的清晰。而后缘涡则具有更大的能量而且很稳定，它会以一种交替的方式 影响下游的流动。在较小的振幅下，雷诺数的增加会加强活动范围和尾迹模式之间 的关系。 1 3 4 转动轴的影响 计算和实验显示，在较高的频率下和较低的振幅下，转动轴的影响会更加明显。 当把转动轴向前缘移动的时候，这会促使尾迹从平行脱离模式向同步脱离模式转 变。这是因为当我们移动转动轴的时候，这会改变前缘涡和后缘涡的转动半径。 如果机翼在靠近后缘的地方振荡，后缘涡会变得比较大。同时，涡交替脱离的规律 会更加明显此外，由于机翼下降而产生的诱导吸力的作用，下游流动对前缘涡束 缚会变小。这两种作用合并起来，就会促使涡的发展模式向同步脱离模式转变。 如果我们采用较高的频率和较高的振幅，转动轴的位置就直接决定了尾迹的模 式是涡重叠模式还是同步分离模式。这就在很大程度上意味着，涡重叠模式的产生 不仅仅是由频率和振幅决定的，而且也是转动速度和前缘旋转半径作用的结果。因 此我们可以推导出：在较高的频率情况下，转动轴的作用可以补偿转动频率的作用。 关于这一点，1 9 8 5 年，哈林和沃克阐述了，也就是，将转动轴向后缘移动会起到和 西北t 业火学颂l ：论文 增加转动速度相同的结果。很明显，这一点和以上的阐述是一致的。 1 3 5 横截面的影响 总的浣来，横截面所起的作用和转动轴所起的作用很相似，现在可以用 n a c a 0 0 1 2 翼型和椭圆形翼型做比较说明。 n a c a 0 0 1 2 翼型具有比椭圆形翼型优越的一点，那就是它可以在圆形的前缘产生 灵活的分离点和方向，就因为这点，使得n a c a 0 0 1 2 前缘点的流体随着振动的分 离更加有规律性。但是对于椭圆形翼型而言，流动的分离常常在振动的顶点有着突 然的加速。因此n a c a 0 0 1 2 翼型的前缘涡一般比较大而且尾迹的发展会得到加速。 此外，在较高的频率下n a c a 0 0 1 2 翼型的较小的前缘曲率带来三个新的效果：l 、当 表面下降的时候，吸力的作用会减少2 、表面上升的时候，顺时针的涡会受到压缩 3 、上表面和下表面的涡流的作用会减少 另一方面，由于n a c a 0 0 1 2 翼型具有尖后缘，后缘流动分离会进一步发展，而且 会更加的稳定，这可能是因为有涡的下表面流动在切向上的分离作用加强，因此后 缘涡比较大而且很稳定。所以，下游的涡尾迹呈现出有规律的交替模式。此外，稳 定的尾涡会沿着上表面产生更多的逆流的能量，这会促使前缘的反时针涡产生分 离。前缘涡在下游的运动也被促进了。 综上所述，和椭圆形翼型比较而言，n a c a 0 0 1 2 翼型有以下三个特点：l 、前缘 涡在下游的运动被加速：2 、后缘涡被进一步发展而且更加稳定；3 、下游流动涡的 交替脱离的规律性增强。这些现象在高频率下，中等振幅下尤为明显。 1 4 本文工作 针对微小型飞机工作在低马赫数、低雷诺数环境的现实，而且由于扑翼的非定 常效应将是微型飞行器发展的必然趋势，作为研究工作的第一步，本文重点研究翼 型在低马赫数、低雷诺数下定常、非定常流场数值模拟方法。主要工作包括以下两 西北工业大学硕士论文 方面： 1 ) 对可压缩的n s 方程进行顸处理，使它能够模拟低马赫数流场。 2 ) 由于在大幅度的俯仰运动情况下，可能会导致翼型尾部网格正交性恶化，因 此要改进0 型网格，主要是通过将一定层数的网格和内边界固连在一起运动 来实现 些苎三兰竺查兰坐二! 堕苎 第三章n - s 方程的求解 3 1n s 方程 我们选择一个随着流体一起运动的微团，然后在这一微团上利用牛顿第二运 动定律，我们就可以得到微分形式的二维n s 方程： 掣+ 掣+ 掣：o【3 1 】 o tc ! xd r 其中： u = ( p ，p u ，p v ，p e ) 7 e = ( p u ，p u 2 + p ，p u v ，p u e ) 7 f = ( p v ，p u v ，p v 2 + p ，p v e ) 7 k e 。= t o ，r 。，t x t ，“f 。+ v r w g ，j 。： 3 2 】 c = 【o ，k ，峨。+ v g y v 一以y ； 铲嚣c z 塞一 铲署c z 考一羚 铲可i t 百。| o u + c 。j ua t 吼一赢面； c 。a t n 一赢面； 现在，我们在空问取一个控制面q ，它的边界线和边界的外法向分别为a q 和； 对控制面内的流体微团随流体跟踪一段时间。当6 t 寸0 时，我们就可以得到n s 方程的积分形式： 导p + 扣：胪。 b ， 昙p u d s q 十( 孑：) 出= ( 一p + k ) ( e ：) 凼+ 卜，( i ：胁 一nm 昙擒十c 赫出2 p 训狮出+ 黔赫出 鲁p 地+ p 赫出= - _ - - + 一 j ( “f “+ v r 廿一g ，) ( n ) d s + j ( “r ”+ v r 卅一叮，) ( f ，n ) d s * 1 以上是四个方程t 但是有五个未知数，那就是p ，p ，“，v ，e ，所以我们还得给出一个 方程，那就是： p ：( z 1 ) p ( e 一半2 )【3 4 】 这样就可以将这个方程组解出。另外我们还可以统一的把上式写为 鲁啪+ 驴：西= 。 b s ， 其中 w = p 删 删 , o e f = p g p u g + ( p f 。) t t 。i ， 斗 i - g + ( p f w ) f ，一i 。 印目+ ( u r 嚣+ v f 掣一g ，) f ，+ ( u r 叫+ v f 一q y ) f y 3 2 矩阵预处理 在低马赫数下，由于速度和音速相比之下变得非常小，方程的几个特征值之 间的差异会变得很大，这样就会使得方程变得很僵硬时间推进方案就会收敛的很 慢，当雷诺数变小的时候，这种情况会更加严重。这种情况一般发生在包含低速区 域的可压缩流场以及极低马赫数的不可压缩流场。因此，必须采取措施减轻方程的 僵硬程度。最早采用预处理矩阵解决这个问题的是c h a r l e sl m e r k l e 他先后发展了 一 一 堕j ! ! 些查竺些主堕壅 无粘的预处理矩阵，粘性预处理矩阵以及非定常预处理矩阵。在这晕，我们采用的 是由p l e t c h e r 和c h e n 发展的非定常粘性预处理矩阵。 没有预处理的非定常n s 方程如下： 型+ 丝+ 竺：堡+ 堡 a。x a y 8 x a y 引入双时间推进方案以后，方程变为： 翌竺q + 旦竺q + c f u l ：0 d fa t 其中，c ( u 1 是通量项，现在再引入一组变量： w = p ，“，v ，r 7 这样方程可变为： r a w q + a c _ 2 q + c f 1 ：o d rd f 其中，t = o u a ，在这里，我们把对时间导数的预处理加上，方程就变为： r d w f 2 + 型q + c ( 卜o d fd f p e t c h e r 和c h e r t 发展的矩阵： f = 一1 00一旦 j 。7 旦p0 一坐 yy 。 三0p一旦 c ，+ 等寺f ) 3 6 】 3 3 空间和时间离散 1 ) 空间离散 在空间离散这部分当中我们用三维的空间离散来作介绍。用贴体曲线网格 娈苎：三些查兰丝土堡苎 来代管绕流物体( 如机翼) 周围的连续流动区域，对其中的任意一个六面体单元， 上式成为 鲁。j 睁+ 。拈。 b ， 其中边界a q 。是这样构成的： 锄讲一书t 哆，t + t ，争+ t 。t 。e 一；+ t 础畦刚 让旁去单元中心的平均值，即 矿”5 去。肜棚 刚公式r 7 、可以写成： 妄( 托。戚m ) 厩舻o 3 9 】 上式表示流动的通量平衡，单元体体积 。可以用三个面矢量以及其主对角线 r 。肚进行计算，如下式： 1 呻_+叶 自z 。t2 专7 r ， 。( s z 一；m t + s , , j q , + sr ， 一；) 一 q 。m 表示流出单元体的净通量： q r - ，* 。f s i + ，一f s l ，一；。+ f s l 3 1 0 】 1 。一f s 1 j 。+ f s1 1 一f s l 1 +il-ji lr j i j , 2 k 一= 上式右端各通量项的计算都是类似的，例如 f s i l= f - ，一= - p ，h ( 可一y 乩一 ( 刚。j ：s ( q 。加u 一争m 斗 - + ( 川。一 峨一抄p 。一辩 ( 删+ u 一扩s 一抄p 。广鼻 ( ) 。一小一j i k 。s2 - 3 1 1 】 西北工业大学顺二e 论文 所有带分数下表的物理量的计算均取相邻量的平均值，例如 f ，r h = y 1 w r ，女+ w u t ) 由于上式使用了中心平均公式，因此上述空间离散格式对均匀直网格具有二阶 精度。对于一般曲线网格，格式的精度取决于网格的光滑程度，如果网格足够光滑， 即网格的形状和大小都不发生突然变化，则格式的精度近似是二阶的。 2 ) 时间离散 我们对实时间导数项采用的是隐式差分，也就是 d u d t ! 竺：二! 竺：竺二 31212at 这种隐式格式具有二阶时间精度 3 4 人工粘性 使用中心平均公式的有限体积格式不具备耗散性，这就是说，某些高频误差分 量在求解过程中是不衰减的，造成这种情况的原因是用差分代替微分时，舍弃了高 阶导数项。如果在有限体积格式中适当的加入一些高阶导数项，就可以避免出现振 荡现象，使数值解完全的收敛到定常状态，以得到光糈的压力分布，这些人为加入 的高阶项，称为人工耗散。 这罩，由于加入了矩阵预处理技术，所以对人工粘性的求法也要加以一定的修 改，因为矩阵预处理的加入改变了方程的特征值，必须将这部分影响加入到人工粘 性当中去。引入人工粘性后，方程变成了如下形式： r 孚q + 掣n + c ( c ，) 一d ( ) ：0 3 1 3 a fd f 在二维情况下，我们可以把人工耗散写为如下形式： d ( 畎，t ) 2 ，- d 州。 - d 。一； 西北工业大学倾士论文 在这个公式当中，对于边( “= 1 ，) 的人工耗散通量为 r ，al ”( 嘭扎，一形，) + i ，7h 亍h i ，f 5 墨( 彬+ ，一3 形扎，+ 3 彬，一彬。，) 1 2 7 为了保证c 一d 完全的守恒，所以我们应该把预处理矩阵r ，引入到人工耗散 h 7 7 当中去。 基于离散化方程最大特征值而导出的系数变量的表达式如下： a 一。= 圭 ( 矾，+ ( 砜u 】 人。+ ；= 圭【( 矾+ ( 砜一 在这里，f 指的是i 方向，叩指的是j 方向，并且有： 石= t m ；石= 厶中。 还有： 咋c ；卟c 争 对于一般的未经预处理的系统来说，它的最大特征值为： = u + c s 对于经过预处理之后的系统来说，它的特征值应为： =r + 1 ) u + 拓i 币再瓦而f 】 3 1 4 】 其中： u = p d 司：s = i d s i 预处理矩阵系数的表达式如下： r ，= 昙 r 。+ f 。， j + = ，z 耗散系数和压敏系数的表达式如下： 西北工业大学埘士论空 ，= 。犯1 m a x ( _ m ，- ) 占 = m a x o ，( 女h 1 一s 旺? ) 】 ”i 。”一 。：区! 二堡! 坠! “7 p 山+ 2 + p ，一i ， 一般情况下：我们取占= 必，k 2 1 = 女4 1 = 1 名8 4 3 5 时间推进方案 积分形式的n s 方程经过了对实时间以及空间的离散并且为了保持方程的守 恒，我们引入了预处理的人工耗散，这样方程就可以写为如下形式： 磊d 畔，+ 只，= o 其中： 只，= 寺 坐卷型坝卟冽 上式p 一就是我们通常所说的残值。这个方程是一个关于时间的常微分方程， 可以采用r u n g e - k u t t a 显式时问多步法加以求解。 现在假设枷1 是蹿次虚时间推进的结果，对于上述方程的经典四阶r u n g e k u t t a 方法可以写成如下形式： 其中 p p p p p 缸一：肿一：打一：缸一： 一 一 一 一 旷 肜 旷矿 f f f f f f f i i i i i 即帅 帅 帅 帅 西北 4 业大学顺七咤义 第四章边界条件 n - s 方程是所有粘性流动所满足的基本方程，只有在一定的边界条件和初始 条件下，方程的解才是唯一的。因此数值求解n s 方程的重要问题之一就是对边界 条件的离散。如果边界条件处理的不好，不仅会给计算带来较大的误差，也可日z - 。i ! 使 计算过程不稳定，甚至迭代不收敛。数值求解绕物体的定常流动时，通常要对物面边 界、远场边界、周期边界性条件和对称边界条件进行边界条件处理。对于周期性边 界和对称边界处理比较简单，本节主要介绍物面和远场边界条件的离散方法。 1 ) 物面边界条件 粘性流流过不可穿透的固壁时需要满足无滑移条件，也就是在壁面上有： “。= 0 ；v 。= 0【4 1 】 当我们采用中心差分的有限体积法求解n s 方程时，所有的变量都应该是格心 点的值，当我们赋贴近物面的第一层格心点的值的时候，可以用贴近物面的那几层 进行插值，也可以直接把物面上的值赋给它。 2 ) 远场边界条件 在绕物体的流动当中物体所产生的扰动可以被认为是传播到无穷远处而没 有被反射，在物理空间进行有限体积离散的时候总要取一个有限远的远场作为 边界。以均匀流场为初场作时间推进的时候，可以认为有扰动波从物面逐层网 格的向远场传播，远场边界条件的处理必须能够不反射这些扰动波，否则时间 推进就有可能不收敛。本节将通过局部准一维特征分析，给出远场无反射特征 边界条件。 在远场边界点o k 建立局部直角坐标系o # r g ，使f 轴与远场法向重合，掌， v 轴在切平面内。沿f ，吁，善方向的速度分别是q l ，q ：，q ，因此远场的法向速 度q 。= 吼。 西北工业火学硕士论文 欧拉方程的微分守恒形式可以写成 型+ 竺+ 一o g + 一o h ：o l 8 0 8 q 其中： g = p 删 印 p w o e 脚2 脚1 q 2 用；+ 尸 p q 2 q 3 p h q 2 f = p q l 期j + p p 孔9 2 p g 【吼 p h q l p q 3 ，卿1 q s p ：( y - i ) p e 一= 1 ( “2 + v 2 + w 3 ) 】 肚n 形 令西= ( n p c l ；，p q ：，p q ，p ) 7 ，上标表示矩阵转置。 现在以u 为参数，则欧拉方程可以写为： 型+ 一型+ b 型+ co _ y u ：0 o t 0 目叩d 奠中： a uo f a a 已ta u b ：型一o g a ua u 3 5 c ：型一o h a ue u 4 2 鲫麻鹏 西北工业丈学倾上论文 c = q 3 0 0 q 3 00 00 0o 0 p 0 o0o 吼0 0 0 吼一1 p 0 r p 吼 因此，c 的五个特征值为： 其中 a = 形 口为音速，q ，是远场边界的外法线方向的速度。 【4 3 定义欧拉方程沿喜轴的准一维特征线是( 手，r ) 平面上的一条曲线工，在这条曲线 上，物理量u 本身的值是给定的而物理量的法向导数沿j 【是不确定的。因此可以 写出对应于五个特征值的特征线方程为： 在( 亭，f ) 平面上，欧拉方程可以化为 o _ q _ v + c 型：o 新 a 毒 也就是： 4 5 4 6 】 44 d 口 一 十 g g g g g i j = = j j i j 如 办 仃 d 口 g g 宁 q = = = = = 蝣一面鸳一研蝣一毋蟛一出蟛一出 西北丁业大学硕士设立 害怕詈+ p 罄= o 鲁怕罄= o西”a 亡 堕+g婴：oat 1 a 芒 鲁怕詈+ 古嚣= 。 竺+ 矿饕+ q 3 罢：o c 3 t 。 a fa f 由上述的第一和第五个方程可以导出： 瓦a s 嘉= o研”a f 其中： s 2 形，它是流体的熵。 因此，对等熵流动，上述方程可以化为： 詈+ 吼嚣+ p 等= o 盟+ 日、塑：0 8 t 。鸽 堕+ 口堕：0 研 a 百a q 3 怕嚣+ 吾嚣一。 鱼c 3 t 怕嘉= o”a f 将上述的第一个方程乘以兰再与第四个式子相加，可以得到 p c 昙心剖势。鲁 因此可得到，方程沿特征线的r i e n a m m 不变量分别是： 【4 7 】 4 8 】 4 9 】 两北工业大学 i ! ；i 十论文 r t 2 q 】 r 2 = q 2 r 3 = 5 r 。：吼一兰 y i 弘”鲁 4 1 0 】 数值计算时，远场边界相对物体较远，可以认为物体产生的扰动沿远场的法向 传播。因此，只要特征线物理量满足上式的r i e n a m m 不变量关系，可以认为物体 的扰动波在远场无反射。根据特征值的负号关系可以导出无反射特征边界条件如 下： 1 亚音速流x ( - a q 。 o ) ，如，也， 都小于零，也大于零，因此，r i , r ：，r ，凡取来流值，而r ，由内场 外插。即在远场边界上的物理量满足关系 其中，上标”m “表示来流值，e 。表示内场外插值。 2 亚音速流出( 0 q ) ，如， 都大于零， 小于零，因