（流体力学专业论文）带壁面抽吸的槽道湍流大涡模拟.pdf

中文摘要 本文采用大涡模拟的方法，对壁面加入周期性抽吸边界条件的槽道湍流进行 了数值模拟。采用高精度紧致有限差分格式与f o u r i e r 谱展开相结合的计算方 法，使用三阶混合显隐格式时间分裂法，求解大尺度量满足的n a v i e r s t o k e s 方 程，小尺度量分别采用四种亚格子应力模型进行模拟，并在涡粘性模型中引入了 壁面函数，得到了雷诺应力、平均速度变形率的空间分布与相位分布，其结果与 直接数值模拟的结果一致。同时讨论了不同亚格子应力模型对模拟结果的影响。 关键词：大涡模拟，亚格子应力模型，紧致有限差分，抽吸 a b s t r a c t l a r g e e d d ys i m u l a t i o nw a su s e dt os i m u l a t et h et u r b u l e n tc h a n n e lf l o ww i t h p e r i o ds u c t i o no nt h ew a l l h i 曲o r d e rc o m p a c tf i n i t ed i f f e r e n c es c h e m ea n df o u r i e r t r a n s f o r m a t i o nw e r ea l s oa d o p t e di no u rm e t h o d t h et h i r do r d e rm i x e de x p l i c i t i m p l i c i t s c h e m ew a se m p l o y e df o rt h et i m ed i s e r e t i z a t i o n l a r g es c a l es t r u c t u r ew a sc a l c u l a t e d f r o mt h en a v i e r - s t o k e se q u a t i o n 1 1 1 es m a l ls c a l es t r u c t u r ew a ss i m u l a t e ds e p a r a t e l yb y f o u rs u b - g r i d s t r e s sm o d e l si nt h i sp a p e r , w h i l et h eb o u n d a r yf u n c t i o nw a si n t r o d u c e d i ns m a g o r i n s k ym o d e t h ed i s t r i b u t i o no fr e y n o l d ss t r e s s ，t h em e a nr a t eo fs t r a i n a n dt h ep h a s ed i s t r i b u t i o no fr e y n o l d ss t r e s sa n dt h em e a nr a t eo fs t r a i nw e r e o b t a i n e d t h e s er e s u l t sw e r ec o i n c i d ew i 山t h o s eo fd i r e e tn u m e n c a ls i m u l a t i o n m e a n w h i l e t h ee f f e c t so ns i m u l a t i o no fd i f f e r e n ts u b g r i d s t r e s sm o d e l sw e r e d i s c u s s e di nd e t a i l k e yw o r d s ：l a r g ee d d ys i m u l a t i o n ，s u b - g r i d - s t r e s sm o d e l ，c o m p a c tf i n i t e d i f f e r e n c es c h e m e ，s u c t i o n 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的 研究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表 或撰写过的研究成果，也不包含为获得叁洼盘堂或其他教育机构的学位或证 书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中 作了明确的说明并表示了谢意。 学位论文作者签名：新大为 签字日期：2 。口5年月，日 学位论文版权使用授权书 本学位论文作者完全了解鑫鲞盘堂有关保留、使用学位论文的规定。 特授权苤注盘堂可以将学位论文的全部或部分内容编入有关数据库进行检 索，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校 向国家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名：郄大为 导师签名： 方一红 签字日期：2 。5 年月，目签字目期：加d 年f 月 f 日 第一章绪论 第一章绪论 自然界存在着层流和湍流两种不同的流动形态，最常见的流动是湍流。湍 流研究上的进展将推动有关技术的进步，推动国民经济与国防建设的发展。同时， 湍流是一种非线性的复杂问题，对其研究将加深人们对广泛存在于自然界和社会 各个领域中的非线性复杂现象的认识。 人们对湍流的研究已有i 0 0 多年的历史。虽然经过包括许多伟大的学者在内 的无数优秀的科学家顽强的努力，对湍流研究取得了相当成果，并且将这些理论 成果与实践相结合解决了很多工程实际问题。但是由于湍流运动的复杂性，其机 理至今仍未被人类完全掌握，湍流也就成为经典物理留下的为数不多的几个难题 之一。因其对科学技术发展的重要性，世界上各个科学大国仍然把湍流研究列为 需要优先发展的重大基础研究课题之。 1 。1 湍流的数值模拟 近代电子计算机规模与速度的飞速发展给人们提供了一种解决湍流问题的 新途径。一般认为包括脉动运动在内的湍流瞬时运动也服从n a v i e r s t o c k s 方 程，而n a v i e r s t o c k s 方程本身是封闭的，不需要另外建立模型。由此提出一种 想法，是否可以不引入任何湍流模型，而用电子计算机数值求解完整的非定常的 n a y i e r s t o c k s 方程，对湍流的瞬时运动进行直接数值模拟，得到各种统计平均 量。这样保证了方程本身是精确的，仅有的误差也只是由数值方法引入的，并且 可以提供每一瞬时流场上所有流动量的全部信息。特别有意义的是此方法能提供 很多在实验上目前还无法测量的量，这就意味着可以提供基础数据，可以利用直 接数值模拟的结果来检验各种湍流模型，发展新的湍流模型。由于数值模拟中流 动条件能够得到精确控制，可以对各种因素单独或交互作用的影响进行系统的研 究，所以成为一种新的湍流预测手段。 湍流数值模拟的方法主要有两种：大涡模拟( l e s ) 和直接数值模拟( d n s ) 。 最早是由0 r s z a g 等人于上个世纪7 0 年代开始从事湍流的直接数值模拟。他 们对谱方法的发展作出了重要贡献，并用此方法计算了湍流转捩问题。k i m 、m o i n 和m o s e r ( 1 9 8 7 ) 用1 9 2 x 1 2 9 1 6 0 个网格点直接模拟了r e = 3 3 0 0 的槽道湍流， 建成k m m 数据库。j i m e n e z ( 1 9 8 7 ) 在k m m 数据库的基础上向人们展示了展向涡 的分布及其分裂和增长，a d r i a n ( 1 9 8 7 ) 再现了流向涡的喷射现象，g u e z e n n e c ( 1 9 8 9 ) 通过象限分裂法分析了k m m 数据库，发现在瞬时流场中近壁区湍流相干 第一章绪论 结构的流向涡大部分并非是反向旋转的涡对，而更多的是强度和尺度不对称的流 向涡，有时甚至只有一个流向涡。r o g e r s 和m o i n ( 1 9 8 7 ) 用1 2 8 1 2 8 1 2 8 个网 格点模拟了均匀剪切湍流，形象地描绘了发卡涡的结构。s p a l a r t ( 1 9 8 8 ) 用 4 3 2 8 0 3 2 0 个网格点模拟了雷诺数r e = 1 4 1 0 的平板湍流边界层，这一数据库 被广泛地用来研究湍流边界层的性质。 1 2 湍流的大涡模拟 大涡模拟最早是由气象学家s a m g o r i n s k y 于1 9 6 3 年提出。大涡模拟是介于 一般模式理论与直接数值模拟( d n s ) 之间的一种方法。其基本思想是把包括脉 动在内的湍流瞬时流动通过某种滤波的方法分解为大尺度运动和小尺度运动两 部分。大尺度运动通过数值求解运动微分方程直接计算出来；小尺度运动对大尺 度运动的影响将在运动方程中表现为类似于雷诺应力一样的应力项，称之为亚格 子雷诺应力。它们将通过建立模型来模拟。大尺度结构强烈依赖于流场的边界条 件和初始条件，并且与平均流之间有强烈的相互作用，可直接从平均流中获取能 量，反之平均流也受到大尺度结构的强烈影响。一般认为小尺度结构是随机的， 其与边界条件、初始条件无关。d e a r d o r f f ( 1 9 7 3 ) 首次把大涡模拟的方法应用 于具有工程意义的槽道流动的数值模拟，并为这一方法奠定了基础。之后， s c h u m e n n ( 1 9 7 3 ) 和g r o t z b a c h ( 1 9 7 7 ) 对他们的工作进行了改进，并推广到较 复杂的流动中。s t a n f o r d 大学的f e r z i g e r 和r e y n o l d s 等人从1 9 7 2 年起对大涡 模拟进行了系统研究：l e s l i e 等人从1 9 7 5 年开始对亚格子模型进行了系统研究； m o i n 和k i m ( 1 9 8 2 ) 第一个用6 3 6 4 2 8 个网格点模拟了槽道湍流；g e r m a n o 和 p i o m e l l i ( 1 9 9 1 ) 提出了一种新的亚格子应力涡粘性模型，取得与实验符合得很 好的结果。 大涡模拟的第一步就是要将一切流动变量划分成大尺度量与小尺度量。这一 过程称之为滤波。如爿f x ，t 1 是任意一个瞬时的流动变量，则其大尺度量可以通过 以下一个在物理空间区域上的加权积分来表示 j ( 列) = g ( i x - x 1 ) a ( x , t ) d v ( 1 1 ) 其中权函数g ( k z i ) 即称为滤波函数。 这样瞬时量与大尺度量之差 彳7 = 彳一一( 1 2 ) 第一章绪论 反映了小尺度运动对a 的贡献，称为a 的亚格子分量，或小尺度量。 1 3 大涡模拟的滤波函数 大涡模拟采用的滤波函数主要有三种：b o x 滤波函数，傅氏滤波函数以及高 斯滤波函数。 13 1d e a r d o r f f 的b o x 滤波 d e a r d o r f f 定义滤波函数g ( 一y ，z ) 为： i = 1 , 2 , 3 ，“5 - ，( 1 3 ) 江1 。72 。73 五l - 缸- 4 - 式中，一为任意网格点坐标，她为第i 方向的网格尺度。大尺度量j 实际上就是 在以t 为中心的长方体单元( b o x ) 上的体积平均值，故这种滤波方法也称为b o x 方法。这种方法很简单，缺点是它的傅立叶变换在某些区间里有负值，并且由于 滤波函数在单元边界上的间断性，难以进行微分运算。 132 傅氏截断滤波器 傅氏滤波器实际上是b o x 滤波器在谱空间的翻版。即在傅立叶展开式中简单 地截去所有波数绝对值高于k 。的分量。 郎) = i k i b k 岛。 而它在物理空间对应的滤波函数： g = 击弘) e - “矾 同样有在某些区间内有负值和难于求微分的缺陷。 ( 1 4 ) ( 1 5 ) 忐。 ，cl f f 、lj x x ，l g 第一章绪论 反映了小尺度运动列a 的贡献，称为a 的亚格子分量，或小尺度量。 1 3 大涡模拟的滤波函数 大涡模捌采用的滤波函数主要有三种：b o x 滤波匝数，傅氏滤波函数以及嘉 斯滤波函数。 13 1d e a r d o r f f 的b o x 滤波 d e a r d o r f f 定义滤波函数g ( x ，y ，z 1 为 i = 1 , 2 , 3 ，i x - x i 一z - ， 。5 1 。7 2 。3 睁一j t a l i 式申，一为任意稠格点坐标，t 为第i 方向的网格尺度。太尺度量j 实际上就是 在以t 为中心的长方体单元( b o x ) 上的体积平均值，故这种滤波方法也称为b o x 方法。这种方法很简单，缺点是它的傅立叶变换在某些区间里有负值，并且由于 滤波函数在单元边界上的间断性，难卧进行微分运算。 3 2 傅氏截断滤波器 傅氏滤波精实际上是b o x 滤波器在谱空间的翻版。即在傅立叶展开式中简单 地截去所有皱数绝对值高于k o 的分量。 4 ( k ，= ；乏 而它在物理空问对应的滤波函数： 1岬 g ”壶g 【引产掀 同样有在某些区n u 内有负值和难于求微分的敢陷。 同样有在某些区间内有负值和难于求微分的缺陷。 ( 1 4 ) ( 1 5 ) 去。 ，。i f f r r l g 第一章绪论 1 3 3 高斯型滤波器 高斯型滤波器滤波函数为 或 珥3l 矛6 1 ，2 e 冲 e f 1 ：e 一。2 内2 4 6 ( t t 。) 2 ( 1 6 ) ( 1 ，7 ) 高斯型滤波函数的傅氏变化也是高斯型函数，在物理空间和谱空间都有很好的性 能，可以微分任意次。虽然高斯型滤波器性能最好，但计算很麻烦，目前用得最 多的还是前两种滤波器，因为它们简单方便，易于实施。 1 4 亚格子模型 其中 将上述的滤波方法运用到n s 方程，即得方程 冬+ 当面：罢v - + ，要+ 挈 ( 1 8 ) 言+ i 吩哆一瓦w 菇+ 茜 r 。= 一只，+ ；毛r 。 r = 珥“j + “羁+ “；“j ( 1 9 ) ( 1 - 1 0 ) 毛称为亚格子雷诺应力，通过建立亚格子应力模型来模拟。 目前常用的亚格子应力模型有涡粘性模型( s m a g o r i n s k ym o d e l ) 、相似性模 型( s i m i l a r i t ym o d e l ) 、动态模型( d y n a m i cm o d e l ) 、结构函数模型、混合模 型( m i xm o d e l ) 、高阶模型等。 1 41 涡粘性模型 涡粘性模型基于“产生项等于耗散项”的假定推导出来，它假定亚格子尺度 第一章绪论 湍流应力与大尺度分量的应变率成正比， 一q 鬲叫f 摹+ 割 其中u ，称为涡粘性系数，具有和运动粘性系数同样的量纲： 坼= ( e ) 2 同 ( 1 1 2 ) 其中c 。是无量纲模型系数，在标准的s m a g o r i n s k y 模型中c 。是一常数：滤波宽 度是最小可求解涡的长度尺度。 涡粘性模型对于简单剪切湍流的模拟很成功，在某些均匀湍流的亚格子尺度 模拟中也同样有效。对于没有平均应变的均匀湍流，可以用涡粘模型预测大多数 的低阶统计量。对于具有应变和剪切的均匀湍流，能量反向迁移很明显，它们从 小尺度涡传到大尺度涡，此时亚格子尺度湍流不总是消耗大尺度涡的能量，因而 建立在小涡从大涡消耗能量基础上的涡粘性模型难以预测。另外，使用涡粘模型 时，必需假定模型常数c 。，该系数和流动的形态、网格分辨率、滤波宽度等因 素有关，做为常数处理，明显不太台理。 f y is d r j c hn i e u w s t a d t ( 1 9 9 4 ) 发现，由于涡粘性模型在预测近壁区的能 量传输机制时的不稳定性而使其在与实验数据比较时存在问题。特别是当涡粘性 模型应用到流动是由一个小扰动叠加到层流上开始的平板流的边界层问题时，由 于来自于平均剪切作用的过分涡粘性而使流体始终保持为层流状态。由此可见， 涡粘性模型的应用在这种情况下是很不成功的。 14 2 动态模型 动态模型( d y n a m i cm o d e l ) 由o e r m a n o 等人于1 9 9 1 年提出，设计此模型在 模拟过程中求解尺度大小并测量模型相关参数。这样就避免了预先设置或调试模 型参数。定义模型参数： 式 ( 吖= 踹 ( 1 1 3 ) 此模型是建立在涡粘性模型基础上。其动态过程基于g e r m a n o ( 1 9 9 2 ) 恒等 第一章绪论 f ( x ，f ) = 铲( x ，f ) r ；( z ，f ) ( 1 1 4 ) 此时滤波尺度为以，一般认为a = 2 ( g h o s a l ，1 9 9 5 ) 。 小( 务i 虱 ( 1 1 5 ) 在考虑到尺度不变( s c a l ei n v a r i a n c e ) 的假设下，可令c ，= c ，这样，就 可以由模拟过程完全决定。 动态模型的系数c ? 在空间的变化很大，且含有相当一部分负值，表现出了 某些数学j 二的不连续性。但允许出现负值是这个模型的一个优点，因为这些负值 表示了在物理空间的一种方向散射。然而过大的负涡粘系数会使数值模拟过程不 稳定，且常常会出现可求解尺度能量的非物理性增加。通常采用的补救措旅是在 时间上和空间上对上式分子和分母都进行平均，但这样做会丧失动态局部方程所 具有的一些概念上的优点。折中的方法是在流体的均匀方向上进行平均。 1 4 3 结构函数模型 结构函数模型基本思想是试图将湍流的局部间歇性考虑到计算中，并减小因 为小尺度涡没有充分发展而对涡粘性的过高估计产生的误差。其表达式为 州础，：掣r 结构模型能够很好的处理自由剪切流、分离流动，如g o n z e ( 1 9 9 3 ) 对卡门 涡街的空间演化和s i i v e i r a n e r o 等后台阶绕流结构的模拟都表明了这一点。 然而，与s m a g o r i n s k y 模型一样，结构函数模型在处理低马赫数的边界层分 离问题上也存在不足。尽管采用平行于壁面的四点公式能够消除壁面上的平均剪 切对涡粘性的影响，但结构函数模型在处理槽道流问题时表现也并非很好。为了 克服这些困难，结构函数模型进一步改进为选择性结构函数模型( s s f ) 和过滤 性结构函数模型( f s f ) 。 1 44 相似性模型和混合模型 相似性模型由b a r d i n a 等人于1 9 8 0 年第一次引入。这个模型强调并基于尺 第章绪论 度不变( s c a l ei n v a r i a n c e ) 这一原则，认为流场速度场中，尺度在以下的结 构与尺度在以上的结构是相似的。这一假设是基于l i u ( 1 9 9 4 ) 得出的 b a n d p a s s f i i t e r e dp i v ( 带通滤波粒子图像测速技术) 测量理论，其指出在 系列连续尺度做出的速度场分析图上，不同尺度的一些特定的结构几乎总是相似 地出现在同一个区域。这样r ? 也必须是个与应力张量相似的量： 矿“c s 。，( 丽一霉霉) ( 11 7 ) 这里就需要引入二次滤波，滤波尺度芦，其中y 1 ，c 。取决于7 。b a r d i n a 建 议取y = 1 ，而l i u ( 1 9 9 4 ) 等人建议取y = 2 ，a k h a v a n ( 1 9 9 9 ) 等建议取y = 4 3 。 但选择一个普适的二次滤波尺度，以及滤波尺度与数值方法之间的关系等一系列 问题，至今仍然没有解决。 根据b a r d i n a ( 1 9 8 0 ) 的结果，采用高斯滤波器或b o x 滤波器，模型的应力 与真实应力相似度可达8 0 ；同时相似性模型也表现出了能量的反向散射特性， 但在仅仅使用相似性模型的模拟中并未能耗散足够多的能量，而使最后的结果失 真。针对于此，b a r d i n a ( 1 9 8 0 ) 建议加入s m a g o r i n s k y 耗散项，得到混合模型： r 护= c ，。( 丽一茸亏) 一2 ( c ：q 2 同瓦 ( 11 8 ) s p e z i a l e ( 1 9 8 8 ) 等人进一步将这类模型推广到了可压缩流动的大涡模拟中。 1 4 5 非线性或梯度模型 相似性模型采用了二次滤波，增t n y 进行第二次滤波的计算支出。解决方法 就是将第一次滤波的结果i ，用t a y l e r 级数展开并进行滤波分析( l e o n a r d1 9 7 4 ， c l a r k1 9 7 9 ，l i u1 9 9 4 ) 。这样得到的就是非线性、梯度的，或涡粘性张量的模 型： e 。( 丽一霉瓦) 兰铋2 詈要 ，。， 这里c ，也取决于y 和滤波方法。 1 5 本文的工作 本文研究对象是具有空间周期性的槽道流动，采用大涡模拟的方法，控制方 第一章绪论 程为滤波后的n s 方程，流向和法向采用高精度紧致差分格式离散，展向采用 f o u r i e r 展开。主要工作如下： ( 1 ) 建立大涡模拟的程序，对槽道湍流进行数值模拟。通过对湍流壁面 率、雷诺应力与平均速度变形率、脉动速度均方根值等湍流重要特 征量的计算，对程序进行验证。 ( 2 )采用经典的s m a g o r i n s k y 模型，引入壁面函数，对下壁面加入周期 性抽吸边界条件的槽道湍流进行大涡模拟，得到雷诺应力与速度变 形率的空间分布、相位分布、相位差等，并与d n s 的结果比较。 ( 4 ) 分别使用相似性模型、混合型模型，对下壁面加入周期性抽吸边界 条件的槽道流场进行大涡模拟，得到雷诺应力与变形率的空间分布、 相位分布以及相位差等结果，并与d n s 的结果比较。 在本文中，我们还建立一种新的亚格子模型，即由d n s 数据中提取 出来的f 模型，用该模型对上述同一问题进行了大涡模拟。并将上 述四种亚格子模型的结果进行对比分析。 第二章数值计算方法与边界条件处理 第二章数值计算方法与边界条件处理 本文研究的是不可压槽道湍流，一般考虑简单的矩形计算域，即设计算域为 0 x x 。，0 y 儿，0 z z 。，如图2 - 1 所示： 定义流向波数口= 0 3 ，展向波数口= 0 4 ，法向长度为2 。流向、展向为周期性 边界条件，其中x 、y 方向采用差分格式进行离散，z 方向使用f o u r i e r 谱展开。 2 1 控制方程 对完全n s 方程： 进行滤波，得到 其中 堕+ 生型：一土鱼+ ，旦旦。( 2 1 1 ) o t 3 x s p 魄 o x jo x j 鲁= 一毒( 画) 一吉善+ 丙2 瓦 c z ，z ， 坼0 = ( 巧+ u ) ( u - j + 形) = 呸巧+ 霹甜j + “巧+ 甜彰 ( 2 1 3 ) 等号右端后三项包含小尺度量，必须建立模型。把这三项之和称为亚格子雷诺应 力 r f = 酗j + “+ “杉 ( 2 1 4 ) 第二章数值计算方法与边界条件处理 和 其中 通常把亚格子雷诺应力张量分解成一个对角线张量与一个迹为零的张量之 嘞= ( 铲三毛如 + l f i u r , 。= 1 + ；岛 c z ， r 。= 一尺。+ ；岛r 。 将对角线张量部分与压力项合并，可定义一个修正的压力 于是滤波后的n s 方程可写成 p = 旦+ ；r 。 p 3 ( 2 1 6 ) ( 2 1 7 ) 亟+ 旦瓦：一堡+ ，堕+ 堕( 2 1 8 ) o t 瓠i j i 瓠t a x i 瓠io x i 其中通过建立模型来模拟。 本文使用b o x 滤波函数，则有如下结论 于是有 亚格子雷诺应力即为 再根据连续性方程 嚣= 瓦和u ，= 0 u i l d j2 u i u j 瓦“j = “羁= 0 r = “鼻； ( 2 ，1 9 ) 第二章数值计算方法与边界条件处理 式( 2 1 8 ) 可化为 堕：o 优 亟瓴亟：一堡+ u 鱼+ 堕 a i a ) i 0 x t 融i 融i 钕i 当：o 血 2 1 1s m a g o r i n s k y 模型 首先采用s m a g o r i n s k y 涡粘性模型： t i = 2 u t s 定义涡粘性系数： q = c c ：，2 l 亨| = c c ，2 c z 品岛) “2 = c c ，2 去 鲁詈+ 鲁 2 “2 其中s m a g o r i n s k y 常数e = o 1 ，滤波宽度= ( 越2 + 缈2 + z 2 ) “2 ，缸、缈和止 分别是在z 、y 及：方向上的网格宽度。 令： 髟一4 2 ( ca ) 2 ， ， 有： l l = k 经典的涡粘性模型一般无法给出准确的近壁区湍流行为，使用时可在壁面附 近进行处理。根据v a nd r i e s t ( 1 9 5 6 ) 的指数粘滞函数，将s m a g o r i n s k y 模型 、 堕魄 + 甄一 ，l q 第二章数值计算方法与边界条件处理 的涡粘系数进行修正，改写成如下形式： 其中d 为近壁衰减函数 其中y + 表示近壁无量纲距离 q = ( c d ) 2 l i i d ：h 冲卜刊 v + ：型一+ ：2 5 u 这样得到新的涡粘性系数： 吲删l 曩i = m ：d a ) ( z 刚1 1 2 , c d a , 2 ( ir 虿o 磊i + 驯” 式( 2 1 1 0 ) 化为 k ：拿( c d ) 2 鲁嗳=一詈十u熹soxs4-+ 毒卜( 善+ 詈 一“，二= 一十u l + ln l 二+ 二l 乱 舐， 反。 舐 缸，【2i 出，晚j 盟：0 融 对上式进行无量纲化 z ，：三，f ，：上，玩：旦，尸，：p - - 5 ，址：生 h h i u m u mu m 2 y 其中h 为半槽宽，v 。为来流最大平均速度。 为方便记述，略去符号上面一横，得到下式 第二章数值计算方法与边界条件处理 鲁m 考一面8 p + 上r e 器+ 土r e 专卜( 考+ 等) 二+ “上= 一+ l + 一lu l 二+ 二i 乱 。玉，魄缸，瓠，缸，l7 【出，钆jj 堕：0 屯 21 2 相似| 生模型 其次我们选用了相似性模型 砖。“= e 。( 丽一霉亏) ( 2 1 1 5 ) ( 2 1 1 6 ) 相似性模型需要进行二次滤波，在本文中采用滤波宽度为o r ，取a = 2 ，其中 为一次滤波宽度： = ( 血2 + 缈2 + 垃2 ) “2 本文取相似性系数：c 。= 1 0 ，将式( 2 1 1 0 ) 化为 | 吾+ 弓詈= 一善+ u 去+ 毒( e 。( 丽一嘱) )f 言柏，茜一百蜘石蠢+ 瓦m 悸吩叫r 叶j j 里：o 【融 对上式进行无量纲化 zr ：兰，r 一：上， k h 。扯篑2 西p 为方便记述，略去符号上面一横，得到下式 鲁托考2 一瓦a p o r e 杀+ 毒( ( 万一砭巧) ) a f。曲，出，苏氖舐“” 7 ， 挚：o 【呸 1 3 ( 2 1 ，1 7 ) ( 2 1 1 8 ) 第二章数值计算方法与边界条件处理 2 13 混合型模型 再者我们选用了混合型模型，它是相似性模型与涡粘性模型的混合。 定义： r = e 。( 雨一茸弓) 一2 ( c ) 2 同墨 ( 2 1 1 9 ) 本文中取s m a g o r i n s k y 常数口= o 1 ，相似性系数e 。= 0 9 ，二次滤波宽度。= 2 。 式( 2 1 1 0 ) 化为： 鲁+ 唔一詈+ u 器+ 毒卜。( 丽一羁) 一q l 善+ 割j 8 t j 8 x j a x 。 瓠j a x i 蚤x i “u 。 。 j 、瓠i 融。) 3 堕：o 觑 ( 2 1 2 0 ) 取同上无量纲化参数得到下式 吉3 u m 蠢3 u = 一筹+ 上r e 杂+ 毒卜。( 瓦一瓦乃) 一去【考+ 警j j a 。孤，舐，a o a 一瓠，l “叭r 。l 舐，缸。jj 丝：o 巩 ( 2 1 2 1 ) 21 4f 模型 在本文研究中，我们还采用了一种新的模型来模拟亚格子雷诺应力。利用 d n s 数据库( 6 4 1 2 8 x 6 4 ) ，对速度场u t 进行b o x 滤波( 3 2 3 2 3 2 ) 。其中流向、 展向使用了s i m p s o n 数值积分公式，法向使用了c o t e s 积分公式。得到大尺度量 玩和小尺度量“0 力 对“进行二次滤波，滤波宽度以，其中口= 1 。由定义得到亚格子雷诺应 进而得到了f 模型 吩= “4 j 第二章数值计算方法与边界条件处理 铲瑚一+ 孝毛如 同样进行无量纲化，得到控制方程： 盟机堕：一堡+ 上! 坠+ 盟 。l j a x i 融t r e 融j a x i 融j 堕：o o x ： 2 ，2 数值方法与边界条件处理 本文在流向、法向均采用差分格式离散，展向使用傅立叶谱展开。流场中的 物理量e ( x ，y ，z ，f ) 可表示为如下形式： 旦一1 2 庐( x ，y ，z ，t ) = 九( x , y ，t ) e “肚 ( 2 2 1 ) 。十生 2 方程中各个物理量展开形式如下 一n l 2 u ( x ，y t ) = ( x 川y ) p 1 邮2 1 一兰 2 旦一l 2 w ( x ，y ，z ，t ) = w m ( x ，y ，t ) e 肛 d 一旦 2 艺( 训，f ) e 一肛 生1 2 p ( x ，舻，t ) = 己( z 川y ) p 1 郇2 m ：i _ 旦 2 其中展向波数口= 2 r e z l ( z l 是展向长度) 。 非线性项佰v ) a 的f o u r i e r 展开可写为： ( 厅v ) 厅= f 陋甲) 厅) f 蛳 ( 2 2 2 ) 其中，f m 何v 弦 称为非线性项( i v ) a 的f o u r i e r 谱展开的系数。 这样方程就简化为我们下边要经常讨论到的准二维问题的方程组 第二章数值计算方法与边界条件处理 其中 鲁+ 蹦厅v ) 啦一v 。+ 面1 v 2 - - + v 。 v 。厅。= 0 ( 2 2 3 ) v 。= ( 昙，昙，一i m , b ) 0 洲 v 等+ 等一内 f o u r i e r 谱展开后，z 方向的周期性边界条件自然满足，x ，y 方向的边界条 件也要进行相应的f o u r i e r 展开。 对于准二维问题的控制方程，采用在时间上的“分裂”格式、空间上的紧致 差分格式进行数值离散。 22 1 时间的离散方法 对于时间离散，我们采用k a r n i a d a k i s 、g e 、i s r a e l i m 和 o r s z a g s a ( 1 9 9 1 ) 建议的混合显一隐相结合的分裂格式，控制方程可写为以下 三个步骤： 学- - v m p ? 1 孥= 面1 v 2 z a ) 其中霸，酩为中间变量，a t ，尾，y o 是适当系数，按要求的精度取不同的值。 时间混合显一隐格式的系数 系数 y oa oo f ia 2舀。卢i卢2 1 阶精度 l lo01 0o 2 阶精度 3 2 2 1 2 02一l 0 3 阶精度 1 1 6 3 - 3 21 3 3 31 柑 置 v 腓 何 l 屏 ： 一 第二章数值计算方法与边界条件处理 为了满足连续性方程v - 。= 0 ，可以推出压力应满足的方程 v 。( v m p = 警 ( 2 2 5 ) ( 注：这里对压力梯度做散度的差分格式应与速度散度的差分格式相容) 。 2 2 2 空间的离散方法及边界条件处理 根据时间分裂格式，将问题分为三个部分分别进行处理：非线性对流项、压 力项和速度粘性项。 2 2 21 非线性的对流项 在计算方程的非线性项 0 7 ”v ) 厅” 时，先在谱空间计算出速度厅”的 导数，再将谱空间的速度及其导数变换到物理空间，在物理空间计算非线性项 ( 厅”。v ) 厅”，然后将非线性项变换为谱空间，得到l ( 面”v ) 舀” 。 可以用迎风紧致有限差分格式计算非线性项的一阶导数。这里给出付德薰、 马延文( 1 9 9 6 ) 建议的五阶迎风紧致格式来计算对流项。 a t 。 以v 兰为例，其五阶迎风紧致差分格式为： 砂 其中 型f + + 型， 22 ( 2 2 6 ) 2 ，+ h + 3 f + ，= ( 一u 。+ 2 + 1 2 u j _ 1 + 3 6 u ，一4 4 u ，l 一3 u f 一2 ) ( 1 2 a y ) ( 2 2 7 ) 3 f 一，+ 2 f 一= 一( - u 2 + 1 2 u + 3 6 u ，一4 4 u 一3 u 。+ 2 ) ( 1 2 a y ) ( 2 2 8 ) i 为y 方向的网格节点标号。 考虑空间周期流场，在流向采用周期性条件。 在边界点上( 江o ，n ) ，可以采用二阶精度的单边插值计算，其格式为 f + o = ( 一“2 + 4 “l 一3 “o ) ( 2 妙) f 一。= 一( - u 。2 + 4 “。一1 3 u 。) ( 2 a y ) 第。章数债计算方法与边界条仆处埋 在邻近边界点f i = 1 ，n 1 ) ，采用二阶精度紧致格式 f + o + 2 ，+ l = ( “2 + 4 “i 5 “o ) ( 2 a y ) 巧+ 2 f ：= ( 5 u 。一4 u 。一l “。一：) ( 2 a y ) 考虑剁流阳j 古j 期性条件，则以“宴为例，其五阶迎风紧致差分格式为 咖 具i j 业f十! 二幽f 2 f + ，i 十3 f + ，= ( 一甜，+ 二+ 1 2 u + 3 6 u j 一4 4 u f _ 1 3 u 卜2 ) ( 1 2 石) 3 f + ，+ 2 f = 一( “i2 + 1 2 u h + 3 6 u 。一4 4 u 一3 u ，+ 2 ) 0 2 a x ) ( 2 2 9 ) ( 2 2 1 0 ) ( 2 2 1 1 ) 计算j _ 域为：卜0 ， 一1 ，在左右边界处考虑( 2 2 1 0 ) 、( 2 2 1 1 ) ， 当i = 0 时( 几边界) ： 2 厂一，l + 3 f + n = ( - u 2 + 1 2 u 1 + 3 6 u ( ) 4 4 u 吖i 一3 u 吖一2 ) ( 1 2 a r ) 3 f 。，十2 fl = 一( - u 埘一? + 1 2 u 埘l + 3 6 u 。一4 4 u l 一3 u ! ) ( 1 2 k x ) 1 1 i = 1 时( 邻近片边界) 2 f n 十3 f + i = ( - u j 十1 2 u 2 + 3 6 u 】一4 4 u o 一3 u m1 ) ( 1 2 a x ) 3 fl + 2 f2 = 一( - u 一i + 1 2 u ) + 3 6 u 1 4 4 u 2 3 u 3 ) ( 1 2 a x ) m 一1 时( 边界) 2 f + ，、+ 3 f 1 m 】= ( 一“i + 12 u o + 3 6 u m 一1 4 4 u m 一? 一3 u m3 ) 0 2 a x ) 3 f ，l + 2 f 1 ) = 一( - u m 一3 十1 2 u m ，+ 3 6 u m4 4 u o 一3 u 1 ) 1 ( 1 2 a r ) i = m 一2l i , ( 邻近右边界) 第二章数值计算方法与边界条件处理 2 ，+ ，一3 + 3 f + 盯一2 = ( 一“o + 1 2 u m i + 3 6 u 吖一2 4 4 u 一3 3 u m - 4 ) ( 1 2 山c ) 3 f 一盯2 + 2 f 一吖一l = 一( - u ，一4 + 1 2 u ，一3 + 3 6 u ，一2 4 4 u 肼一l 一3 u o ) ( 12 缸) 由此形成如下的系数矩阵 32 3 2 2 2 23 23 32 3 2 3 氏 曩+ r 这两个方程用直接消元法进行求解即可。以此类推，非线性项中的其它项可以同 样处理。 对于亚格子应力模型项，使用与非线性项相同的伪谱方法进行处理。 2 2 22 压力h e i m h o f t z 方程的计算 压力满足的方程为： v 轨= 孚 令b = m 2 f 1 2 以吨，= 等，这样压力函数妒满足的方程为h e l h o l t z 方 程： 33 ， 3 2 第二章数值计算方法与边界条件处理 馨+等叫：厂0 x 。加2 采用紧致格式计算速度散度，以考为例，令，= 瓦0 u ，在内点采用六阶对称 紧致有限差分格式： f + 1 + 3 f + 只一。= 【( “，+ 2 一u i _ 2 ) + 2 8 ( u h u h ) ( 1 2 a y ) ( 2 2 1 4 ) 在邻近边界点上( f - l ，n 一1 ) ，用四阶精度的对称紧致格式 f “+ 4 f + e l = 3 ( u 一u h ) 6 y 在边界上( f _ o ，n ) ，采用三阶精度的紧致格式： e + 5 f , = 7 ( u 2 一u o ) + ( u 】一“3 ) 1 ( 2 y ) + 5 f 一i = _ 7 ( u 。一2 一u 。) + ( u 。一j u n 一3 ) 3 ( 2 妙) 对于宴，需要让其满足流向周期性条件，所以采用六阶对称紧致有限差分格式： o x 当i = 0 时 当i = 1 时 f + l + 3 鼻+ f l = ( “m u h ) + 2 8 ( u 一u - 1 ) ( 1 2 a x ) ( 2 2 1 5 ) 曩+ 3 f o + 一1 = ( “2 一u 一2 ) + 2 8 ( u 1 一“w 1 ) ( 1 2 a x ) 五+ 3 鼻+ f o = ( “3 一u m 1 ) + 2 8 ( u 2 一u o ) ( 1 2 k x ) 当i = 2 ，m 一3 时 ( + l + 3 ( + e 一1 = ( “。2 一“，2 ) + 2 8 ( u + l u ，1 ) ( 1 2 a x ) 第二章数值计算方法与边界条件处理 当i = m 一2 时 f _ 一l + 3 f m 一2 + f m 一3 = ( “o u ，一4 ) + 2 8 ( u ，一l u ，一3 ) 1 ( 1 2 6 x ) 当i = m 一】时 + 3 f 一l + f m 一2 = 【( “l u m 一3 ) + 2 8 ( u o u m - 2 ) 1 ( 1 2 6 x ) 形成如下的系数矩阵 3l 13 l l3 1 1 31 13 r f e ： 昂一， 一： r 一 得到的方程可用直接消元法进行求解。 对压力方程采用了九点四阶精度紧致有限差分格式，如图2 2 所示，最终得 到的离散方程： 哑纽纽! 二塑! ! 纽! 缝! ! 垒= ! 垒! ! 纽! 二! 坠! + ( z ) ! ! 竺! 壁生! 生! 二! ! 二! ! 生! ! ! 壁! = ! ! ! 生! j 受! = i 生! ! 生! = ! ! 二! ! 生! ! ( 每) 2 b ( 8 0 f + h + l 上l + 办一2 + 疵+ 2 ) = b y , + + l + 一1 + 丘一2 + 工+ 2 第一章数值计算方法与边界条件处理 11 吗2 两+ 丽4 刮3 则式( 2 2 1 6 ) 可化简为： 其中 1 02 ， 1 02 ， 丽一丽山2 丽一丽山 g i ，j = 鼋 l 。j 七 i “jr 七j t i 七j 。_ 七 q 。 压力方程为椭圆型方程，此类方程解的情况几乎完全取决于边界条件。本文 设压力方程边界条件为： ( 1 ) 流向周期条件； ( 2 ) 根据k a r nj a d a k i s ，i s r a e l i 和o r s z a g 的结论，在上下壁面处，压力 满足：害i = 0 。 g n l v = on 实际计算结果表明，这样的边界条件是能够满足精度要求的。计算区域如图2 3 所示： 卜厂上上 珈 务、 ，)g = 互 + j 2 一( 西q + r +谚 口 + 卜 谚2 口 + ， +旃 口 + 破o d+ j 一谚 口 + 一 j 吐2 4 + +谚 口 + 广 一垂4 a 第二章数值计算方法与边界条件处理 边界点 内点 同样，为满足流向周期性条件，式( 2 2 1 7 ) 在边界处有： 当i = 0 时： a 4 吮h ，一1 + 如西，一i + a 2 晚，一l + q 丸吨j + a o # o ，+ 4 1 卉，+ 口2 蟊，+ 1 + 口3 破j + 【+ 口4 丸一1 j + l = g 0 当i = m 一1 时 这样方程系数矩阵就变为一个加边的块三对角矩阵，可用超松驰( s o r ) 进行求 解。 对于典型线性方程：丘v = b ，其中a 为系数矩阵，x 为解向量，b 为方程右 端项。方程解的情况取决于系数矩阵a 的性态。在这里使用二阶精度单边差分对 压力场上下壁面边界条件进行离散： f 望 ：二! ! q 垫二丝：o ， l 砂。 2 d y 塾