一种强耦合Spalart-Allmaras湍流模型的RANS方程的高效数值计算方法.pdf

航空学报 A c t aA e r o n a u t i c ae tA s t r o n a u t i c aS i n i c a S e p 2 52 0 1 3V 0 1 3 4N o 9 2 0 0 7 - 2 0 1 8 I S S N1 0 0 0 - 6 8 9 3O N1 1 1 9 2 9 V h t t p ：dh k x b b u a a e d u c n h k x b b u a a e d u c n 一种强耦合S p a l a r t A l l m a r a s 湍流模型的 R A N S 方程的高效数值计算方法 杨小权1 2 ，杨爱明2 一，孙刚2 1 上海飞机设计研究院，上海2 0 1 2 1 0 2 复旦大学力学与工程科学系，上海 2 0 0 4 3 3 摘要：在工程实际中，一方程湍流模型或两方程湍流模型的求解通常和雷诺平均N a v i e r S t o c k e s ( R A N S ) 方程的求解 是解耦的，也称之为松耦合求解。在松耦合求解过程中，R A N S 方程和湍流模型方程通常采用不同的数值方法异步求 解。这种求解方式很容易产生因两者计算精度不一致而引起的额外数值耗散。为了消除这种耗散，将R A N S 方程与 S p a l a r t A l l m a r a s 模型方程耦合成一个系统方程强耦合R A N S 方程，并发展了一种用于求解该系统方程的高效强耦 合算法，其中对流项离散采用了R o e 格式，时间项的离散采用了隐式L U S G S ( L o w e r - U p p e rS y m m e t r i cG a u s s S e i d e l ) 格 式，为了提高计算效率，采用了三层V 循环多重网格方法。通过翼型机翼和振荡翼型机翼等算例验证了本文发展的 强耦合算法不仅具有较好的收敛性，而且计算精度明显优于松耦合算法，特别对于阻力的预测，强耦合算法更加准确。 关键词：S p a l a r t A l l m a r a s 湍流模型；R A N S 方程；强耦合算法；多重网格方法；双时间方法 中图分类号：V 2 1 1 3文献标识码：A文章编号：1 0 0 0 6 8 9 3 ( 2 0 1 3 ) 0 9 2 0 0 7 1 2 湍流问题至今仍然是困扰整个流体力学界的 一个难题，湍流的基本机理至今还没有完全弄清， 这就决定了各种湍流研究方法必然有各自的局限 性。湍流的数值模拟大致可分为三类：直接数值 模拟( D i r e c tN u m e r i c a lS i m u l a t i o n ，D N S ) 、大涡 模拟( L a r g eE d d yS i m u l a t i o n ，L E S ) 和求解雷诺 平均N a v i e r S t o c k e s ( R e y n o l d s A v e r a g e dN a v i e r S t o k e s ，R A N S ) 方程。鉴于前两种方法需要 极大的计算资源以及方法本身的问题，目前还无 法适用于高雷诺数工程问题的模拟，因而求解 R A N S 方程自然成为工程实际中的最佳方法。 对于一方程和两方程湍流模型而言，其求解 方式可分为两类：第一类是“松耦合”求解方式，即 在同一次时间推进中R A N S 方程和湍流模型方 程的求解相对独立；第二类是“强耦合”求解方式， 即在同一次时间推进中R A N S 方程和湍流模型 方程耦合求解。目前第一类求解方式在工程实际 中得到了广泛的应用；而第二类求解方式还比较 少见，但国内外学者已进行了一些有价值的研究。 K u n z 和L a k s h m i n a r a y a n a 1 研究了显式耦合两 方程肛湍流模型的R A N S 方程的稳定性，发现 显式强耦合算法和显式松耦合算法的收敛性大体 一致。L i u 和Z h a n g 2 1 发展了基于多重网格的显 式强耦合两方程肛叫湍流模型的R A N S 方程的 定常可压缩算法，计算结果发现强耦合求解方法 的收敛性优于松耦合解法。B a r a k o s 和D r i k a k i s 3 3 对比了耦合两方程k - e 和是一c U 湍流模型的 R A N S 方程的3 种求解方法( 即：隐式强耦合解 收稿日期：2 0 1 2 - 1 0 - 3 0 ；退修日期：2 0 1 2 - 1 1 - 2 2 ；录用日期：2 0 1 3 - 0 1 0 4 ；网络出版时间：2 0 1 3 0 5 2 31 1 ：2 4 网络出版地址：W W W c n k i n e t k c m s d e t a i l 1 1 1 9 2 9 V 2 0 1 3 0 5 2 3 1 1 2 4 0 0 2h t m I 基金项目：国家自然科学基金( 1 1 1 7 2 0 7 0 ) ；复旦大学研究生创新基金( E Y H 2 1 2 6 0 2 2 ) * 通讯作者T e l ：0 2 1 - 5 5 6 6 5 5 3 4E - m a i l ：a m y a n g I f u d a n e d u c n 引搠罄武tY a n gXo ，Y a n gAM ，S u nG A ne f f i c i e n tn u m e r i c a lm e t h o df o rc o u p l i n gt h eR A N Se q u a t i o n sw i t hS p a l a r t - A I I m a r a st u r b u l e n c e m o d e le q u a t i o n A c t a A e r o n a u t i c ae t A s t r o n a u t i c a S i n i c a 2 0 1 3 3 4 ( 9 ) ：2 0 0 7 - 2 0 1 8 杨小权。杨爱明t 孙嘲一种强耦合S p a l a r t - A I I m a r a s 湍流模型的R A N S 方程的高效数值计算方法航空学报2 0 1 3 ，3 4 ( 9 ) ：2 0 0 7 - 2 0 1 8 万方数据 航空学报S e p 2 52 0 1 3V 0 1 3 4N o 9 法、显式和隐式松耦合解法) ，研究发现隐式强耦 合求解方法最好。L e e 和C h o i 4 研究了耦合两方 程q - w 湍流模型的R A N S 方程的A F A D I ( A p p r o x i m a t e d F a c t o r i z a t i o n A l t e r n a t i v eD i r e c t i o n I m p l i c i t ) 解法，用于求解可压缩翼型和机翼的定 常黏性绕流问题，研究中还讨论了控制方程中的 湍动能项对收敛性的影响。Y a n g 5 1 和H u a n g 6 等研究了强耦合S p a l a r t A l l m a r a s 湍流模型的 R A N S 方程的五阶W E N O ( W e i g h t e dE s s e n t i a l l yN o n O s c i l l a t o r y ) 格式和预处理技术，用于求解 定常气动问题。J a h a n g i r i a n 和H a d i d o o l a b i 7 1 发 展了基于非结构网格的显式强耦合两方程肛e 湍 流模型R A N S 方程的可压缩算法，用于求解二维 俯仰振荡翼型的黏性绕流。这些研究工作大多集 中在定常问题的求解中，有关非定常问题的研究 还比较少见。湍流模型与R A N S 方程的耦合方 程，针对流动的输运特性，能够更精确地描述湍流 运动，同时采用统一的数值方法同步求解，减少了 数值耗散，提高了计算精度。相比于其他复杂的 湍流模型，S p a l a r t A l l m a r a s 8 1 模型的数值鲁棒性 较好，对于分离流动的精度较高，在航空领域得到 了广泛应用。S p a l a r t A l l m a r a s 湍流模型的控制 方程相对简单，便于实现和R A N S 方程的耦合求 解，且耦合求解的数值鲁棒性更好。 为了消除S p a l a r t A l l m a r a s 模型和R A N S 方 程在松耦合求解中因采用不同的数值方法异步求 解时引起的额外数值耗散，本文发展了基于S p a l a r t A l l m a r a s 模型方程和R A N S 方程相耦合的高 效强耦合求解方法，通过定常非定常算例分析了 强耦合求解方式的计算精度和收敛性能，对比了 两种不同耦合求解方式对翼型和机翼黏性绕流流 场气动力特性的影响。 1控制方程 在笛卡儿坐标系下将积分形式的可压缩 R A N S 方程与可压缩S p a l a r t A l l m a r a s 模型方程 进行耦合，可以得到如下的系统方程： f 譬d V + 书，。n d s 一书厂。n d s j v | I l d V o ( 1 ) 式中：V 为控制体体积；n 为控制体法向矢量；d s 为面积微元；守恒变量为 w 一 f Dp u l p u z p u 3 l D Ey 1 ( 2 ) 式中：p 和E 分别为密度和总能量；v 为涡黏性系 数相关量；U ，、U 。和“。分别为流体沿坐标轴z 、Y 和z 的速度分量。 在运动过程中网格是刚性的，没有发生变形。 考虑网格的运动速度批。，对流项通量和黏性通量 的分量厂。：和，。分别为 fc ：一 l 。t = p ( “：一U b i ) p ( u l U b i ) U l + P & I D ( “：一U b i ) “2 + P 3 i 2 p ( “。一U b i ) u 3 + P & 3 p ( U i U b i ) H + p u b i ( “。一) ； 。 0 盯l l C r i 2 盯，R ( 3 ) 式中：“i 为速度分量；H 为总焓；P 为压力；盯为热 传导系数；T 为温度；z i 为笛卡儿坐标系下的坐 标分量。源项矗一 o 0000 h 。 T 。 式( 3 ) 中的黏性应力鳓和源项分量h 。分别 定义为 = 肛( 芸+ 差) + 硪警 c 4 ， h 。一一垫巧V2 ；+ C b 。( 1 一f 。：) n ；+ C b l ( 1 一f t 2 ) f v z + 去一c 。) ( 号) 2 ( 5 ) 式中：肛为黏性系数；v 为分子黏性系数；d 为到壁 面的最近距离；8 i 为克罗内克占符号；n 为涡量。 为了封闭控制方程式( 1 ) ，引入状态方程： P = ( y 一1 ) p ( E 一虿1 州：) ( 6 ) 式中：y 为比热容比。方程式( 1 ) 采用有限体积法 在任意网格单元( i ，J ，最) 上进行离散，则可得空间 半离散控制方程 v 辈+ R ( ，) 一0 ( 7 ) 笪弧 塑兆“ K C + + 嘶 卜 1 一仃 万方数据 杨小权等：一种强耦合S p a l a r t A l l m a r a s 湍流模型的R A N S 方程的高效数值计算方法 嘉是；纂篇姜蓑萎吉是麓兹赇2 数值计算方法解，半离散控制方程式( 7 ) 引入伪时间导数项 一 V 娑+ V 挲+ R ( w ) 一0 ( 8 ) 2 1 双时间方法 d rO t 。 式( 8 ) 在伪时间方向的定常解即为物理时I 司上的 数值解。方程( 8 ) 中的物理时间项的求解采用隐 式L U S G S ( L o w e r U p p e rS y m m e t r i cG a u s s S e i d e l ) 方法。 可压缩流动的S p a l a r t A l l m a r a s 湍流模型方 程嘲的表达式为 8 _ E + 掣：S ，+ S 。+ D ( 9 ) 3 t a z ， 1。“。 式中：S P 、S D 和D 分别为S p a l a r t A l l m a r a s 湍流 模型生成项、破坏项和耗散项，其表达式分别为 S ，一C b 。( 1 一f 。) n ； ( 1 0 ) s 。一 c “ ( 1 一，t z ) f v z + 厂t z 方一 C w x ( 詈) 。 ( 1 1 ) D 一了C b 2v - 褰+ 盯o Z ，盯 a x i E V + ( 1 + c b 2 ) ； 亳) ( 1 2 ) 湍流黏性系数口。与变量；的关系式为 严t I D V 厂，1 ( 1 3 ) 模型中其余相关参数定义为 一杰，x 三詈 一，一去 一g ( 耀) “6 g r + c 。z ( r 6 一r ) ，r m i n ( 虿砉孑，1 0 ) 一川+ 赤，v z 2 一C 。3e x p ( - - C 。4 2 2 ) ( 1 4 ) 模珏! ! 中相关常数如T ： C b l 一0 1 3 55 ，C b 2 0 6 2 2 口一妻，k 一0 4 1 ，C ，1 7 1 ( 1 5 ) C 。，一百C b l + 盟堕，C 。2 一o 3 惫巧 C 。3 2 ，C 。3 1 1 ，C 。4 2 对控制方程式( 8 ) 的物理时间项进行二阶向 后差分离散，即 V 掣+ V 鲨之娑x - 型o r 厶 + R ( w 井1 ) 一0 ( 1 6 ) 将式( 1 6 ) 中的伪时间项进行一阶向后差分离散，即 V 等+ V 垫芝兰二A 型t +r Z R ( w 计1 ) 一0( 1 7 ) 对残值进行线化处理，即 ( 恙+ 老) H ( 善) “ A w 一 一! ! ：二垒丝：! 二：一R ( w m ) ( 1 8 ) 2 A t 、。 式中：r 为伪时间；t 为物理时间。式( 1 8 ) 可以采 用如下L U S G S 格式 9 3 进行求解 ( L + D ) D - 1 ( D + 【，) A w 一 一翌二竺笪一R ( ，m ) ( 1 9 ) 2 A t 。、 L U 分解矩阵L 、D 和【，的表达式分别为 f L 一一( A 三。女+ B 去，_ l ， + 融卅一) D _ 訾+ 等州甄+ 卯+ O c ， J l U 一( A ；1 _ + B i J - 1 ， + C 0 ，抖1 ) ( 2 0 ) 式中：O A 、叮B 和盯c 分别为J a c o b i a n 矩阵A 、B 和C 的谱半径；口为一个大于1 的常数，在本文计算中 取1 2 。口取值越大，格式稳定性越好，但收敛速 度会降低。 式( 1 9 ) 求解分三步进行，即：向前扫描、向后 扫描和更新变量W 。 在式( 2 0 ) 中的J a c o b i a n 矩阵只考虑了对流 通量，但是黏性通量和源项的存在，对控制方程的 收敛性有一定的影响，特别是源项，它的存在对控 制方程的收敛性有较大影响。因此，在线化处理 过程中要对J a c o b i a n 矩阵的特征值进行修正，现 以e 方向为例，给出其具体表达式： 盯：5G A + 劳：。+ S ； ( 2 1 ) 0 r T 式中：P r 为普朗特数；S ；的表达式为 万方数据 航 空 学报S e p 2 52 0 1 3V o l3 4N o 9 婪+ 拿1 v f 聋+ 一a S p1 o a va v a ya v ( 等+ 等) 。 ( 2 2 ) 对于定常计算只需令物理时间相关项为0 即可。 2 2空间离散 空间离散采用格心有限体积法。单元界面处 的对流通量采用基于R o e 的近似R i e m a n n 解的 G o d u n o v 通量函数近似 1 ，即 1 f 。，一( f 。，( W x ) + f 。，( w - 。) ) 一 厶 1 lAI ( W R W L ) ( 2 3 ) 厶 式中：w 。，、W 。分别为控制体界面左、右两侧的守恒 变量；A 为R o e 平均J a c o b i a n 矩阵。为抑制非物 理振荡，对J a c o b i a n 矩阵的特征值进行熵修正， 这里采用H a r t e n 的熵修正 1 1 。为了提高格式精 度，界面两侧守恒变量w 。和w 。的值通过 M U S C L 插值得到，为了避免插值过程中引入非 物理的振荡，采用了v a nL e e r 通量限制器口2 | 。 2 3 初始条件和边界条件 1 ) 初始条件 流场的守恒变量的初始值由自由来流给定， 与湍流模型相关的变量v 和湍流黏性系数口。的 初值为 护“一( 2 4 ) l z t 。t l D o 。V o 。f 。1 2 ) 边界条件 速度在物面满足无滑移边界条件，变量v 和 湍流黏性系数肚。在物面为0 ，压力满足法向压力 梯度为0 。流场守恒变量远场边界采用R i e m a n n 不变量处理，变量v 和口。在远场通过出入流条件 来确定，即 卜- 1 3 4 1 9 4 6 ( 2 5 ) I 岸。一0 0 0 9 2 4 多重网格方法 多重网格方法 1 3 - 1 4 3 就是利用一系列粗细不 同的网格来求解同一离散问题，将迭代误差分为 高频波动误差和低频光滑误差，在细网格上利用 有效的迭代方法消除高频波动误差；利用粗网格 消除低频光滑误差。通过插值算子在粗网格与细 网格间传值。其主要的两个过程就是限制和插值 的过程，限制过程就是细网格向粗网格传值的过 程；插值过程就是粗网格向细网格传值的过程。 在限制过程中本文采用体积加权平均法，插值过 程采用直接赋值方法。 3算例分析 为了验证本文发展的S p a l a r t A l l m a r a s 湍流 模型和R A N S 方程相耦合的强耦合算法在求解 可压缩黏性绕流中的收敛性和计算精度，选择了 以下4 个算例，计算条件如表1 和表2 所示。表 中：M a 为马赫数；O l 为迎角；R e 为雷诺数；口。为平 均迎角口。为俯仰运动迎角的幅值。 表1R A E2 8 2 2 翼型和O N E R AM 6 机翼计算条件 T a b l e1F l o wc o n d i t i o n sf o rRA E2 8 2 2a i r “Ia n dO N E R A M 6w i n g 表2A G A R D 振荡翼型实验数据 T a b l e2A G A R D e x p e r i m e n t a ld a t ao fo s c i l l a t i n ga i r f o i l 在采用双时间方法进行计算时，每个时间周 期分为3 6 0 个时间等份，每个物理时间步的最大 迭代次数为3 0 ，共进行4 个周期的计算。 所有计算是在D E L L 台式电脑上进行，该机 器性能参数为：I n t e l ( R ) C o r e ( T M ) i 7 2 6 0 0k C P U3 4G H z ，1 6G B 内存。 翼型和机翼的俯仰运动规律为 a ( ) 一口。+ O t os i n ( t ) ( 2 6 ) 式中：为俯仰振荡角速度。是由无量纲的减 缩频率k 定义，k 与的关系为 k 一f ( 2 U 一) ( 2 7 ) ， 0 ，J，、l 一 万方数据 式中：C 为翼型的根部弦长；U 。为自由来流速度。 3 1R A E2 8 2 2 翼型跨声速定常黏性绕流 算例采用C 型网格，网格大小为2 0 1 4 9 ，其 中翼型上分布1 6 0 个网格单元，第一层网格离壁 面的距离为1 0 1 0 一，远场距离区在大约1 5 倍 弦长处。 图1 给出了强耦合算法和松耦合算法计算的 R A E2 8 2 2 翼型压力系数C 。分布与实验值的对 比。从图1 可知，强耦合算法计算的压力系数要 比用松耦合算法计算的压力系数更加接近实验 值，且捕捉到的激波位置更为准确。 图1 翼型表面压力系数计算值与实验值的对比 F i g 1C o m p a r i s o no fc o m p u t e da n de x p e r i m e n tp r e s s u r e c o e f f i c i e n t sa ta i r f o i ls u r f a c a 表3 给出了强耦合算法和松耦合算法计算的 力系数与实验值的比较。从表3 可知，强耦合算法 计算的升力系数C L 和力矩系数G 与实验值更加 接近，松耦合算法计算的阻力系数C D 与实验值更 加接近些，但两种方法计算的阻力系数相差不大。 表3R A E2 8 2 2 翼型的力系数计算值与实验值的比较 T a b l e3 C o m p a r i s o no fc o m p u t e da n de x p e r i m e n tf o r c eC O - e f f i c i e n t so fR A E2 8 2 2a i r f o i l 图2 给出了强耦合算法和松耦合算法计算的 R A E2 8 2 2 翼型的残值l g ( R E S R E S o ) 和升、阻力 系数( C 。和C D ) 收敛曲线。 l t e r a t i o n 1 0 3 ( a ) C o n v e r g e n c eh i s t o r yo f r e s i d u a l s 1 t e r a t i o n 1 0 ( b ) C o n v e r g e n c eh i s t o r yo f l i f tc o e f f i c i e n t s I t e r a t i o n 1 0 3 ( c ) C o n v e r g e n c eo f d r a gc o e f f i c i e n t s 图2 残值和升、阻力系数收敛曲线 F i g 2C o n v e r g e n c eh i s t o r yo fr e s i d u a l s ，l i f ta n dd r a gc o e f f i c i e n t s 万方数据 航空学报 S e p 2 52 0 1 3V 0 1 3 4N o 9 从残值收敛图可知，强耦合算法要达到收敛 标准，单层网格需要约70 9 0 次迭代；三层多重网 格需要约7 8 0 次迭代；松耦合算法要达到收敛标 准，单层网格需要69 1 0 次迭代，三层多重网格需 要7 4 0 次迭代。就强耦合算法的计算时间而言， 单层网格的每步计算时间约为0 0 3 1S ，三层多重 网格的每步计算时间约为0 0 4 3S ，在达到收敛时 单层网格需要近2 2 0S ，三层多重网格只需要 3 4S ，加速比约为6 5 ；就松耦合算法的计算时间 而言，单层网格每步需要0 0 2 9S ，三层多重网格 每步需要0 0 4 0S ，在达到收敛时单层网格需要近 2 0 1S ，三层多重网格只需要3 0S ，加速比约为 6 7 。从升、阻力系数收敛图可知，多重网格方法 对升、阻力系数的加速收敛效果非常明显，强耦合 算法的升、阻力系数收敛速度要略快于松耦合 算法。 3 2O N E R AM 6 机翼跨声速定常黏性绕流 算例采用C _ H 型网格，网格大小为2 0 1 4 9 6 1 ，其中翼型上分布1 6 0 4 0 个网格单元，第一层 网格离壁面的距离为1 0 1 0 一，远场距离区在 大约1 5 倍弦长处。 图3 给出了强耦合算法和松耦合算法的 O N E R AM 6 机翼表面压力系数分布计算结果与实 验值的对比。由图3 可知，强耦合算法比松耦合算 法计算的激波位置更接近实验值，特别在机翼展向 z 6 6 5 和z 6 8 0 处，强耦合算法能够更加 准确地捕捉到机翼上表面的A 激波。 ( d ) z b - = ：$ 0 图3O N E R AM 6 机翼表面压力系数的计算值与实验值 的对比 F i g 3C o m p a r i s o no fc o m p u t e da n de x p e r i m e n tp r e s s u r e c o e f f i c i e n t sf o rO N E R AM 6w i n gs u r f a c e 图4 给出了强耦合算法和松耦合算法计算的 O N E R AM 6 机翼残值和升力系数的收敛曲线。 从残值收敛图可知，强耦合算法在达到收敛标准 时，单层网格需要约75 2 0 次迭代；三层多重网格 万方数据 杨小权等：一种强耦合S p a l a r t A l l m a r a s 湍流模型的R A N S 方程的高效数值计算方法 需要约9 1 0 次迭代；松耦合算法在达到收敛标准 时，单层网格需要约76 4 0 次迭代，三层多重网格 需要约8 3 0 次迭代。就强耦合算法而言，单层网 格每步平均时间为1 2 3S ，三层多重网格每步平 均时间为1 6 4S ，在达到收敛时单层网格需要约 2 5 7 0h ，多重网格只需要4 1 5h ，加速比约为 6 2 ；就松耦合算法而言，单层网格每步平均时间 为1 1 6S ，三层多重网格每步平均时间为1 5 5S ， 在达到收敛时单层网格需要约2 4 6 2h ，多重网 格只需要3 5 8h ，加速比约为6 9 。从升力系数 的收敛图可知，无论是否采用多重网格方法，强耦 合算法的升力系数收敛速度要略快于松耦合算法。 I t e r a t i o n 1 0 3 ( a ) C o n v e r g e n c eh i s t o r yo fr e s i d u a l s I t e r a t i o n 1 0 3 ( b ) C o n v e r g g n c eh i s t o r yo f l i f tc o e f f i c i e n t s 图4 残值和升力系数收敛曲线 F i g 4C o n v e r g e n c eh i s t o r yo fr e s i d u a l sa n dl i f tc o e f f i c i e n t s 组，计算量上略有增加。强耦合算法和松耦合算 法的单步计算时间差异不大，但强耦合算法的总 体计算时间略有增加。 3 3 二维翼型俯仰振荡黏性绕流 算例采用的计算网格和算例与3 1 节一致。 图5 给出了N A C A0 0 1 2 翼型俯仰振荡时，强耦 合算法和松耦合算法计算的升、阻力系数收敛历 程。由图5 可知，无论是升力系数还是阻力系数 均需要迭代4 个周期才能得到完整的周期收敛 解。强耦合算法计算的升力系数及阻力系数均低 于松耦合算法，阻力系数的这种差别尤为明显。 I t e r a t i o n 1 0 4 a ) C o n v e r g e n c eh i s t o r y o f l i f tc o e f f i c i e n t s I t e r a t i o n 1 0 4 ( b ) C o n v e r g e n c el l i s t o r yo f 血a gc o e f f i c i e n t s 图5N A C A0 0 1 2 翼型升力系数和阻力系数收敛曲线 F i g 5C o n v e r g e n c eo fl i f ta n dd r a gc o e f f i c i e n t so fa no s c i l l a t i n gN A C A0 0 1 2a i r f o i l 图6 给出了强耦合算法和松耦合算法的升力 就计算量而言，松耦合方式求解的是一个五 系数和力矩系数迟滞环计算值与实验值的对比。 元方程组，而强耦合方式求解的是一个六元方程 由图6 可知，两种方法计算得到的升力系数差别 万方数据 航空学报S e p 2 52 0 1 3V 0 1 3 4N o 9 不大，均与实验值吻合得很好，而计算的力矩系数 的差别特别明显，强耦合算法计算的力矩系数与 实验值吻合得更好，而松耦合算法计算的力矩系 数与实验值差别比较大。 川。) ( b ) M o m e n tc o e f f i c i e n tl o o p s 图6N A C A0 0 1 2 翼型升力环和力矩系数环的计算值与 实验值的对比 F i g 6 C o m p a r i s o no fc o m p u t e da n de x p e r i m e n tf o rl i f t a n dm o m e n tc o e f f i c i e n tl o o p so fa no s c i l l a t i n gN A C A0 0 1 2a i r f o i l 3 4 三维机翼俯仰振荡黏性绕流 算例采用的计算网格和算例与3 2 节一致。 图7 给出了强耦合算法和松耦合算法计算的俯仰 振荡L a n n 机翼升力系数和阻力系数的收敛历 程。由图7 可知，计算需要进行4 个周期才能获 得完整的周期解。就力系数而言，强耦合算法计 算的升、阻力系数均小于松耦合算法的计算值，这 和二维振荡翼型计算的结果一致。 I t e r a t i o n 10 4 ( b ) C o n v e r g e n c eo fd r a gc o e f f i c i e n t s 图7L a n n 机翼的升力系数和阻力系数收敛曲线 F i g 7C o n v e r g e n c eo fl i f ta n dd r a gc o e f f i c i e n t so fa no s c i l l a t i n gL a n nw i n g 为了更加准确地分析计算所得压力系数分布 的周期性，将压力系数进行傅里叶展开 1 弘20 | ，且 只保留一阶分量，即 C p C p M 。+ C p R 。s i n ( c o t ) + C p I 。c o s ( w t ) ( 2 8 ) 式中：C ，M e a n 、C ，。和C ，。分别为压力系数的平均 值、同相分量和正交分量。 图8 、图9 和图1 0 分别给出了强耦合算法和 松耦合算法计算的L a n n 机翼在不同展向位置的 压力系数分布平均值、同相分量和正交分量与实 验值的对比。由图8 图1 0 可知，两种算法计算 得到的不同展向位置的压力系数分布平均值与同 相分量差别不大，均与实验值吻合得较好，主要的 差别就在于压力系数分布的正交分量，强耦合算 法的计算值更接近实验值。 万方数据 杨小权等：一种强耦合S p a l a r t A l l m a r a s 湍流模型的R A N S 方程的高效数值计算方法 2 0 1 5 壶 U 4 x l c ( a 1 ：b - 2 0 0 0 20 40 60 81 0 x l c ( d ) z l b = 6 5 x c ( b l d b = 3 25 OO 2O 4O 6O810 00 20 40 60 8100 x l c ( e ) z b = 8 2 5 图8L a n n 机翼俯仰振荡压力系数平均值与实验值的对比 O 20 40 608l0 x l c ( f ) z l b - - 9 5 F i g 8C o m p a r i s o no fc o m p u t e da n de x p e r i m e n tm e a np r e s s u r ec o e f f i c i e n t so fa no s c i l l a t i n gL a n nw i n g 工C ( a ) 爿b - 2 0 x l c ( d ) z b = 6 5 岂 U 4 工C ( b 1 b = 3 25 鲋c ( e ) z b = 8 2 5 X ，c f c l ：b = 4 75 图9L a n n 机翼俯仰振荡压力系数同相值与实验值的对比 x l c ( f ) z l b - - 9 5 F i g 9C o m p a r i s o no fc o m p u t e da n de x p e r i m e n tr e a lp r e s s u r ec o e f f i c i e n t so fa no s c i l l a t i n gL a n nw i n g 万方数据 2 0 1 6 航 空 学报S e p2 52 0 1 3V o l3 4N o9 X f ( a ) z b = 2 0 X i c ( d ) z b = 6 5 O 2 01 一O1 O 2 乏迭茹抟磷一睁 。S t r o n g l yc o u p l e d Looselyc o u p l e d 。E x p e r i m e n t a ld a t a ( L o w e r ) o E x p e r i m e n t a ld a t a ( U p p e r ) O O 2 三 U0 - 0 毯一 C ( e ) z b = 8 25 O 2O 40 6080 x l c ( c ) z b = 4 75 X C ( Dz l b = 9 5 图1 0L a n n 机翼俯仰振荡压力系数正交量与实验值的对比 F i g 1 0C o m p a r i s o no fc o m p u t e da n de x p e r i m e n ti m a g i n a r yc o e f f i c i e n t so fa no s c i l l a t i n gL a n nw i n g 4 结论 为了消除S p a l a r t A l l m a r a s 模型和R A N S 方 程在松耦合求解中因采用不同的数值方法异步求 解时引起的额外数值耗散，发展了基于S p a l a r t A l l m a r a s 模型方程和R A N S 方程相耦合的强耦 合算法，用于求解翼型和机翼的可压缩黏性扰流， 通过算例对比分析了强耦合算法的收敛性和计算 精度。计算结果表明： 1 ) 对于跨声速定常流动问题而言，强耦合算 法具有很好的收敛性，三层多重网格加速收敛的 效果十分明显，加速比约为6 2 ；松耦合算法的加 速比约为6 7 。强耦合算法预测的升、阻力系数 略微偏大与实验值吻合得更好，强耦合算法的计 算精度明显优于松耦合算法。 2 ) 对于周期性非定常流动问题而言，强耦合 算法的计算精度明显优于松耦合算法，预测的升、 阻力系数略微偏小且更加准确，力矩系数迟滞环 与实验值吻合得更好。 1 2 3 3 4 5 6 参考文献 K u n zR F ，L a k s h m i n a r a y a n aB S t a b i l i t yo fe x p l i c i tN a v i e r - S t o k e sp r o c e d u r e su s i n gk - ea n d a l g e b r a i cR e y n o l d s s t r e s st u r b u l e n c em o d e l s J o u r n a lo fC o m p u t a t i o n a lP h y s i C S ，1 9 9 2 ，1 0 3 ( 1 ) ：1 4 1 1 5 9 L i uF ，Z h a n gX As t r o n g l yc o u p l e dt i m e - m a r c h i n gm e t h o df o rs o l v i n gt h eN a v i e r S t o k e sa n dk - c ot u r b u l e n c em o d e l e q u a t i o n sw i t hm u l t i g r i d J o u r n a lo fC o m p u t a t i o n a lP h y s i C S ，1 9 9 6 ，1 2 8 ( 2 ) ：2 8 9 3 0 0 B a r a k o sG ，D r i k a k i sD I m p l i c i tu n f a e t o r e di m p l e m e n t a t i o n o ft w o - e q u a t i o nt u r b u l e n c em o d e l si n c o m p r e s s i b l e N a v i e r - S t o k e sm e t h o d s I n t e r n a t i o n a lJ o u r n a lf o rN u m e r i c a lM e t h o d si nF l u i d s ，1 9 9 8 ，2 8 ( 1 ) ：7 3 9 4 L e eS ，C h o iDW O nc o u p l i n gt h eR e y n o l d s - a v e r a g e d N a v i e r - S t o k e se q u a t i o n sw i t h t w o - e q u a t i o n t u r b u l e n c e m o d e le q u a t i o n s I n t e r n a t i o n a lJ o u r n a lf o rN u m e r i c a l M e t h o d si nF l u i d s ，2 0 0 6 ，5 0 ( 2 ) ：1 6 5 1 9 7 Y a n gJ ，H s i e hT ，W a n gC I m p l i c i tw e i g h t e de s s e n t i a l l y n o n o s c i l l a t o r ys c h e m e sw i t ha n t i d i f f u s i v ef l u xc o m p r e s s i b l e v i s c o u sf l o w s A I A AJ o u r n a l ，2 0 0 9 ，4 7 ( 2 ) ：1 4 3 5 1 4 4 4 H u a n gJ ，L i nH ，Y a n gJ I m p l i c i tp r e c o n d i t i o n e dW E N O s c h e m ef o rs t e a d yv i s c o u sf l o wc o m p u t a t i o n J o u r n a lo f 万方数据 杨小权等：一种强耦合S p a l a r t A l l m a r a s 湍流模型的R A N S 方程的高效数值计算方法 2 0 1 7 E 7 3 8 9 1 0 3 1 1 1 2 1 3 1 4 15 C o m p u t a t i o n a lP h y s i c s ，2 0 0 9 ，2 2 8 ( 2 ) ：4 2 0 4 3 8 J a h a n g i r i a nA ，H a d i d o o l a b iM A ni m p l i c i ts o l u t i o no ft h e u n s t e a d yN a v i e r - S t o k e se q u a t i o n so nu n s t r u c t u r e dm o v i n g g r i d s I C A S - 2 4 ，2 0 0 4 S p a l a r tPR ，A l l m a r a sSR Ao n e e q u a t i o nt u r b u l e n c e m o d e lf o ra e r o d y n a m i cf l o w s L aR e c h e r c h eA e r o s p a t i a l e ， 1 9 9 4 ，1 ( 1 ) ：5 - 2 1 Y o o nS 。J a m e s o nA Am u h i g r i dL U S S O Rs c h e m ef o r a p p r o x i m a t eN e w t o ni t e r a t i o na p p l i e dt Ot h eE u l e re q u a t i o n s N A S A C R 一1 7 9 5 2 4 。1 9 8 6 R o ePL A p p r o x i m a t eR i e m a n ns o l v e r s ，p a r a m e t e rv e c t o r sa n dd i f f e r e n c es c h e m e s J o u r n a lo fC o m p u t a t i o n a l P h y s i c s ，1 9 8 1 ，4 3 ( 2 ) ：3 5 7 3 7 2 H a r t e nA ，H y m a nJM S e l fa d j u s t i n gg r i dm e t h o d sf o r o n e d i m e n s i o n a lh y p e r b o l i cc o n s e r v a t i o nl a w s J o u r n a lo f C