NS方程激波计算的摄动有限差分方法.pdf

N S 方程激波计算的摄动有限差分方法N S 方程激波计算的摄动有限差分方法申义庆高智杨国伟中国科学院力学研究高温气体动力学实验室，北京，1 0 0 0 8 0摘要：摄动有限差分( P F D ) 方法从～阶迎风差分格式出发，将差分系数展开为网格步长的幂级数，通过提高修正微分方程的逼近精度来获得更高精度的差分格式由于格式基于一阶迎风格式，因此具有迎风效应、 网格节点少等特点本文首先通过对B u r g e r s 方程的摄动差分格式的推导，将摄动有限差分格式引入时间相 关法的计算，并构造了守恒形式的摄动有限差分格式，然后推广到一维N a v i e r - S t o k e s 方程组的计算数值比较研究表明：本文构造的N s 方程摄动有限差分格式具有比～阶迎风较高的精度和分辨率，而且保持了一阶迎风格式的无振荡性质O 引言高精度、高分辨率差分格式是计算流体力学数值方法研究的目的之一一般来说，差分格式的精度越高，使用的网格点越多，而对间断解的高分辨率性质则通过采用附加条件的限制来实现，如二十年来提出和发展的T V D 、E N O 、N N D 等格式u 3 0 ，但由于所涉及的网格点较多，给边界的处理带来一定的困难。

