（流体力学专业论文）不可压缩流体边界层条纹结构无粘稳定性的数值研究.pdf

摘要 在工业设计中，包括飞行器设计与发动机等的设计，特别是在机翼与涡轮叶 片的设计中，流动发生转捩的位置是重要的设计参数，这就需要对流动的转捩性 质以及流场的稳定性特性加以分析和认识。对背景湍流度很大的流动，边界层中 有时会产生沿展向周期分布的流向速度的条纹结构。这时在预测转捩发生的位置 时，传统的预测方法会失效，这就需要对边界层的条纹结构的稳定性进行描述。 为了研究不可压缩边界层的条纹结构在边界层转捩过程中所起到的作用，本 文采用时间模式的数值模拟方法，利用无粘的线性扰动方程进行研究。法向采用 c h e b y s h e v 配置点法，流向和展向采用f o u r i e r 谱方法，时间上采用三阶精度的 差分格式进行计算。首先在流场中放入一些由b l a s s i u s 剖面所得到的线性解的扰 动，经过足够长的时问，就可以得到不稳定的波，并可以求得其增长率、分布特 性、以及随流向波数的变化。 通过对基本流场的计算和分析，我们得到以下几点结论： ( 1 ) 使用伪谱c h e b y s h e v 配置点法计算边界层是有效的。 ( 2 ) 边界层条纹结构中存在着增长波，并且其随时间呈指数增长，其扰动的 传播呈现明显的周期性变化。 ( 3 ) 对不同的流向波数，其特征值随着流向波数的增加表现出先递增后递减 的性质，在口= 0 4 附近具有最大的增长率，其对应的频率为( 0 2 1 8 6 + 0 0 0 1 4 5 7 i ) 。 ( 4 ) 分别给出了对应于不同流向波数的特征函数分布。 关键词：条纹结构稳定性c h e b y s h e v 配置点法数值模拟 a b s t r a c t i nt h ee n g i n e e r i n gd e s i g n i n g ，c o n s i s t i n go ft h ed e s i g no ft h ea i r c r a f t sa n dt h e e n g i n e s ，t h et r a n s i t i o nl o c a t i o no ft h ef l u i di san e c e s s a r yc r i t i c a ld e s i g np a r a m e t e r , e s p e c i a l l yo na l la i r p l a n ew i n go rat u r b i n eb l a d e i tn e e dh a v eag o o da n a l y s ea n d g o o dk n o w l e d g et ot h ec h a r a c t e ro ft h et r a n s i t i o na n dt h es t a b i l i t yo ft h ef l o w f o r h i 曲l e v e l so ff r e e - s t r e a mt u r b u l e n c e ，s t r e a k ys t r u c t u r e si nb o u n d a r yl a y e rw i t ht h e s p a n w i s ev a r i a t i o n sa l t e r n a t i n gw i l lb eg e n e r a t e d a se s t i m a t i n gt h et r a n s i t i o nl o c a t i o n ， t r a d i t i o n a lm e t h o di su s e l e s s t h e r e f o r e ，i tn e e dh a v i n ga d e s c r i p t i o nt ot h es t a b i l i t yo f t h es t r e a k ys t r u c t u r e si nb o u n d a r yl a y e r an u m e r i c a li n v e s t i g a t i o ni sm a d ef o rt h ee f f e c to f s t r e a k ys t r u c t u r e si nt r a n s i t i o n b yi n v i s c i dl i n e a rd i s t u r b a n c ee q u a t i o nw i t ht e m p o r a lm o d ec o n s i s t i n go ff o u r i e r e x p a n s i o n s i nh o r i z o n t a l d i s c r e t i z a t i o n ，c h e b y 7 s h e vc o l l o c a t i o np o i n ti nn o r m a l d i s c r e t i z a t i o na n dt h i r d o r d e rd i f f e r e n c es c h e m e si n t e m p o r a ld i s c r e t i z a t i o n a s i n i t i a l i z i n gt h ed i s t u r b a n c ew i t ht - sw a v e ，t h ei n s t a b l ew a v e , i t sg r o w t hr a t e ， d i s t r i b u t i o n ，a n dt h ev a r i a t i o nw i t ht h es t r e a m w i s ew a v en u m b e r sc a nb er e c e i v e df o r a l o n ge n o u g ht i m e t h ef o l l o w i n gc o n c l u s i o n sa r ef o u n d ： ( 1 ) c h e b y s h e vc o l l o c a t i o np o i n ti nn o r m a ld i s c r e t i z a t i o nf o rb o u n d a r yl a y e ri s v a l i d ( 2 ) t h ew a v e se x p o n e n t i a lg r o w i n ga n dp e r i o d i ca l t e r n a t i n ga r ei ne x i s t e n c ei n s t r e a k ys t r u c t u r e si nb o u n d a r yl a y e r ( 3 ) f o rd i f f e r e n ts t r e a m w i s ew a v en u m b e r s ，e i g e n v a l u e sa r ei n c r e m e n t a la tf i r s t a n dd e c r e m e n t a ll a t e rw i t ht h es t r e a m w i s ew a v en u m b e ri n c r e a s i n ga n dt h em a x i m u m r a t eo fg r o w t hi si ne x i s t e n c ea tn e a ra l p h a = 0 4 ；i t sf r e q u e n c y i s ( o 218 6 + 0 0 0 1 4 5 7 i ) ( 4 ) t h ed i s t r i b u t i o n so fe i g e n f u n c t i o nw i t hd i f f e r e n ts t r e a m w i s ew a v en u m b e r s a r er e c e i v e d k e yw o r d s ：s t r e a k ys t r u c t u r e s ，s t a b i l i t y , c h e b y s h e v c o l l o c a t i o n ， n u m e r i c a ls i m u l a t i o n 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的 研究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表 或撰写过的研究成果，也不包含为获得鑫盗盘堂或其他教育机构的学位或证 书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中 作了明确的说明并表示了谢意。 学位论文作者签名 签字日期：知7 年多月，日 学位论文版权使用授权书 本学位论文作者完全了解苤壅盘鲎有关保留、使用学位论文的规定。 特授权丕鲞盘堂可以将学位论文的全部或部分内容编入有关数据库进行检 索，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校 向国家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名导师签名： 签字嗽砷年多月r 乒日 签字吼 彤z 矽中乡月，仁 第一章绪论 1 1本文的目的和意义 第一章绪论 流体在某种条件下，流动状态由层流变为湍流，我们称其为转捩。在流体力 学的研究领域，转捩问题是一个古老而又重要的问题，它被认为是流体力学中剩 下的难题之一。它曾吸引了不少著名的力学和物理学家参与研究，但至今仍然是 流体力学中有待进一步研究的最重要的问题之一【l 】- 【3 】。 对于工业中的大多数流动来说，其边界层所产生的摩擦力在全部阻力中占有 相当大的比例。由于层流边界层所产生的摩擦阻力远小于湍流边界层所产生的阻 力，因而，人们希望通过优化，能够控制边界层的流动，使其保持层流状态或者 推迟进入湍流状态，以达到使得摩擦阻力更小的效果。此外，在诸如发动机中， 为了促进燃料与空气充分混合，希望流动是湍流状态，这样可以节省能源和燃油， 其意义也是显而易见的。对于机翼或者涡轮叶片来说，从层流到湍流的转捩位置 是一个重要的临界设计参数。这些具有重大实际意义的具有挑战性的课题都与边 界层流动的稳定性及其转捩问题的研究密不可分。 一般认为由层流到湍流的转捩的发生是由于基本流的早期不稳定引起的。而 从一种状态到另一种状态，是由于原先流动中的扰动自然演化的结果，这是由基 本流的一些特性所决定的。分析基本流的这种特性和研究扰动的演化过程，也就 是稳定性研究的主要工作。 传统的对于转捩的预测主要是基于幅值的指数增长理论的e n 方法，这是一 种半经验方法，只适用于预测所谓自然转捩问题，即由小扰动引起的转捩，而不 适用于背景湍流度较大的所谓“b y p a s s ”转捩情况。对于高湍流度引起的“b y p a s s ” 转捩问题，我们就需要从较高背景湍流度引起的基本流中寻找不同于传统理论的 新的不稳定性机理。实验和理论都发现，某种外界扰动下，将会使得流场产生沿 展向分布的条纹结构，对这种条纹结构的不稳定性特征及其在转捩中的作用的研 究，有助于我们更好的理解和描述这种在较高背景湍流度情况下的不稳定波增长 以及转捩的特性。 因此，本文尝试对由较高背景湍流度在边界层流中引起的条纹结构的进行稳 定性分析，以期得到一些有意义的结论，以便了解这类条纹结构在“b y p a s s ”转 捩中的作用。 第一章绪论 1 2流动稳定性问题概述 一般认为，在背景湍流度很小的情况下，从层流到湍流的转捩大致经历如下 过程：首先由外界的扰动影响到边界层内部，引起初始的扰动；然后，由于边界 层的不稳定性特性，扰动在向下游传播的过程中逐步增长，开始的扰动振幅很小， 扰动波之间在开始阶段基本上是独立演化的，就称之为“线性稳定性阶段”；随 着扰动的逐步增长，达到了有限幅值状态，它们相互之间的非线性影响已经不能 忽略，此时流动就到达了“非线性”的阶段；进一步的，非线性作用的逐步增强， 从而从层流转捩到湍流。 1 2 1 线性稳定性 线性稳定性理论又称作小扰动线性化理论，其主要目的在于确定流动在什么 样的条件下对于无穷小的扰动是稳定的，在什么样的条件下，扰动是不稳定性的。 在r e y n o l d s 通过实验证实圆管流动中存在层流和湍流两种不同的状态之后不久， 就有人提出了是由于层流失稳导致转捩的理论猜想。 1 9 世纪，就已经有了无粘流的稳定性理论。从无粘稳定性方程出发，r a y l e i g h 成功的得出了重要的关于层流速度剖面稳定性的定理，它指出具有拐点的速度剖 面是不稳定的【4 】。 2 0 世纪，o r r 和s o m m e r f e l d 针对平行平板间的流动建立模型，假设波的振 幅足够小，忽略非线性项，建立起了研究平行流稳定性的小扰动方程，即 o r r - s o m m e r f e l d 方程( 以后为求简便，称其为o s 方程) 。h e i s e n b e r g 第一个从 理论上证实在r e y n o l d s 数很大的时候有不稳定解，但他没有具体算出临界 r e y n o l d s 数。t o l l m i e n 和s c h l i c h t i n g 进一步发展了h e i s e n b e r g 的渐进方法，具体 计算了小扰动的解，预测了小扰动会以进行波的形式出现，以后这类问题的解就 被称为t o l l m i e n - s c h l i c h t i n g 波( 以后简称t - s 波) 。对于求解o s 方程，在数学 上遇到了很大的困难，经过h e i s e n b e r g ，t o l l m i e n ，s c h l i c h t i n g 以及林家翘等人 的相继努力，到1 9 4 5 年才使解答该特征值问题的渐进匹配法有了坚实的基础。 从提出方程到最终求解该方程，中间历时2 0 年的时间，足见该问题的复杂性 5 】。 线性稳定性理论不能预测层流转变为湍流的全部过程，但它可以指出何种速 度剖面是不稳定的，分辨出增长最快的小扰动频率，并指出如何改变控制流动参 数以推迟转捩的发生。在不可压缩流动中，稳定性理论对求解不稳定波的特点， 以及边界层失稳的细节，发挥了十分重要的作用。 2 第一章绪论 1 2 2 非线性理论 到2 0 世纪5 0 年代，有关流动稳定性的线性理论已经成熟，但还不能解决转 捩问题，仍然需要非线性理论。 1 9 4 5 年，l a n d a u 已经提出一种失稳而导致转捩的非线性机制。1 9 6 0 年，s t u a r t 提出了有名的弱非线性理论，将l a n d a u 的设想具体化。随后有人又做了改进。 周恒和赵耕夫 6 】在1 9 8 2 年将非线性振动理论中的k b 方法应用于流动稳定性的 研究中，其实质和s t u a r t 的方法一样，但其形式更带有普遍性，更便于应用。但 后来周恒【7 8 】通过详细的研究发现：弱非线性理论对某些问题，特别是对转捩 或自由剪切流中涡的演化研究的结果不能令人满意。流动稳定性弱非线性理论中 存在的问题，使其不能很好的研究扰动的演化，最主要的问题是处理平均流修正 的方法不正确，导致不能处理一般的初值问题。他并且提出了改进方法，其结果 也得到了数值模拟的很好的证实。经过不懈的努力，弱非线性理论得到了完善和 发展。 在从理论上对非线性问题进行研究的同时，人们从实验中发现了转捩过程中 重要的非线性现象。k l e b a n o f f 等人在1 9 6 2 年首次对平板边界层的转捩做了系统 的研究，发现虽然引入的扰动是二维的，但在沿流向演化时，不可避免地要出现 三维扰动。虽然后来发现，实际上引入的初始扰动严格来说并不是二维的，但按 照线性理论，一般来说二维波更不稳定，所以三维波的出现无法用线性理论解释。 他们所发现的三维波在流向的波数和引入的二维扰动的波数是一样的。后来 k n a p p 等人在实验中又发现了流向波数为引入的二维波波数的一半的三维波，即 具有亚谐波的性质。s a r i e 等人用可视化方法直观地显示出了三维扰动的结构， 发现了三维波和引入的二维波的波数有三种不同的比例关系。一种就是 k l e b a n o f f 所发现的那种，其结构如图1 2 1 ( 口) 所示。、另两种三维波则都是亚 谐波，其结构都如图1 2 2 ( 6 ) 所示，而三维波的展向和流向波长之比则不同。一 种可由c r a i k 的三拨共振理论解释，另一种可由h e r b e r t 的二次失稳理论解释。 因此三种不同的三维波分别被称为k ( k l e b a n o f f ) 一型，c ( c r a i k ) 型和h ( h c r b c r t ) 型 5 】。 图1 2 3 烟线演化示意图 3 第一章绪论 转捩过程r f l 重要的非线性现象在实验室的发现又促进了非线性理论的发展。 c r a i k 在1 9 7 1 年提出了共振三波理论：b e n n e y 等人在1 9 8 1 年提出了直接共振理 论：h e r b e r t 在1 9 8 3 年提出了二次失稳理论；k a c h a n o v 在1 9 8 7 年提出了更一般 的共振理论。k a c h a n o v 及其同事们【9 】通过实验验证了亚谐共振理论，并证明了 k 一型转捩和n 型转捩有相同的亚谐共振机制。通过比较不同的实验数据， k a c h a n o v 1 0 认为，不管近壁剪切流的形式相同与否，也不管扰动的初始条件如 何，在转捩的突变( b r e a k d o w n ) 过程中，起主导作用的扰动演化历程是普遍相 同的。b a k e 1 1 】等人的实验结果表明，不同类型的转捩，其后期过程的表现形式 是一致的。 另一方面，不得不需要指出的是，在二十世纪7 0 年代与理论分析的同时， 开始利用计算机开展直接求解n s 方程的数值模拟工作。f a s e l 在一个增长的边 界层中得到了一个二维解。后来，o r s z a g 和k e l l s 与o r s z a g 和p a t e r a 1 2 得到了 具有周期边界条件的三维d n s 的结构，在一个壁边界剪切流中捕捉到k l e b a n o f f 型三维变形。通过对d n s 的数值分析，证实了线性机制的存在，并揭示出：三 维模态随时间以指数规律增长。在此推动下，基于f l o q u e t 理论的二次稳定性分 析才得到了发展，并且从这个分析得到的结果成功证实了实验结果，人们才对 k l e b a n o f f 的实验有了真正的理解。 1 3 由条纹结构( s t r e a k ys t r u c t u r e s ) 引起的不稳定性问题的研究 进展 机翼或者涡轮叶片的流向方向上的从层流到湍流的转捩位置是一个重要的 临界设计参数。e n 方法是一种经常被采用的转捩预测方法，它是一种半经验方法， 通过流动稳定性理论，利用振幅的指数增长来判断转捩位置。该方法使用于判断 由小扰动引起的自然转捩，然而，通常在背景湍流度较高的流动之中，初始扰动 足够大时，边界层流中转捩的位置可能提前，甚至发生在不稳定区域之前。因此， 人们开始注意研究这类不同于自然转捩的所谓“b y p a s s ”转捩。 实验表明在自由流湍流度较大的情况下，边界层中的流向速度会沿展向产生 周期性的快慢条纹结构 1 3 】。b r a d s h a w 1 4 澳1 量了表面剪切力的展向变化，并且 讨论了这种展向变化与风洞阻尼网的特性之间的联系。此后，c r o w 1 5 使用渐进 匹配展开法从数学上得到了二者之间的更严谨的联系。 此外实验表明，在自由流湍流度较大时，包括p o i s e u i u e 流和边界层流在内 的大多数流动，发生转捩的临界r e y n o l d s 数都小于线性理论预测的结果。 k e n d a l l 1 6 】，g u l y a e v 1 7 ，以及a l f r e d s s o n 和m a t s u b a r a 1 8 】做的边界层流动的可 视化实验，发现在背景湍流度较大的时候，边界层在发生转捩之前未观察到t - s 4 第一章绪论 波的存在。b a k c h i n o v 1 9 等于1 9 9 5 年的实验中，通过控制粗糙度产生了沿展向 周期变化的条纹结构，研究结果认为条纹结构失稳引起的转捩与哥特( g o t h e r ) 涡 和横流不稳定性产生的转捩的情况类似。 事实上，通过外界的湍流度和壁面的粗糙度都可以在边界层中产生沿展向周 期性变化的条纹结构。对于这类条纹结构，因为速度分布在空间上存在拐点，因 此可能存在无粘不稳定性。许多人开始研究这种条纹结构在转捩中所起的作用。 1 9 9 9 年，h e n n i n g s o n 与a n d e r s s o n ，b e r g g r e n 等人 2 0 】利用理论分析的办法， 在一定的假设下，给出了由自由流湍流度引起的展向条纹结构的流向速度分布， 见图1 3 1 。 纠獬 图1 3 1 展向条纹结构的流向速度分布( 实线为理论结果，符号为实验结果) 2 0 0 1 年，a n d e r s s o n ，b r a n d t ，h e n n i n g s o n 等人 2 1 禾l j 用无粘方程研究了该条 纹结构的无粘稳定性问题及不同条纹幅值对稳定性的影响。w a l e f f e 2 2 2 3 】， r e d d y 2 4 等人对于槽道流中条纹结构的稳定性进行了研究，发现滑移或无滑移 边界条件对稳定性的研究没有很大的不同。b o t t a r o 与k l i n g r n a n n 2 5 1 9 9 6 年也对 这种条纹结构的无粘稳定性及粘性所起的作用进行了研究。吴雪松和罗纪生 2 6 】 对流场中所存在的某种“修正t - s 波 ( m o d i f i e dt sw a v e s ) 的情况进行了研究， 结果表明，对于小幅值的条纹构，修正的t - s 波的增长率要大于由b l a s s i u s 解得 到的t s 波的增长率。 研究表明，如果条纹结构的幅值足够大，将存在不稳定波，它的发展将有可 能引起提前转捩。人们也将这类由条纹结构引起的不稳定称为条纹结构的“二次 失稳”。 5 第一章绪论 1 4本文的主要工作 目前，人们的研究中所采用的条纹结构大多是简化的理论模型，其中用了许 多假设。其实，可以根据边界区方程，直接计算出自由流中湍流度在边界层中引 起的条纹结构，这样可以使理论更完整。本文采用的条纹结构正是利用这样的方 法得到的，该条纹结构是由吴雪松教授提供的。在此基础上，本文研究了边界层 中条纹结构在边界层转捩过程巾所起到的作用。 通过分析可以发现带条纹结构的基本流的速度存在拐点，因此可能存在无粘 不稳定性。基于这些考虑，本文将采用时间模式的数值模拟的方法，利用无粘的 线性扰动方程，对带条纹结构的基本流进行不稳定性的研究，希望能得到不稳定 的扰动，并对不稳定波的增长率，不稳定波沿壁面法向的分布规律，以及流向波 数对不稳定波的影响等问题进行研究。 数值方法上，法向采用c h e b y r s h e v 配置点法，流向和展向采用f o u r i e r 谱方 法进行计算，时间上采用三阶精度的差分格式进行计算。 在第二章中给出基本控制方程的推导，以及数值计算所需要的具体格式和边 界条件。 第三章中，通过数值模拟的方法研究流场的稳定性。首先在流场中先放入一 些由b l a s s i u s 剖面所得到的线性解得到的扰动，计算其随时间的演化。经过足够 长的时间，就可以得到不稳定的波。然后就可对其增长率、分布特征、以及随流 向波数的变化进行分析。 第四章给出主要结果的分析与讨论。 6 第二章控制方程与计算方法 第二章控制方程与计算方法 2 1计算选用的坐标系 本文的研究对象为平板边界层，故选用如下坐标系( 见图2 1 1 ) ： u z 图2 1 1 本文计算所使用的坐标系 其中流向方向为x 方向，法向方向为y 方向，展向方向为z 方向。 2 2计算用控制方程 2 2 1 基本控制方程 将不可压缩边界层中流体的速度、压力写成基本流和扰动的形式： z ，= u + “，= v + 矿，w = w + w ，p = p 十p ，考虑到基本流速度具有流向条 纹结构，设u = u ( y ，z ) ，v = 0 ，矽= 0 ，扰动速度和压力记为“，矿，w ，p ， 带入连续性方程和n s 方程，可以得到基本流连续性方程和n s 方程，略去扰 动的非线性项和粘性项，可以得到有流向条纹的不可压缩边界层中扰动所满足的 线性方程( 以下为求简洁，略去表示扰动量的上标) ： 7 第二章控制方程与计算方法 d 甜o vo w 一 + + = 0 出砂 比 锄。，锄a u ，a u 一一劲 + c 一十v + w = 一二 西 缸 砂 g z 0 x ( 2 1 ) a v + u 尘：一望 一 、， 趺瓠 印 丝+ u 坐：一望 选取来流速度虬为特征速度，边界层b l a s s i u s 解定义的长度单位 万= t 丽为特征长度，其中x 为考察位置距离平板前缘的距离，沩流体的 动力粘性系数。用上述特征量对方程进行无量纲化。上标“”表示无量纲量， 则有如下关系式： x * 。= x 8 ，, y ，。= 。y 8 , z * = z 蹦2p 慨哦 ( 2 2 ) u = u 虬，t = t l ( 8 虬) 、 雷诺数定义为：r e = u 。8 1 v 对方程( 2 1 ) 进行无量纲化，得到无量纲的无粘线性扰动方程，无量纲化后的 方程形式与( 2 1 ) 一样，为方便起见略去上标“，。 2 2 2基本方程的化简与整理 方程( 2 1 ) 是一个含有u , v , w , p 四个未知数的方程组，为了求解和分析的方便， 我们对方程进行化简。注意到泸吣，z ) ，p 保留，化简起来比较方便。 三个方向上的动量方程分别对x ，y ，z 求导数，然后再相加，可以得到： 皤+ u 昙) 沁+ 。+ 心) + 2 q 屹+ 2 址嵋= 一卸 ( 2 3 ) 考虑连续性方程，有： 卸+ 2 u 匕+ 2 也嵋= 0 ( 2 4 ) 然后，利用关于1 ，和w 的动量方程对x 求导数，z 1 0 匕和叱，9 v l ：, f 弋a ( 2 4 ) ， 最后得到一个关于p 的偏微分方程： 8 第_ 章控制方程与计算方法 ( 昙+ u z v - 2 u y - 2 u z 攻= o ( 2 5 ) 该方程就是带条纹结构的无粘扰动的线性方程，可以通过数值模拟求解该方 程研究扰动的演化。1 为了方便进行求解，我们将基本流【，( y ，z ) 写成u ( y ，z ) = ( y ) + a y , ( y ，z ) 的 形式，( y ) 为边界层的b l a s s i u s 解，4 u ( y ，z ) 为由于外界湍流度在边界层内引 起的条纹结构。并将其代入到方程( 2 5 ) o e ，整理，可以得到： 其中， ( 昙州y ) 昙陟2 吲y ) 茜p = 4 q ( 2 6 ) q = 2 u t ，( y ，z ) 妥p + 2 吣妒) 兰p u ( y ，z ) 昙卸 ( 2 7 ) y ( y ，z ) 蒜p + 2 u ：( 妒) 蠢p u ( y ，z ) 盖卸 ( 2 7 离散方程进行计算时，等式左边项可以用隐格式进行计算，等式右边项用通 常处理非线性项的方法进行处理。 2 3计算方法 在本文中，基本流在展向是周期的，因此，本文在流向和展向上采用傅立叶 谱方法进行计算；在法向上，使用c h e b y s h e v 配置点法；在时间上，则采用三阶 精度的差分格式进行计算。 2 3 1 流向和展向的计算方法 设基本流u ( y ，z ) 在展向上可以用傅立叶级数展开。即可以把基本流写成如 下形式： h u ( y ，z ) = ( y ) + 4 u ( y ，z ) = ( y ) + 4 ( u 。( j ，) + u 。( y ) ( 秒+ e - 棚肚) ) ( 2 8 ) m = l 其中，是外界湍流度在边界层内引条纹结构的基本展向波数，m 为条纹 结构的展开项数。根据基本流的特征，扰动在流向和展向上采用傅立叶谱方法进 行求解，在流向和展向上做傅立叶展开，于是控制方程中的p 和非线性项，也就 可以写成如下形式： 9 第二章控制方程与计算方法 n ：i l p ( x ，y ，z ，t ) = e i a 。p ( y ，坍，t ) e t m 庳 他 2 肥1 2 q ( x ，y 石f ) = p 胁o ( y ，聊，f ) 抄 n z 2 ( 2 9 ) ( 2 1 0 ) 展向的计算域为：t = 2 z r p 设解为( 2 9 ) 的形式的原因是希望得到在流向为行波形式的解，口为给定的参 数。 对于方程右端项q 的计算，采用伪谱方法，即先将谱空间的3 ( y ，m ，t ) 变换到 物理空间配置点上，和条纹结构的速度乘积求出q 的值，再把它变换回谱空间中。 2 3 2 法向的计算方法 对于边界层流动来说，其物理区域是半无限区域，本文法向采用c h e b y s h e v 配置点法进行计算，它可以保证计算中的系数矩阵的阶数比较低，而且在法向方 向上有很高的计算精度。由于c h e b y s h e v 多项式的定义域为s 卜l ，1 】，这就需要 做坐标变换，将物理区域y 0 ，佃) 映射到计算域s 【一l ，1 】中。这可以通过不同 的数学变化来进行，有人对不同的变换及参数下，通过对给定条件的o s 方程的 求解进行了比较，验证了代数变换比指数变换有更高的精度。因此，本文采用代 数变换进行计算，其变换表达式如下： 1 一j 1 ，= 一 。 1 + s 这样就将y 0 ，佃) 映射到计算域s “一1 ，1 】中。s = l 处对应于壁面，s _ - 1 处 对应于无限远处。采用记号d = 善，就有如下网格变换中的导数关系： 凼 杀= g d ；万0 2 = 9 2 d 2 + g 警。 ( 2 1 1 ) 西鲫一黜 其中g = 害= 一半，塞_ g ，一( 川) 。 在法向的计算域上布置配置点，其位置为： s n = c o s ( x p m ) ，刀= o ，以 在这些配置点上，使用c h e b y s h e v 多项式对其进行展开，有： 3 ( y ，m ，f ) = p ( 刀，m ，f ) 7 ：，【s ( y ) 】 牙( y ，m ，f ) = 盆( ，聊，f ) c 【s ( y ) 】 1 0 ( 2 1 2 ) ( 2 1 3 ) 第二章控制方程与计算方法 2 3 3 时间离散 为求解方程( 2 6 ) ，将方程左边看作线性部分，而将方程右边看作非线性部分 进行计算。线性部分采用半隐格式，非线性部分使用显式格式。 方程( 2 6 ) 可以统一写成如下的一般形式： 掣+ l ( ”) ：m ) a t 、。 其中l ( u ) 为甜的线性项，f ( u ) 为甜的非线性项。于是，时间上采用的差分 格式如下： 2 r o l o ( u 斛1 ) + a t l ( u 肘1 ) = 口，厶 纠) + a t f l j ( u ”7 ) 】 ( 2 1 4 ) i = 0 系数选取如下： 表2 3 1 外推法的参数选择 口28 。屈及 一阶格式1lo0100 二阶格式 3 221 2021 0 三阶格式 1 1 633 21 3331 其中，缸表示时间步长。 在求解的时候，先使用小步长，如a t 1 0 ，第一步使用一阶格式计算，第二 步使用二阶格式计算，第三步及其以后的各步使用三阶格式计算，直至计算出正 常时间步长所需要的前两个时间层上的解，然后再使用正常时间步长及三阶格式 进行时间上的计算。 这样，将空间上和时间上离散格式都代入到方程( 2 6 ) 中，利用展开级数的正 交性，可以得到方程在配置点上的离散为： ( 4 ，一最一c 卅) + 1 ( 胛，m ，f ) = 见p ( 少，m ，f ) h + 鹕4 盆( y ，m ，矿一) ( 2 1 5 ) p = o 1 = 0 其中 4 ，= ( t o + 姚f ) ( 9 2 d 2 + 昭d ) 乙( s ) 吃= ( + f 必f u o ) ( 口2 + m 2 2 ) 瓦( s ) q = 2 i a a t u o ，g d 五( j ) 见= 9 2 d 2 + g g 臼一( 口2 + 聊2 2 ) 。 第二章控制方程与计算方法 2 4边界条件 因为给定流向波数，问题可以归结为一个二维问题。因此，计算域选取为： ( y ，z ) ( 0 ，+ ) ( 一t 2 ，t 2 ) 在展向上采用周期性边界条件：p ( y ，z + t ，f ) = p ( y ，z ，f ) ，其中t 为展向的 基本波的波长。 由于在壁面和法向无穷远处，法向速度，为零，因此对于压力边界条件有： 在尸0 处， 在无穷远处， 动n 一2 u 砂 至：o ， 砂 即孚：o o s 即孚：o 倦 将空间上和时间上离散格式代入上式，在y = o 处和无穷远处分别为： n vn y 挖2 声( 刀，m ，f ) = o ，( 一1 ) ”刀2 p ( 刀，m ，f ) = o ( 2 1 6 ) n = on - - o 至此，边界条件( 2 1 6 ) ，加上方程( 2 1 5 ) 在，一1 个配置点上提供的m - 1 个方 程，构成脚1 个方程，求解该方程组，可以求得p ( 胛，m ，f ) 。这样就可以通过它 求得压力p 在弦平面上的分布。 2 5速度与压力的关系 由前面推导的离散方程，求解可以得到压力p 的分布，但同时，我们还希望 知道速度的分布状况，因此，就需要推导速度求解的方程以及对其进行离散。 根据流动稳定性理论，流动的瞬时量可以分解成基本量和扰动量，扰动量则 可以使用行进波的形式进行表示： ，= 缈( y ，z ) e k a 耐( 2 1 7 ) 其中，对于时间模式，口和都为实数，c o 为复数。我们又可以将缈写成 c o = 戤+ i o j t ，这里蛾表示波传播的频率，劬表示波的增长率，其为正时，表示 波会随时间增长，波是不稳定性的，反之，则波会衰减，是稳定波。 将速度与压力都按如上形式展开，有： 材u ( y ，z ) d “一耐 ，v ( y ，z ) “一科 w ，= w ( y ，z ) d “一耐 p p c y ，z ) e 。埘 将( 2 1 8 ) 代入到控制方程( 2 1 ) 中，得n - 1 2 ( 2 1 8 ) 第_ 章控制方程与计算方法 整理，得到： 磁”+ 竺+ 坐：0 砂 瑟 一砌+ 妇沈+ v 型+ w 型：一圾p ，、a y 良 1 ( 2 1 9 ) 一i r o v + 妇u v ：一皇巳 砂 一泐w + 娩u w ：一皇巳 印 ，：f 塑 一缈+ a u ( 2 2 0 ) 望 w = f 监 ( 2 2 1 ) 其中，国可由上面求解p 的时候得到的增长率以及圆频率得到，而p 沿展向 的导数可以表示为： 鱼：n , 2i m f l b ( y 妙 ( 2 2 2 ) 0 z ，身，2 7、7 求得v ，w 后再由流向的动量方程，可以给出甜的分布： 0u a u - t a p v i 一一w - = - 甜：皇，_ 堕 ( 2 2 3 ) 1 3 第三章不可压缩流体边界层条纹结构无粘稳定性的数值研究 第三章不可压缩流体边界层条纹结构无粘稳定性的数值研究 3 1本文所使用的基本流 本文计算所使用的基本流是由吴雪松教授提供的。该流场是考虑外界湍流 度，通过求解一个非线性的边界区方程得到的。本文主要是对该流场使用时间模 式进行稳定性分析，因此所使用的流场为在某一特定流向位置x 处的流场。因为 流向条纹结构沿展向排布，故可以设基本流的表达式为： u ( y ，z ) = v o ( y ) + 4 u ( y ，z ) ，其中u t y ，z ) 表示基本流场，u o ( y ) 表示b l a s s i u s 解， u ( ) ，z ) 表示由条纹结构所引起的基本流的变化，彳。则表示条纹幅值。图3 1 1 给出了基本流在y ，z 平面上的等值线。我们注意到该流场结构在展向上是关于 z = o 的位置轴对称的。图3 1 2 和图3 1 3 分别给出了速度沿法向和展向的导数 等值线图。我们可以注意到，在壁面附近速度梯度很大。同时，图3 1 - 4 图3 1 5 、 图3 1 6 中给出了速度的二次导数等值线图，可以看到速度剖面中具有拐点。 z 图3 1 - 1 平均流洮z ) 在弘平面的等值线图 图3 1 2 _ a u 在弦平面上的等值线图 d y 1 4 第三章不可压缩流体边界层条纹结构无粘稳定性的数值研究 图3 1 3 掣在弦平面上的等值线图 必 图3 1 - 4 等在弦平面上的等值线图 图3 1 - 5 害在归平面上的等值线图 第二章不可压缩流体边界层条纹结构无粘稳定性的数值研究 图3 1 6 霎娶在弘平面上的等势线图 c s , o z 对于本文中所计算的流场，其基本流表达式对应于4 = 1 的情况。其基本参 数如下所示： - 静2 序m j 叹堋，例0 1 ， 对于基本流的r e y n o l d s 数定义为r e = ( 2 r r h q ) 2 3 2计算参数的选取及其验证 在本文中，需要给定流向波数。这里准备计算流向波数口从0 1 到0 8 的值。 计算总时间长度不定，一直计算到增长波呈指数增长状况为止。 计算时，y 和z 方向上的计算区域为： ( y ，z ) ( 0 ，佃) ( 屯：2 ，t 2 ) ， 其中厶= 2 n 。选取的网格为：在y 方向，1 ，= 1 2 8 个网格；在z 方向， - = 1 6 ；时间步长f = o 0 2 。 为了验证以上计算参数的选取，计算格式的选择是否得当。本文采用如下的 方式进行验证：使用文中的空间和时间离散方法数值模拟计算o s 方程，将其求 得的结果与线性稳定性理论的结果进行对比，以此验证方法的可行性。 图3 2 1 为对b l a s s i u s 解为基本流时，不稳定波的幅值沿时间变化曲线与线 性稳定性理论的对比，图3 2 2 则给出了数值模拟求得的特征函数的分布与线性 稳定性理论的比较。计算所选参数为r e = 8 0 0 ，口= 0 2 2 1 2 1 6 ，= 0 。n s 表示 数值模拟的结果，l s t 表示线行稳定性理论的结果。 1 6 第三章不可压缩流体边界层条纹结构无粘稳定性的数值研究 t 图3 2 1 扰动幅值随时间变化的对比 1 0 8 0 6 e 兰 o j o 工 o y 图3 2 - 2 特征函数分布 得到的结果重合的很好，这在某种程度上证明本文的计算方法还是可信的。 本文还对法向的网格加密一倍的结果进行了对比。图3 2 3 给出了初始扰动 幅值为o 0 1 ，t = 1 0 0 0 的特征函数比较图，从图中我们可以看到，结果是完全重 合。因此我们可以认为法向上网格的选取是合适的。 y 图3 2 3 不同法向网格数下计算的特征函数对比 为了得到展向有条纹结构流场的不稳定扰动波的解，我们可以在初始流场中 1 7 第三章不可压缩流体边界层条纹结构无粘稳定性的数值研究 加入任意初始扰动，随着时间演化，稳定的特征波将逐渐衰减，趋于零，不稳定 的特征波将随时间增长。如果对于给定的参数，只有一个不稳定的特征波，则经 过长时间的演化，流场中的主要成份将剩下该不稳定的特征波，且以指数形式增 长。这样，我们就可以得到该不稳定波的各种波参数和扰动沿法向的分布。 为了方便，本文在第一次计算时，采用与基本流为b l a s s i u s 解的o s 方程得 到的特征函数线性作为流场的初始扰动，即取： 3 ( y ，0 ，0 ) = a j b ( y ) 其中， 3 ( y ) 为由b l a s s i u s 解为基本流的o s 方程得到的特征压力，a 为扰 动的幅值，本文中口取l o 。 在以后的计算中，我们只需要用在该参数附近得到的特征函数作为该参数的 初始条件即可，这样可以更快的求得该该参数下的特征函数，可减少计算时间。 第四章条纹边界层结构的稳定性计算结果与分析 第四章条纹边界层结构的稳定性计算结果与分析 事实上我们注意到，基本流场沿展向和法向有变化，而沿流向无变化，同时 计算方程使用的是线性方程，因此只需要计算在不同流向波数下的情况即可。 本文分两方面给出计算结果和分析： ( 1 ) 同一流向波数下的计算结果与分析 ( 2 ) 不同流向波数下之间的比较 4 1口= 0 4 时的情况 4 1 1 波的增长率 我们将特征函数p ( y ，z ) 沿展向展开，写成如下形式： p ( y ，z ) = 多( y ) 纱励 ( 4 1 ) m = - 7 计算流向波数口= 0 4 ，入上所述，只在展向谱中m = 0 的位置加入给定初值， 初始幅值为1 0 。经过一段很长时间的计算( f 1 0 4 ) ，我们得到如下的m = 1 波的 幅值的随时间的增长变化： t 图4 1 1 扰动p 沿法向最大幅值随时间变化( 朋= 1 ) 将其取自然对数，