（流体力学专业论文）嵌边法出流条件在可压缩流直接数值模拟中的应用.pdf

中文摘要 直接数值模拟是目前研究流动问题的有效方法。由于计算域总是小于实际流 动区域，选择恰当的边界条件是保证计算可靠的必要条件。对可压缩流，较常用 的是基于特征关系的边界条件。但在非定常流的计算中，这种方法有且寸并不能给 出令人满意的结果，特别是要求无反射的出流边界条件。其主要原因是，从特征 关系角度看，只要是亚音速区，就总有入射特征波。除非能事先给出此波之值， 否则总会有误差。而事先给出其值一般来说是不可能的。 在实现无反射出流条件的方法中，有一种所谓的嵌边法( f r i n g e m e t h o d ) ， 也称缓冲区法、海绵区法等，其主要思想是：在出流处增加一段特殊区域嵌 边区，在该区域内，经过修正的n a v i e r s t o k e s 方程可将出流扰动衰减至足够小以 尽可能地减弱出流反射波的强度。该方法对不可压缩流的有效性己在实际应用中 被证实。本文采用直接数值模拟的方法，以嵌边法为出口边界条件试算了亚、超 音速平板边界层扰动波的演化、超音速流中涡的传播、亚音速平板边界层扰动在 尾缘的绕射问题，证实了嵌边法在可压缩流的计算中，比一般基于特征波分析的 方法，可以取得更接近于真正无反射边界条件的结果。特别是在亚音速情况更是 如此。本文还提出了一种衡量边界条件是否有效的评判标准。最后讨论了嵌边法 的参数选取问题。 关键词：直接数值模拟，嵌边法，边界条件，可压缩流 a b s t r a c t d i r e c tn u m e r i c a ls i m u l a t i o ni sa ne f f e c t i v em e t h o dt os t u d yf l u i df l o w s b e c a u s e t h ec o m p u t a t i o n a ld o m a i ni sa l w a y ss m a l l e rt h a nt h ed o m a i no fr e a lf l o w s ，a p p r o p r i a t e b o u n d a r yc o n d i t i o n ss h o u l db ea p p l i e dt oe n s s u r et h ec o r r e c m e s so ft h ec o m p u t a t i o n a l r e s u l t s f o rt h en u m e r i c a ls i m u l a t i o no f c o m p r e s s i b l ef l o w s c h a r a c t e r i s t i cb o u n d a r y c o n d i t i o n sa r eo f t e nb e i n gu s e dt oe n s u r et h a tt h ed i s t u r b a n c e sc a l lp a s st h r o u g ht h e b o u n d a r yw i t h o u tc a u s i n gr e f l e c t i o n s h o w e v e r ，s o m e t i m e s ，t h e yd on o ty i e l d s a t i s f a c t o r yr e s u l t sf o ru n s t e a d yf l o w s o nt h eo t h e rh a n d ，“f r i n g em e t h o d ”，a l s oc a l l e d b u f f e rz o n em e t h o d ”。h a sb e e n p r o v e dt ob ee f f e c t i v ei ng u a r a n t e en o n - r e f l e c t i o nf o ri n c o m p r e s s i b l ef l o w s i t sm a i n i d e ai st h a taf r i n g e ，o rb u f f e rz o n e ，i sa d d e dt ot h eo r i g i n a lc o m p u t a t i o n a ld o m a i n ，i n w h i c ht h eg o v e r n i n ge q u a t i o n sa r e m o d i f i e d ，s u c ht h a td i s t u r b a n c e sa r ea l m o s tt o t a l l y d a m p e di np r o p a g a t i n gt h r o u g ht h i sb u f f e rz o n e ，s ot h e r ew i l lb en od i s t u r b a n c et ob e r e f l e c t e da tt h ee x i to f t h eb u f f e rz o n e i nt h i sp a p e r ，b ya p p l y i n gt h e “f r i n g em e t h o d t os e v e r a lu n s t e a d yc o m p r e s s i b l ef l o w s ，s u c ha st h ee v o l u t i o no f d i s t u r b a n c e si n s u b s o n i ca n ds u p e r s o n i cb o u n d a r yl a y e r s ，v o r t e xc o n v e c t i o ni nas u p e r s o n i cf l o w ，a n d t h es c a t t e r i n go ft sw a v e sa tt h et r a i l i n ge d g eo fas u b s o n i cb o u n d a r y l a y e r so naf l a t p l a t e ，w es h o w e dt h a tt h ef n n g em e t h o di sa l s oe f f e c t i v ef o rc o m p r e s s i b l ef l o w s w e h a v ep r o p o s e dac r i t e r i o nf o rt h ee f f e c t i v e n e s so f t h em e t h o d ，a n da l s od i s c u s s e dt h e e f f e c t so f t h ep a r a m e t e r si nt h ef r i n g em e t h o d k e yw o r d s ：d i r e c tn u m e r i c a ls i m u l a t i o n ，f r i n g e m e t h o d ，b o u n d a r yc o n d i t i o n ， c o m p r e s s i b l ef l o w 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的 研究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表 或撰写过的研究成果，也不包含为获得盘生态鲎或其他教育机构的学位或证 书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均己在论文中 作了明确的说明并表示了谢意。 学位论文作者签名：影彩；移矿签字日期：弘西年6 月印日 学位论文版权使用授权书 本学位论文作者完全了解盘鲞盘鲎有关保留、使用学位论文的规定， 特授权鑫盗盘堂可以将学位论文的全部或部分内容编入有关数据库进行检 索，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校 向国家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名：絮剜砌导师签名 遑南 签字日期：埘年占月刁日签字日期：砷。作g 月2 7 日 第一章绪论 第一章绪论 数值模拟是研究流动问题的有效方法。当问题本身遵循的规律比较清楚，所 建立的数学模型比较准确并为实践证明能反映问题本质时，数值模拟具有很大的 优越性。一般来说它耗费少，时间短，条件易于控制，便于优化设计，比实验研 究更自由、更灵活。特别是随着计算机及计算方法的发展，这些优点将更明显。 钱学森先生曾认为，将来力学问题绝大部分都可以在弄清其物理本质并得到可靠 的数学方程后通过计算而解决，实验的作用主要是为弄清物理本质和验证计算结 果。对于湍流和转捩的研究也是如此。 做数值模拟首先要建立反映问题本质的数学模型，然后寻找高效率、高精度 的计算方法，其中包括微分方程的离散及求解方法，边界条件的处理等等。计算 方法的研究是当前数值模拟面临的重要问题之一，其中，边界条件的处理更是难 点之一。由于计算域总是小于实际流动区域，也就是说计算域外总有与计算域内 流动相互作用的流动存在，选择恰当的边界条件是保证计算可靠的必要条件。早 期，边界条件的重要性并未得到广泛的认识，而实际上，计算边界条件的给法不 仅影响数值稳定，还会影响到计算结果的精度。 1 1 边界条件概述 关于边界条件的先驱工作可以追溯到上世纪5 0 年代。k a w a g u t i ( 1 9 5 3 ) 试图 研究边界条件对计算结果的影响。他用台式计算机，每周工作2 0 小时，用了一 年到一年半的时间，才得到了雷诺数为4 0 的绕圆柱流的一个解。因而无法重复 这个计算以试验三种不同的流出边界条件0 1 。在电子计算机出现之前做这种计算 试验是不实际的。 c h e n g ( 1 9 7 0 ) 通过他对b u r g e r s 方程的数值实验指出，对于网格雷诺数大于 1 的情况，由于边界条件的处理，高阶格式可能还不如一阶格式精确，而且边界 误差可能是内点截断误差的二倍。他注意到，对于同样的圆柱绕流问题，j e n s o n 和h a m i e l e c 等“3 ( 1 9 6 9 ) 的计算之间存在着很大的矛盾，这种矛盾是由边界误 差引起的，甚至当缸寸0 时，这种误差仍然保持着。从某种意义上说，在计算 流体动力学中，边界条件被看成具有支配作用的重要性。1 。 出流边界条件是计算流体力学中最常要处理的边界条件之一。例如，在水波 问题中，水波可以传得很远。但研究水波问题的计算域却不能很大。如果处理不 好，计算域中的水波在穿越边界时就可能在边界产生反射波，使得在计算域中原 第一章绪论 本并不复杂的波系变得越来越复杂，而这显然与实际不符。因此，如何提边界条 件，使得计算时不出现反射波就是一个重要问题。在转捩及湍流的直接数值模拟 中，如何使计算域中的扰动在边界上无反射，也同样是保证计算可靠的前提。不 少人对于出流边界条件采用外推的方法。对于有些问题，外推方法可以给出相当 好的结果，特别是对定常流。但也常常会导致计算不稳定，或计算结果不准确。 f r o m m ( 1 9 6 3 ，1 9 6 4 ) ”。5 在对不可压缩非定常流体的计算时使用了周期的流入和流 出边界条件，其优点在数学上是显而易见的。但是它的适用情况受到很大的局限。 对于很多工程技术问题，往往要求在边界上不会产生多余的扰动，即所谓的 无反射。因此在2 0 世纪7 0 年代以后，很多学者开始致力于无反射边界条件的研 究工作。 1 2 无反射边界条件( n o n r e f l e c t i n gb o u n d a r yc o n d i t i o n s ) e n g q u i s t m a j d a ( 1 9 7 7 ) 最先提出了无反射边界条件的理论基础，并提出 了一种构造近似无反射边界条件的方法，该方法使得在边界上的某一方向的波被 吸收掉”l 。k w a k 成功地利用该方法针对小扰动方程计算了机翼附近的非定常的近 音速流”1 。b l a s c h a k 和k r i e g s m a n n 也做过类似的研究工作”1 。k r e i s s 分析了对 于双曲型系统初始边界值混合问题的适定性”1 。这表明，使用无反射边界条件是 有效的，它可以使计算域大大的减小。实际上，“无反射”并不是在边界完全没 有反射波传入计算域，而是指采用某种近似的方法，使得在边界上不产生显著改 变上游流场的反射波。这样，计算域不必取得很大，计算量可以大大减小。 迄今为止，对于无反射边界条件已经有了较多的研究和认识。r b e r t h e t o a s t r u c ( 2 0 0 3 ) “把无反射边界条件分为非局部的边界条件和局部的边界条 件。非局部的边界条件在边界上由某种积分关系式导出m 1 ，或者用物理量的傅立 叶变换来得到边界条件“”3 。局部边界条件包括辐射边界条件。”1 ( 基于波方程 解的渐进展开) ，基于特征关系的边界条件“”2 2 1 和嵌边法“”1 ( 缓冲区法) 。 1 2 1 非局部边界条件 早期关于边界条件的研究的一个方向是建立精确的边界条件( 至少在离散精 度的意义上) ，精确的边界条件满足外部无限域内的所有方程和物理边界条件， 包括无穷远辐射条件。因此，其基本思路是解析地研究无限域内的波动问题。然 后由解析解或其控制方程组导出边界条件。由此导出的边界条件在一般情况下以 非局部性为特征，称为非局部边界条件。它们可以在空间上非局部，即所有边界 点的运动相互耦合，也可以在时间上非局部，即边界点在某时刻的运动与过去所 第一章绪论 有时刻的运动藕联，也可以同时在空间和时间上具有非局部性。采用边界积分方 程导出非局部边界条件是较常用的方法，但它只适用于线性问题。 很多学者做了不少有关非局部边界条件的研究工作”。“1 。j ，b k e l l e r d a n g i v o li ( 1 9 8 9 ) 对于简化的波动方程和有限域中的l a p l a c e 方程使用有限元的方 法在人工边界上得到了一个精确的d t n ( d i r i c h t e tt on e u m a n n ) 边界条件1 。 m i c h a e lb g i l e s ( 1 9 9 0 ) 利用傅立叶变换寻找入射波，提出了一个通用的理论用 于构造欧拉方程的无反射边界条件“”。他构造的定常边界条件( s t e a d y s t a t e b o u n d a r yc o n d i t i o n s ) 在线性假设下，误差为二阶非线性项，非定常边界条件 是理想“无反射”的二阶近似，而准一维边界条件( q u a s i - o n e d i m e n s i o n a l b o u n d a r yc o n d i t i o n s ) 是一阶近似ig j i n m b r a z a ( 1 9 9 3 ) 把波动方程得 到的无反射边界条件( d i r i c h l e t 型边晃条件) 应用到了椭圆型非定常剪切流的 完整的n a v i e r s t o k e s 方程的数值模拟中o “。 在r b e r t h e t d a s t r u c “的文章中提到，非局部边界条件是精确的边界 条件，这里“精确”指的是导出的边界条件的提法是精确的，它精确地满足外部 无限域内的所有方程和物理边界条件，由于计算量的庞大求解时也必须采取某种 程度上的近似，实际上得到的仍是相对近似的结果。没有任何一种边界条件是精 确的、完美无缺的，除非计算域取得比实际流动区域还大，但这是不可能的。 由于非局部边界条件将所有边界点的运动耦合或者将全部过去时刻的运动 与现时刻的运动耦联，因此，就近场波动在时空域的数值模拟而言，即使是对于 线性问题，大计算量的非局部边界条件也难以应用。从6 0 年代开始，对于局部 的边界条件逐渐成为研究的焦点。 1 2 2 局部边界条件 局部边界条件具有“局部性”的特点，即一个边界点在某一时刻的运动仅与 其相邻点在相邻时刻的运动有关，这样可大大简化数值计算。局部边界条件尽管 不能严格的模拟无限域对近场波动的影响，但能够在满足实际需要的精度下模拟 这一影响。常用的方法有：辐射边界条件、基于特征关系的边界条件和嵌边法。 辐射边界条件通过对单向波方程的近似，将波方程展开，模拟外向波在边界 上的被吸收过程( 实际是传出过程) ，也称为吸收边界条件。随着近似的程度提 高，作为吸收边界条件的偏微分方程的阶数也要提高，而高阶边界条件的稳定性 不易证明，即使能够证明。将它用于数值计算时的差分格式也相当复杂。很多工 作致力于寻求形式简单而近似程度高的边界条件“”。但是，辐射边界条件存在 着很大的缺陷：它只适用于由双曲型方程控制的线性系统。 基于特征关系的边界条件在不可压缩和可压缩流中广为应用“”2 。它的基本 第一章绪论 有时刻的运动藕联，也可以同时在空间和时间t 具有非局部性。采用边界秘分方 程导出非局部边界条件是较常用的方法，但它只适用于线性问题。 很多学者做了不少有关非局部边界条件的研究工作”。”。j b k e l l e r d a n o i v o l i ( 1 9 8 9 ) 对于简化的波动方程和有限域中的l a p l a c e 方程使用有限元的方 法在人工边界上得到了一个精确的d t n ( d i r i c h l e tt on e u m m m ) 边界条件1 。 m i c h a e lb g i l e s ( 1 9 9 0 ) 利用傅立叶变换寻找入射波，提出了一个通用的理论用 于构造欧拽方程的无反射边界条件”。他构造的定常边界条件( s t e a d y s t a t e b o u n d a r yc o n d i 2 i o n s ) 在线性假设下，误差为二阶非线性项，非定常边界条件 是理想“无反射”的二阶近似，而准一维边界条件( q u a s i o n e d i m e n s i o n a l b o u n d a r yc o n d i t i o n s ) 是一阶近似。g j i n m b r a z a ( 1 9 9 3 ) 把波动方程得 到的无反射边界条件( d i r i c h l e t 型边界条件) 应用到了椭圆型非定常剪切流的 完整的n a y i e r - s t o k e s 方程的数值模拟中”。 在r ，b e r t h e t d a s t r u e “的文章中提到，非局部边界条件是精确的边界 条件，这里“精确”指的是导出的边界条件的提法是精确的，它精确地满足外部 无限域内的所有方程和物理边界条件，由于计算量的庞大求解时也必须采取某种 程度上的近似，实际上得到的仍是相对近似的结果。没有任何一种边界条件是精 确的、完美无缺的，除非计算域取得比实际流动区域还大，但这是不可能的。 由于非局部边界条件将所有边界点的运动耦合或者将全部过去时刻的运动 与现时刻的运动耦联，因此，就近场波动在时空域的数值模拟而言。即使是对于 线性问题，大计算量的非局部边界条件也难以应用。从6 0 年代开始，对于局部 的边界条件逐渐成为研究的焦点。 12 2 局部边界条件 局部边界条件具有“局部一阵”的特点。即一个边界点在某一时刻的运动仅与 其相邻点在相邻时刻的运动有关，这样可大大简化数值计算。局部边界条件尽管 不能严格的模拟无限域对近场波动的影响，但能够在满足实际需要的糟度下模拟 这一影响。常用的方法有：辐射边界条件、基于特征关系的边界条件和嵌边法。 辐射边界条件通过对单向波方程的近似，将波方程展开，模拟外向波在边界 上的被吸收过程( 实际是传出过程) ，也称为吸收边界条件。随着近似的程度提 高，作为吸收边界条件的偏微分方程的阶数也要提高，而高阶边界条件的稳定性 不易证明，即使能够证明将它用于数值计算时的差分格式也相当复杂。很多工 作致力于寻求形式简单而近似程度高的边界条件“。但是，辐射边界条件存在 着很大的缺陷：它只适用于由双曲型方程控制的线性系统。 基于特征关系的边界条件在不可压缩和可压缩流中广为应用”2 2 1 。它的基本 基于特征关系的边界条件在不可压缩和可压缩流中r 为应用“1 。它的基本 第一章绪论 思想是根据特征线理论，分析穿过边界的波。对于流出计算域的波，其值完全由 计算域内的值确定，可以用单边差分格式求得边界上的空间导数；而对于穿过边 界流入计算域的波，其值则通过由物理边界条件导出的入流及出流波的关系由出 流波算出。p o i n s o t s k t ：e l e ( 1 9 9 2 ) “”建议了一种方法，称为卜s 方程特征 关系边界条件( n a v i e r s t o k e sc h a r a c t e r i s t i cb o u n d a r yc o n d i t i o nn s c b c ) 。 在边界上，采用局部一维无粘假定，即忽略粘性项和平行于边界的对流项，这时 流动的控制方程成为一维波动方程，从而可用特征线理论对传出、传入的波按上 述思想进行分析。但数值模拟的实践表明，在具有大幅值扰动的流场中，特征关 系边界条件常常是不准确的“。 在实现无反射出流条件的方法中，有一种所谓的嵌边法( f r i n g em e t h o d ) ， 也称缓冲区法、海绵区法等”2 5 1 。其主要思想是：在出流处增加一段特殊区域一 一嵌边区。在该区域内，通过在控制方程中引入人工阻尼，使将出流扰动衰减至 足够小，以尽可能地减弱出流反射波的强度。该人工阻尼可以在控制方程中显式 地增加一项”1 ，也可以是仅通过网格密度的改变而引入的数值阻尼。“。这一方法 对不可压缩流的有效性已在实际应用中被证实。本文的主要工作就是研究该方法 在可压缩流中是否也有效，以及该方法中的一些参数应如何选取。 1 3 本文的工作 上文曾经提到，没有任何一种边界条件是完美的。非局部边界条件繁琐复杂， 计算量庞大：辐射边界条件仅适用于线性系统；基于特征关系的边界条件在大扰 动的流场中有可能失效。而基于通常所用的特征关系边界条件所导出的无反射边 界条件实际上并不是人们所要求的条件。其主要原因是，从特征关系角度看，只 要是亚音速区，就总有入射特征波。除非能事先给出此波之值，否则总会有误差。 而事先给出其值一般来说是不可能的。 本文将通过嵌边法对若干可压缩流的试算，考察其作为出流边界条件在可压 缩流中的有效性，并提出了一种衡量边界条件是否有效的评判标准。最后研究了 有关参数的选取问题。 1 4 嵌边法原理 嵌边法的主要思想是：在出流处增加一段特殊区域嵌边区，在该区域内， 经过修正的n a v i e r s t o k e s 方程可将出流扰动衰减至足够小以尽可能地减弱出 流反射波的强度。 第一章绪论 在原计算域要应用此方法的边界外增加一段计算域，在此域内对n s 方程增 加一项特殊的源项h 。经过修正的守恒型n s 方程可以表示为 型+ 旦堡+ 望+ 塑：堕+ 盟+ 里+ 日( 1 1 i ) 一一 li-j o ta x 却 娩瓠 却 a z 日= 五( 坝u e - u ) ，【，= ( 户，删，胆小e ，= 灭箸百+ u2 + v 2 + w 2 ) 式中，罢、a f + j a i g 为对流项，誓、孕和罢l 为粘性项，日为嵌边 c 珥oy(圮(蕊o y ( ” 法对n s 方程的修正项。中，u 是要求解的物理量，玑是预先指定的嵌边区 出流处应取的值。z ( x ) 在嵌边区内总为非负数。 嵌边函数丑( x ) 的选取要保证求解的物理量u 在嵌边区内平缓地趋近以，其 解析表达式为 撒m s ( 专，1 ) 书尹+ 1 _ 1 ) 】 2 ) 式中形函数s ( x ，) 为： s ( x ，矿) = k 。，z 。，石。d ，。，m 的含义如图l 所示。 删 图1 嵌边函数旯( 工) 的典型形状 ( 1 - 3 ) 倒譬 o u + 半 o 一 ， 柚 l ：r l 一2 第二章控制方程及数值计算方法 2 1 控制方程 第二章控制方程及数值计算方法 2 ，1 。1 二维守- 匾型n - $ 方程 本文算例所采用的控制方程均为二维可压缩流动的守恒型n a v i e r - s t o k e s 方 程。在直角坐标系下，该方程写为： 丝+ 丝+ 堑：堡+ 笪 + + 一= o + 2 a f 舐砂 苏 砂 ( 2 一1 ) 其中，u 为要求解的物理量；_ a e 利_ a f 为对流项，压力项包含在内；孕、孥 o x o y 出 o y 为粘性项，包含了热传导项。而： e 。= p p “ p e 。 e = p u p u + p d “v ( 加。+ p - 铲弘，一) = 秤瓦a u 考 驴铲u y + v x ) = 聘+ 期， r 。= ；( - “，+ z b ) = 詈( 一西a u i + z 万a v ) e = ，= 、fllif吖p 神 伽咖 g 一 rv+ o 砀 g 一 掣 fp十 。k 第二章控制方程及数值计算方法 e = 三o2 + v 2 ) + 两p i g 。= 一圾2 一i 1 耳1 瓦万p 习面a t 毋2 一螺2 一i 1 虿1 万丽p ) 万o t 式中“，v ，p ，t ，p ，k ，q ，q 。分别为流体的流向、法向速度，压力，温度，密 度，粘性系数，热传导系数，流向、法向热通量。 状态方程写为： p = p r r ， 其中r 为气体常数。 2 1 2 方程的无量纲化 ( 2 2 ) 以尾缘处的边界层排移厚度6 ，来流速度“。，无穷远处的密度肛，温度k 和粘性系数。为无量纲化时的参考量，则相应的无量纲物理量为( 式中下标“o 。” 表示无穷远处的参数，上标“$ ”表示无量纲的量) ： x2 否y2 ；“2 瓦”2 瓦 ，= i t ，p + 风p 乩：，p = 瓦p ， ( 2 3 ) + 2 老，f 2 瓦t ，肚m = z l ( ，_ - 1 ) r e p r 无量纲数：马赫数m 。：生，雷诺数r e ：垦型，音速：属， a 。 p 。 粘性系数满足s u t h e r l a n d 公式 ( 篙) ， t ， 其中c = 1 1 0 4 k t ，p 为p r a n d t l 数，定义为p = c ，k ，c ，是定压比热，只 在计算中取为常数0 7 2 。y 为比热比，取为1 4 。k 为热传导系数。 第二章控制方程及数值计算方法 将n s 方程无量纲化。为方便起见略去上标“木”，则所得方程形式上与( 2 - 1 ) 一样，且u 、e 、f 的形式不变，但占。，e 中的粘性应力项要除以雷诺数r e ， 具体形式如下： e v = o t 。r e r ，。r e ( u f 。+ v g y ) r e q x 热通量表达式为 e = o r m r e 1 口r e ( u i + v l v y y ) r e 目， 旷越i 一面蒜濮- x ，y 无量纲化后的状态方程变为： 2 2 计算方法 2 2 1 通量分裂 p = p r 彬：。 ( 2 5 ) 为保证计算稳定，对流项应根据j a c o b i a n 矩阵特征值的正负进行通量分裂， 并使用迎风差分格式进行黼胁方向的等黼设j a c o b i a n 矩阵堡o u 吼 其特征值为) ，i = l ，4 ，则可找到特征向量s 使a = s a s 。于是通量可作 以下分裂： 丝：丝型：4 型：塑：坐：旦：塑坠尘壁：旦 融o u3 x 一苏o x苏苏 ：3 s a + i s 一- t u + o s a - i s 一- i u ：o e 一+ + _ o e = _ - ，( 2 - 6 ) 0 x僦0 出 a 的四个特征值分别为 = h 一口，如= 如= “，厶= “+ a 。人= d i a g ( ，丑：，l ，以) 令辟掣，a + = d i a g ( 臂怎鸳确则e 可以分裂为鹃如时】， 第二章控制方程及数值计算方法 其中i = l ，2 ，3 ，4 向 e 2 = 芪= p 咭 考= p 时“+ 7 7 ；a ) ：p叩v(2-7) 爵= p ( 昙( “2 + v 2 ) 砰+ 旌“n + 疗a 2 ( ，一1 ) ) f 印= 篱+ 0 j 一2 置+ 笃) ( 2 y ) 落= k 一岔) ( 2 y ) 【玎；= k + 嚣) ( 2 y ) 按照同样的办法可以得到f + ，与e 2 具有类似的形式。 分裂后，原控制方程( 2 - 1 ) 变为： 一0 u + 堡+ 堡+ 一o f + + 一0 f - ：里+ 堡。a彘 缸却 却m咖 2 2 ，2 差分格式 22 2 1 空间差分格式 对流项采用五阶精度的弱迎风紧致格式： ( 2 8 ) 2 量+ 3 多2 去k 一1 2 瓯+ 3 四- 4 4 站。) ， 。， i s 等+ z 警= 瓦1 、e 。- 一1 2 e l _ 。一3 6 耳+ 4 4 瓦，+ s 甄：) 。 在邻边界点处，则要降为三阶精度的弱迎风紧致格式 l 警+ z 等= - l ( e ：, + 4 e 2 a x一训i 融舐 7 【z 等+ 警= 上k - 4 e 2 a x i 喙，) l 融苏 、” ” 粘性项采用六阶精度的中心型紧致格式 9 ( 2 1 0 ) 第二章控制方程及数值计算方法 1 2 孥+ 3 6 0 。f + 1 2 0 f ，+ , ：塑k 掣出生二越，( 2 - 1 1 ) 幽出出缸 2 2 。2 2 时间差分格式 时间上采取二阶龙格一库塔( r u n g e k u t t a ) 差分格式。控制方程可写作 o = - u ：震辑) ，这里r 是空间导数的离散算子。则二步二阶精度的r u n g e - k u t t a 法 d f 形式如下： 2 3 坐标变换 = “+ 出( 帮) “”+ i = = ”- 等( j z “8 + ) 2 1 2 在计算有关边界层的流场时，由于边界层内物理量在y 方向梯度大，为了保 证计算精度而又不增加很多工作量，y 方向可采用变网格，使边界层内，特别是 靠近壁面处网格较密而边界层外网格较稀。具体方法为用某种变换将计算坐标上 的均匀网格变为物理坐标上的非均匀网格： 设对应于物理坐标( x ，y ) 的计算坐标为( 孝，r ) ，在计算坐标中为均匀网格。 设x 方向的变换为工：x ( 善) ，y 方向的变换为y = y ( 玎) ，则导数_ o f ，荽与导数罢， 以o k o g 兰的关系为 笪：笪丝：堑堡 i 玩8 0 x8 | a 翌：笪塑：望至 谚a no y8 q | 8 h 本文采用了以下形式的坐标变换： i ：_ 也丑卞，n ：4 - 丽- 3 ( 2 - 1 3 ) l + a ( a 一叩) 4 ， ， 其中，= 是上边界处一个网格宽度和壁面处一个网格宽度之比并取为2 0 0 。 砂o 7 7 是从o l 均匀变化。在本文3 ，1 节计算亚音速平板边界层扰动演化的算例中。y l 为y 方向的计算域长度；在3 3 节亚音速平板边界层扰动在尾缘的绕射的问题中， y f 为y 方向计算域的一半。 第二章控制方程及数值计算方法 2 3 边界条件 本文首先通过若干算例来研究嵌边法作为出口边界条件的优缺点。作为比 较，也用了目前在可压缩流直接数值模拟中常用的基于特征关系的边界条件。本 文只在计算域下游出口处用了嵌边法，其他边界( 固壁除外) 则都采用了基于特 征关系的边界条件。 采用基于特征关系的边界条件时，对于流出计算域的波，其值可以完全由计 算域内的信息确定，即用单边差分格式求得；而对于穿过边界流入计算域的波， 其值原则上应由计算域外的信息确定而不能用计算域内的点差分求得。这时要从 物理上考虑来解决这一问题。p o i n s o t & s k l e l e 建议了一种方法，称为n s 方 程特征关系边界条件( n a v i e r - s t o k e sc h a r a c t e r i s t i cb o u n d a r yc o n d i t i o nn s c b c ) 。 其基本思想是：忽略n s 方程中的粘性项及与边界平行的对流项后应用特征线理 论的方法，因此所得关系式又称局部一维无粘关系l o d i ( t h el o c a l 0 n e d i m e n s i o n a li n v i s c i dr e l a t i o n s ) 。 他们将对流项做如下特征分解： 妒l 妒2 仍 5 = 0 钺詈一胪 碘2 塞一菪d 亡d 亡 o v 琵 九( 塞+ 肛 对流出计算域的项，可用计算域内的点做单边差分。但如果丑 0 ，则纪是流入 计算域的，应由外部的状态确定。考虑到计算域外的压力往往是事先就可知道的， 所以他们建议令纯= 2 , k ( 1 一m ：) ( p 一凡) ，其中k 为反射系数。这种边界条件将 无穷远处压力的信息带回亚音速边界，反映了其对计算域的影响。k - - - - o ，即被认 为是完全无反射边界条件。本文中作为比较的算例就采用这样的边界条件。 第三章嵌边法的应用算例 第三章嵌边法的应用算例 我们以嵌边法为出口边界条件算了三个可压缩流的算例，分别为亚、超音速 平板边界层扰动波的演化；超音速流中涡的传播；亚音速平板边界层中的扰动在 尾缘的绕射，并与采用基于特征关系边界条件所得到的结果作了比较。 3 1 亚、超音速平板边界层扰动波的演化 3 1 1 亚音速平板边界层扰动波的演化 3 1 1 1 基本流的计算 当计算域较平板短时，就有扰动波是否能顺利从计算域下游出口传出去的问 题。 取来流马赫数为o 7 ，无穷远处的气体参数值相当于高空1 0 0 0 米处的标准大 气值，即咒= 2 8 1 7 足，c 。= 3 3 6 4 m s ，p 。= o 8 9 8 7 6 1 0 5 p a ，p 。= 1 1 1 l 昭卅3 ， 。= 1 7 5 8 x1 0 。p a - s 。取某处边界层排移厚度巧作为长度无量纲参考量，其它 物理量无量纲化时的参考量分别为：无穷远处流体速度u 。；无穷远处流体温度 l ：无穷远处流体密度成；压力p 。u 。2 ；无穷远处流体粘性系数u 。后 面文中提到的最，如没有特别说明均采用无量纲的。这样得到的雷诺数 璺盟：9 0 0 0 ，对应的排移厚度为占：6 0 4 8 l o - 4 m ，它相当于自平板前缘起 。 1 4 7 4 m 处由相似性解得到的排移厚度。 计算基本流，可首先求得相似性解作为初始流场，然后保持入口不变，求解 n a v i e r s t o k e s 方程，直到定常为止。这样就得到亚音速平板边界层的基本流场。 3 1 1 2 嵌边法算例和参考算例 我们所研究的问题是边界层内的扰动演化。由于计算域有限，扰动将穿过下 游边界。问题是何种边界条件可以保证扰动穿过边界时无反射，即不会产生额外 的扰动。 在平板基本流场的入口处引入t - s 波。将流场分解为基本流和扰动两部分如 第三章嵌边法的应用算例 下式 “= u o + “ v = + v p = 只+ p t = 瓦- i - t ( 3 1 ) 其中u ，v ，p ，t 为流速、压力及温度，砜，昂，五为基本流解，“。，v ，p ，丁。为 扰动。t s 波扰动在入口处的表达式为： “= a ( u ，( y ) c o s ( 0 1 ) + 虬( y ) s i n ( a 】t ) ) ”2 口( _ ( y ) 。0 8 ( c a ) + v j ( y ) s i n ( o x ) ) ( 3 - 2 ) p = a ( p ，o , ) c o s ( c a ) + p ，( y ) s i n ( o x ) ) 丁。= 口( z ( y ) c o s ( 耐) + 正( y ) s i n ( 耐) ) 其中a 为幅值，“，( y ) ，u ，( y ) ，v ，( y ) ，q ( y ) ，p ，) ，a ( y ) ，r a y ) ，正( 力分别为对应 于入口处基本流的频率为的归一化后的t s 波的速度、压力与温度特征函数的 实部和虚部，归一化的条件为当a = 1 时，h 的模的最大值为1 。 取计算域长度为1 0 0 ，另加长3 0 的嵌边区。在嵌边区内采用修正的n s 方 程进行计算。计算足够多整周期，在扰动波已穿过嵌边区后，即停止计算。 为了考察嵌边法的有效性，又计算了一个参考算例与之比较：取足够长的计 算域，使得计算同样多整周期后，扰动波虽已穿过有效计算域出口( 相对于嵌边 法算例的嵌边区入口) ，但下游还有足够长的计算域，以致扰动还未到达出口， 从而下游的边界条件对计算结果不会有什么影响。比较两个算例在有效计算域 ( 不含嵌边区) 内两者的结果，即可断定嵌边法是否行之有效。参考算例和嵌边 法算例示意图见图3 - 1 。 嵌边法算例b s 参考算例t - s 匮强 图3 1 参考算例和嵌边法算例示意图 第三章嵌边法的应用算例 在入口处所加t - s 波的振幅为0 0 l ，频率为o 1 5 4 。嵌边法算例中计算域( 包 含嵌边区) 的流向长度法向长度为1 3 0 2 0 ，网格数为1 3 0 x 2 0 0 ，法向采用 2 3 节式( 2 1 3 ) 将，7 上的均匀网格变为y 上的非均匀网格，越靠近壁面网格越 密。 嵌边函数的选取如图3 - 2 所示。有关参数为：嵌边区长度上，= 3 0 ，a 。= 1 0 ， a = 5 1 6 l ，a 俐= 3 1 6 l r ，出口边界条件采用式( 3 3 ) 所示的外推格式计算边界 点的导数。对于函数r ， ( 乱= ( 筑，一s ( 乱+ ，( 筑 ( 3 3 ) 其中，i n 、i n 一1 、i n 一2 、i n 一3 是离散点沿z 方向的点号，i n 是出口边界上点的点 号。 3 1 1 3 数值模拟结果 图3 - 2 计算中采用的嵌边函数 x o “ 第三章嵌边法的应用算例 x ( b ) 图3 - 3 两种算法所得横向速度v 的比较 图3 3 ( a ) 、( b ) 分别给出了距壁面1 8 5 和3 9 l 处扰动速度v 的分布。图 中标有r e f e r e n c e 的表示参考算例结果，标有f r i n g e 的表示嵌边法算例结果( 其 它图同) 。可以看出，二者在有效计算域中吻合得很好，而在嵌边法算例的嵌边 区内，扰动波有明显的衰减，扰动值逐渐改变至接近于零。 我们试用文献 t 9 中的无反射出流条件计算该问题，反射系数在0 与0 3 之 间调整，但总不能得到令人满意的结果。实际上，该条件是利用出口压力与环境 压力差来确定入射特征波值，但在非定常情况中无法保证其值与真实情况一致。 3 1 2 超音速平板边界层扰动波的演化 来流马赫数为4 5 ，气体参数取高空5 0 0 0 米处的标准大气值，无量纲化情 况同亚音速平板边界层问题，则雷诺数为9 0 0 0 0 。嵌边区长度为6 ，其他参数选 取同3 ，1 1 ，嵌边区出e l 用特征关系边界条件给出。可以得到与亚音速流相同的 结论。把扰动幅值增大至0 0 5 ，或把扰动波换成振幅、频率各不相同的2 - - 3 个 t s 波之和，都能得到令人满意的结果。见图3 - 4 和图3 - 5 ( 图中3 0 3 6 为嵌边 区) 。 第三章嵌边法的应用算例 图3 4 扰动幅值增大至0 0 5 的结果比较 x 图3 - 5 扰动波为3 个t - s 波的结果比较 实际上，如出口处全是超音速，则不会有入射特征波( 不考虑粘性时) ，因 而出口边界条件很容易处理。而超音速边界层近壁面处总是亚音速，但所占区域 很小，其影响随着马赫数提高而减小。因而超音速边界层出口边界条件的处理比 亚音速边界层要容易一些。 3 1 3 衡量边界条件的评判标准 上面提到，比较嵌边法算例和参考算例在有效计算域内的结果，可以断定嵌 边法是否有效。我们认为，不可能事先知道出口处入射特征波之值是无法做到完 全无反射的根本原因。因此比较出口处从计算结果反推的入射波和实际应有的入 射波是鉴别方法好坏的更有效指标。而实际应有的入射波则可以通过加长计算 域，使得在同一计算时间内扰动尚未到达出口处，再反推原出口处的负通量而得 到。下面将从这一角度比较嵌边法和特征关系边界条件。 第三章嵌边法的应用算例 ( a ) 正通量( b ) 负通量 图3 - 6 在x = 3 0 处，三种方法所得正负通量的比较 图3 - 6 中给出了3 1 2 节算例中第4 5 个无量纲时间实际出流位置处三种算 法，即用加长计算域( 参考算例) 、嵌边法及文献 1 9 中方法( 结果用n s c b c 表 示) 压力反射系数取为o 时算得的x 方向对流项的正负通量沿y 的变化。可以看 出，用特征关系边界条件得到的结果与参考算例相差较大，而嵌边法算例和参考 算例的正负通量值在图中几乎没有差别。 3 2 超音速流中涡的传播 不。 文献 1 9 中有一个涡向下游传播穿过出流边界的算例。初始的涡如图3 7 所 + ， 0 一， _ “0 一厂弋、 2 z _ 、弋 图3 7 初始涡位置 平均流为超音速流，m a u1 1 。t = o 时，速度场用不可压缩无粘流的流 c 函数( 坐标原点在涡中心) 来定义： 第三章嵌边法的应用算例 ( u 。f l u 。o j + 万1 d v