紧致格式H 1 利用较少的节点能获得较高的精度，但在计算激波时会带来非物理振荡，而且由于格式中利用了节点上的导数值，对于复杂边界的情况也难以给定数值摄动方法基于一阶迎风差分或二阶中心差分格式，对差分方程的非微商项进行数值摄动展开，通过提高修正微分方程的逼近精度来获得高精度的格式，由于摄动格式只利用了三个节点，具有迎风性等特点，因此在计算流体力学中将具有广阔的应用前景数值摄动有限差分格式的高精度、高分辨能力已在众多的模型方程和不可压缩流动的计算中得到了应用和验证㈣但对可压缩流动的计算尚需进一步的发展本文首先通过对B u r g e r s 方程的摄动差分格式的推导，将摄动有限差分格式引入时间相关法的计算，并构造了守恒形式的摄动有限差分格式，然后推广到一维N a v i e r —S t o k e s 方程组的计算数值比较研究表明：本文构造的N S 方程摄动有限差分格式比一阶迎风具有较高的精度和分辨率，而且保持了一阶迎风格式的无振荡性质1 模型万程搬动有限爱分万法对B u r g e r s 方程要+ u 罢：去窘㈤1一十一= ———_L ，m孤R e 舐2其在定常状态时为“耋：土垒( 2 )“面2 面萨蟛’对方程( 2 ) 中一阶导数采用一阶迎风格式，二阶导数采用二阶中心格式，有如下差分方程 一l + o t “，竺[ 堕+ 坐z ，，燮：土坚二靼( 3 )2o缸2。

缸R e缸2其中，Q = s i g n ( u ，) 设对差商系数U ，进行如下摄动展开“ 彦“ - 1 ／J ( 1 + 么l △x + 么2 A x 2 + 彳3 △x 3 + 么4 △x 4 + ⋯)( 4 )5 1第十一届全国激波与激波管学术会议将( 4 ) 代入( 3 ) 后可得修正微分方程为 嘣詈卜善燕碧恤卜去r 窘”善赢器归2 Ⅲ5 ，在J 点对方程( 2 ) 进行局部点凝固系数，可得关系式≥沪( 眦广1 罢I ，在( 5 ) 中对缸”利用待定系数法，则可求得系数4 丽1 ( 口R e u j ) ”，刀= o ，1 ，2 ，⋯由此，将“归代替( 3 ) 中的差商系数”／，则可得n 阶得摄动差分格式：半“归警+ 半“印警= 丙1 笪等塑如对方程( 2 ) 利用时间相关法求解，即求方程( 1 ) 的定常解，则n 阶格式为 詈”半“，訾+ 半“，警= 南其中，“』+ l 一2 u 』+ 材／一对时间项詈I ／的求解可采用R u n g e —K u t t a 方法 6 )利用一阶导数的守恒型格式，即将( 9 ) 式中左端的空间导数的离散用迎风守恒型格式代替，并对右端作相应变化，可以将( 9 ) 推广到守恒型摄动有限差分方法，塑1 + 堕! [ 生竺：—L f 堑二竺一竺兰幽 讲I ，’A xR e A x ”A ，+ l ，2A ，一l ，2 ’( 1 0 )本文中，忍州，：采用R o e 型迎风格式，即乃Ⅳ，：= 圭( 乃+ 。

‘一I 口州，zI △j + v 2 u ) ，其中乃= m ，) = 丢即2蚧∥氆H 飞一肌= j 等A j + l ／2 ≠：0 ，4 ，+ l ，2 = 1 + 4 l A x + A 2 △x 2 + 4 3 缸3 + 4 4 缸4 + ⋯，么 丽1 @ R e ‰z ) ”，力= 0 ’l ，2 ，⋯·a ( u J )A 』“，22 02 一维N a vie r —S t o k e s 方程的摄动有限差分格式一维N a v i e r ．S t o k e s 方程无量纲化后可写为守恒形式：5 2缸4∑瑚I l4N S 方程激波计算的摄动有限差分方法或非守恒形式：( 1 l a )( 1 l b )设A = S - 1 A S ，其中A = 讲昭( ^ ，允：，丑) 和s = ( s l ，S 2 ，S 3 ) r 分别为A 的特征值矩阵和为左特征向量矩阵，关于方程中有关表达式可参见文献【9 】若对对流项采用一阶迎风格式，粘性项采用二阶中心格式，则有 詈”s ‘1 抄s 鲋渺半∥抄s 龇凇掣一E J + l ，2 一E 卜1 ，2 A x类似第2 节，利用数值摄动方法，构造了如下的摄动有限差分格式 ≯O U + J ‘i A 叶s 鲋渺半肌n 旺= G∥拉s 鲋邪掣一lE J + l ，2 一只』埘G m = 噻志( “ 纫∽矿】- 1 ，( 1 2 )( 1 3 )G f = 噻石‰( 气s 柳( ¨门一，氏= 萼旦堡笋，气= 三笋鱼笋。

同样有守恒型摄动有限差分格式如下： 詈”学= 去( G 川- 1 属j + ] 1 2 - - G ；! m ‰：)( 1 4 )本文计算中，_ 川：= 丢( 乃+ ，+ ‘一s 君IA j + I ／2 [ S j + I ／a A j + I ／，_ U ) ，％雌1 1 2 卜川z 《扣一川小州～，G 『j + l 陀= 【丢N 丽面1 ( 气哟玎( “j + l 彪) ) ”】q ，氏= T 3 R eP j + 1 1 2 ”』+ l ，2 缸5 3等等=y —x》譬沁～击沁一甜OO qOG0字半耻3 数值算例( 1 ) ．对方程( 2 ) ，有解析解为U = t a n h ( - R e x ／2 ) ，一L ≤x ≤L数值计算中，取边界条件为：“( 一L ，f ) = t a n h ( R e L ／2 ) ，u ( L ，f ) = t a n n ( 一R e L ／2 ) ，时间离散采用二阶R u n g e —K u t t a 格式，收敛准则：专薯l 甜尹一“／nI ／△f ≤．5 ×l ～，摄动有限差分格式、守恒型摄动有限差分格式与一阶迎风格式的计算结果分别见．表1 、表2 。

表1 ．非守恒型摄动有限差分格式( P F D ) ( n = 4 时) 与一阶迎风格式( 卜u w ) 计算结果的误差比较计算条件格式最大误差平均误差R e = 1 0 0 ，N = 8 01 一U W0 ．3 7 1 4 8 9 e 一10 ．1 2 3 5 4 1 e 一1L = O ．2 ．R w = O ．2 5P F D0 ．1 0 6 9 3 5 e 一20 ．3 2 6 2 6 8 e 一3R e = 1 0 0 ，N = 8 01 - U W0 ．2 0 3 0 4 3 E + 00 ．8 4 0 2 5 2 e 一2L = 2 ．R l f ：2 ．5P F D0 ．4 3 7 0 9 8 e 一10 ．1 1 8 4 4 6 e 一2R e = 1 0 0 0 ，N = 8 01 - U W0 ．2 0 2 9 9 6 e + 00 ．8 3 9 8 9 1 e 一2L = O ．2 ．R 一= 2 ．5P F D0 ．4 3 7 6 6 2 e 一10 ．11 8 6 5 l e 一2R e = 1 0 0 0 ，N = 8 0l —U W0 ．3 9 9 9 0 5 e 一10 ．1 0 2 7 9 2 e 一2L = 2 ．R 。

2 5P F D0 ．6 2 4 9 0 7 e - 40 ．2 1 8 3 9 1 e 一5R e = 1 0 0 0 0 0 ，N = 8 01 一U W0 ．4 0 0 4 9 2 e 一20 ．9 9 8 6 51e 一4L = O ．2 ．‰= 2 5 0P F DO ．5 1 1 6 8 2 e 一50 ．7 1 9 3 3 9 e 一6R e = 1 0 0 0 0 0 ，N = 8 01 - U W0 ．4 4 8 6 8 0 e 一30 ．i 1 6 5 1 5 e 一4L = 2 ．R F 2 5 0 0P F D0 ．5 2 9 5 2 6 e 一40 ．1 9 2 11 3 e 一5表2 ．守恒型摄动有限差分格式( P F D ) ( n = 4 时) 与守恒一阶迎风格式( 卜u w ) 计算结果的误差比较计算条件格式最大误差平均误差R e = 1 0 0 ，N = 8 01 一U W0 ．2 1 8 8 5 5 e 一10 ．7 4 3 9 0 3 e 一2L = O ．2 凡= O ．2 5P F D0 ．9 3 0 8 5 8 e 一30 ．2 4 6 1 6 7 e 一3R e = 1 0 0 ，N = 8 0l —U WO ．1 7 1 2 9 8 E + 00 ．7 1 8 5 5 2 e 一2L = 2 。

R N = 2 ．5P F D0 ．5 5 4 9 3 7 e 一10 ．1 5 8 0 2 2 e 一2R e = 1 0 0 0 ，N = 8 01 一U WO ．1 7 1 2 5 5 e + 00 ．7 1 8 1 8 7 e 一2L = O ．2 ，凡= 2 ．5P F D0 ．5 5 4 4 2 8 e 一10 ．1 5 7 8 1 2 e 一2R e = 1 0 0 0 ，N = 8 01 一U W0 ．3 9 2 5 1 3 e —l0 ．1 0 0 9 1 7 e 一2L = 2 ，R ．= 2 5P F DO ．17 6 4 7 9 e 一30 ．5 2 7 0 6 2 e 一5R e = 1 0 0 0 0 0 ，N = 8 0l - U W0 ．3 9 9 7 2 2 e 一20 ．9 9 7 3 9 3 e 一4L = O ．2 ．R M = 2 5 0P F D0 ．4 9 9 9 7 5 e 一50 ．7 4 0 1 7 9 e 一6R e = 1 0 0 0 0 0 ，N = 8 01 - U W0 ．4 5 0 8 3 6 e 一30 ．1 2 0 5 7 5 e 一4L = 2 。

P 、- - = 2 5 0 0P F D0 ．5 0 1 6 5 5 e - 40 ．2 1 4 3 0 3 e 一5从表1 中可看出，P F D 格式的最大误差要比一阶迎风格式的最大误差最大可小到2 个数量级，如最大网格雷诺数R 2 ．5 时当网格雷诺数R ～D ( 1 ) 时，P F D 格式最大误差相比于 一阶迎风格式提高不是很大，但也能提高3 倍以上，如R e = 1 0 0 0 ，N = 8 0 ，L = O ．2 ，R 庐2 ．5 的情形从表中可看出，当网格雷诺数介于0 ( 1 0 ) ～O ( 1 0 0 ) 的量级时，摄动格式的优势较为明显， 最大误差和平均误差都能降低两个数量级 2 ) ．对一维N a v i e r - S t o k e s 方程( 1 0 ) ，边界条件取u 叫5 4N S 方程激波计算的摄动有限差分方法p ( O ，f ) = u ( O ，f ) = r ( o ，f ) = 1 ，u O ，} )一2 ／。