非均匀时间步长显式龙格库塔积分方法刘黎

1、热机气动热力学编号：092080中国工程热物理学会学术会议论文非均匀时间步长显式龙格库塔积分方法刘黎1，李晓东1,胡方强2(1.北京航空航天大学流体与声学工程实验室，北京100191，中国)(2.奥多米尼大学数学与统计系，弗吉尼亚23529,美国)(Tel:E-mail: lixd)摘要：显式龙格库塔方法是数值模拟常微分方程解的一类常用单步积分方法，广泛应用于计算流体 力学等领域。传统的显式龙格库塔方法要求时间步长全场一致，对于计算流体力学、计算气动声学等 研究领域中常常遇到的涉及不规则几何形体的流动或复杂多尺度流动问题的数值模拟，往往导致昂贵 的计算资源需求。为了

2、提高数值模拟复杂流动问题的计算效率，本文发展了非均匀时间步长耦合积分 方法，使显式龙格库塔格式适用于局部步长积分。文中给岀了该耦合方法的一般线性公式，并推广至 几种常用显式龙格库塔方法，与间断有限元空间离散和一种有限差分离散相结合的数值验证均表明该 方法的有效性和稳定性。关键词：显式龙格库塔积分，非均匀时间步长，计算流体力学0前言龙格库塔方法(Runge-Kutta method, RK)是数值模拟常微分方程(Ordinary Differential Equations, ODE)解的一类常用单步积分方法，分隐式和显式两种。其中，显 式龙格库塔方法由于执行相对简单，精度较高，广泛应用于计算流

3、体力学，计算气动声 学和计算电磁学等计算科学领域，如经典 4阶RK格式10，2阶和3阶的总体变差减小 RK 格式(Total Variation Diminishing Runge-Kutta shemes, TVD RK)叫以及各种基于 稳定性范围或频散耗散特性的优化格式何13】。传统的显式龙格库塔方法要求积分步长均匀，即整个计算域内步长处处相等。对于计算流体力学(Computational Fluid Dynam-ics, CFD)和计算气动声学领域(Computational Aeroacoustics, CAA )多涉及的具有不 规则几何形状或复杂多尺度流动问题的数值模拟，出于计算时间

4、和局部要求高分辨率的考虑，数值模拟往往采用分块的非均匀网格。然而，根据时间稳定性条件(Coura nt-Fri-edrichs-Lewy Condition, CFL)，均匀时间步长受最小网格尺寸限制，导致昂贵的计算 资源要求，尤其在非定常数值模拟中；这也是数值模拟方法应用于实际复杂工程问题的 主要障碍之一。因此，为了解决这一实际应用的难题，我们希望实现非均匀时间步长显 式RK推进，即在局部网格块上采用满足局部CFL条件的局部时间步长，从而减少粗网格上积分时间，提高整场计算效率。近五年来，关于局部时间步长推进的研究工作开始出现。Tam等人基于Adams-Bas-hforth型时间积分方法在笛卡

5、尔坐标系下发展了适用于2:1比率的变时间步长时间推进格式，并与高阶有限差分离散相结合进行了验证11。此后，Garrec等人将该格式进一步基金项目：国家自然科学基金(50890181和50676003)资助项目其中，U是未知的时间相关向量， 塔时间积分是求解该方程常用方法之一。 假设tn时刻的解Un已知，tn 1 tk2him kp FUnunUn 1 Unta?ik it(apikiap2k2 川 a pp 1k p 1)ptbjkji 1推广至曲线坐标系下5。然而，该格式中各点加权系数依赖于时间步长比率，不同比率 情况需要重新推导，实际应用存在一定局限。近期，Lin, Jia ng & Li

6、提出了一种优化的 Hermit差值技术，可直接应用于任意时间步长比率的Adams型积分。Munz等人通过引入Cauchy-Kovalevskaja (CK)方法(通常也被称为 Lax-Wendrof时间离散方法)，在 利用时空偏导数的任意高阶格式(Arbitrary high order schemes using DERivatives, ADER )中实现了任意时间步长比率的局部时间推进1415。同样基于CK方法，Gassner等人发展了称为时空展开的间断伽辽金格式(the Space-Time Expa nsion Disc on ti nu ous G-alerkin scheme,

7、STE-DG)，也适用于局部步长的时间离散问。尽管，在线性问题中C-K方法被证明是高效的，对于非线性问题，CK方法必须借助Leibniz广义法则实现，计算量因此增加。最近，本文作者及合作者从一般p阶显式积分形式出发，发展了一种非均匀时间步长显式RK积分的高阶线性耦合方法，并与间断伽辽金有限元方法(Discontinuous Gal-erkin method, DG)相结合，进行了必要的数值稳定性分析和数值验证。该方法适用于任意时间比率情况，可方便的推广至各种显式RK格式，因此具有普适性。本文工作在已有基础之上，进一步将所发展非均匀时间步长积分方法推广至计算流体力学(CFD)模拟中常用的几种显式

8、 RK格式，包括经典 4阶显式RK格式10，2阶和3阶精度的T VD RK8，线性强稳定性保持格式( Lin ear Stro ng Stability Preservi ng Run ge-Kutta schemes, SSPRK linear)。此外，还与4阶频散关系保持格式相结合，探索了该耦合方 法在有限差分离散中的应用关键技术，并进行了初步数值验证。1非均匀时间步长显式龙格库塔积分本节首先将简单回顾非均匀时间步长显式RK积分方法，包括线性耦合公式及其简单推导。然后，介绍经典4阶RK格式等几种常用格式的推广。1.1线性耦合公式在数值模拟中，经过空间离散后有待时间积分的偏微分方程可简化成下

9、列形式：(1)卫FUtF表示空间算子，和空间离散格式相关。显式龙格库 这里，首先考虑最一般的p阶显式RK格式。t时刻的解un 1可通过以下公式求得：kiunCununun tinif Jrvel :I十iDtEnutrdifiie liiiit t=lmA-CfAt-2 TCoaise meshFine mesh图1非均匀网格的非均匀时间步长积分示意图为方便阐述，图1给出了一个一维非均匀有限元网格示意图。该网格由两块均匀网 格组成，网格尺寸分别为任意尺寸Xc和Xf ；下标c和f 分别代表 粗网格(coars-e mesh)和 细网格(fine mesh)。单元0和单元1是分别与网格交界毗邻的粗

10、/细网格单位，彼此相邻。若全场采用一致的 CFL数，两块网格的非均匀积分时间步长将满足关系式：tc : tfXc： Xf(3)当两块网格分别以各自的步长开展显式RK时间积分时，与网格交界毗邻的粗/细网格单元的空间离散要求相邻细/粗网格单元的信息，包括相同参数条件下的中间阶量匕和相同时间层上的未知变量U。因此，非均匀时间步长积分过程中，为了保证粗、细网格间数值信息的正确传递，单元0和单元1上的中间阶量ki及未知变量U需要进行适当的耦合处理。为了非均匀时间步长显式RK积分的执行，本研究需要考虑两种情况下的耦合关系：1)相邻粗/细网格位于相同时间层上，称为公共时间层(common time leve

11、l)”如图1中tn 1 tm 1时间层(m, n均为整数);为了使两块网格能够以各自的步长进行RK时间推进，需要通过相邻细、粗网格单元上已知量kfj、Uf和kc,i、Uc计算对应的耦合量02、Uf和kci、Uc，且必须同步进行。2)粗/细网格位于不同时间层上，称为 中间时间层(in termediate time level) ”如图1中tm和tn时间层；为了当前粗/细网格的R K时间积分，需要通过相邻细/粗网格单元上已知的kf/c,i ,Uf/c计算对应的耦合量kf/cjUf/c。这里，V表示向量v对应的耦合向量。事实上，推导结果表明，这两种 情况下得到的耦合公式满足统一形式。(9)对于线性

12、问题，数学分析发现，RK时间积分的中间阶量ki和未知量U对于时间kpJcp2cp3CppJl000p 1tpU tp t tnVV的各阶偏导数满足下列线性关系式：k110001000U tk21c22000t002U t21c32c33000t20（4）CPt其中，Cj与rk系数a.j相关。根据（4）式，不难建立情况1）所需耦合中间阶量K ki与已知量K ki之间 的关系：K f CPtcPtQKf Tf Kf（5）Kc CP tfPt1cC 1KC Tc Kc（6）其中，Tc和Tf分别为相邻粗、细网格单元所需的耦合系数矩阵。它们互为逆矩阵，均 为下三角型矩阵；这对于公共时间层上相邻粗/细网格

13、单元上耦合程序的同步执行十分重要。得到耦合中间阶量，再利用（2）式即可获得公共时间层上耦合未知量。对于中间时间层上的耦合关系式，可以借助泰勒级数展开实现。在线性问题中，中 间时间层上积分所需 K通过对已知临近时间层上的K进行p阶泰勒级数展开得到。因而，粗网格单元0在tm上的RK积分和细网格单元 1在tn上的RK积分所需要的耦合中 间量可分别通过（7）式和（8）式求得，进而得到中间时间层上耦合未知量。K f CP tcB 2Pt1fC 1K f Rf K f（7）Kc CP tf B 尸 t；C 1Kc Rc Kc（8）其中，B 1和B 2是由泰勒级数展开得到的上三角型系数矩阵，1和2则分别是中

14、间时间层与已知临近时间层的时间差，如图1中所示。Rc和Rf即为中间时间层上 RK积分必需的耦合系数。当 0时，B 为单位矩阵。那么，（7）式和（8）式将分别还原成 （5）式和（6） 式。因此，公共时间层上耦合关系式可视为中间时间层上耦合关系式的一个特例。换句 话说，（7）式和（8）式即为适用于任意时间步长比率的非均匀时间步长显式RK积分所需的一般耦合关系式。关于耦合关系式的更多细节，可参看文献4。1.2几种常用显式RK格式的推广1.1中的非均匀时间步长显式 RK积分的线性耦合公式广泛适用于各种显式 RK积分 格式，且对于不同格式区别仅仅在于（ 4）式中系数矩阵 C。在参考文献4中，已经介 绍了

15、一些高阶RK格式的应用。本文将进一步给出几种 CFD中常用显式RK格式的推广。1.2.1经典4阶显式RK格式具有4阶精度的经典显式 RK格式10可以表示成（2）式，RK系数为a21932a313411 2, 94313420b! b4 16,b2 b3 1 3因此，对应系数矩阵C可利用与（4）式中C同样的方法得到，即10 0 0(10)11/2 0 0 C1 1/2 1/401112141.2.2总体变差减小RK格式2阶和3阶精度TVD RK格式8是CFD和CAA中常用格式之一，尤其对于含有物 理间断的问题。它们通常表示为如下形式：u（0） unUi1 i u 0i U i 1tF U i 1

16、(11)i 1,lIl,mn 1mUU .其中，1为常系数。对于2阶（m2 ）和3阶(m 3)TVDRK格式，系数分别为：whe n m 2,1 1, 21 2;whe n m 2,1 1, 21 4,32 3.为了得到耦合系数矩阵C，（ 11）式被重写成和（2）式相似的一般m阶RK形式，获得对应RK系数aij ;从而求得2阶和3 阶 TVDRK的耦合系数矩阵C为：10whe nm 2,C11100(12)whe nm 3,C 110 .11 21 41.2.3线性强稳定性保持 RK格式SSPRK格式是一类具有较强稳定性的RK格式【7，即允许较大的 CFL数。尽管方法设计的出发点主要是针对非线

17、性问题，事实证明，它对于常系数的线性问题同样有效，且具有高阶精度。对于线性问题，m阶SSPRK格式在保证CFL数为1.0的前提下，数值精度最高可达到 m阶（linear SSP RK （m, m）；如果牺牲一定精度，稳定性范围还 将增大。因此该方法在CFD中应用广泛。这里，我们以线性SSPRK （m, m）格式的耦合 系数表达式推导为例，其它线性SSPRK （m, p）（ p m ）格式的耦合系数可类推得到。线性SSPRK （m, m）格式的表达式如下:0UUiUnUi1tFUi1i 1,|,m 1,(13)m,kU ktF U m1k 0Un1U m其中， m,k为常系数，满足下列关系式1,

18、00,1m,k m 1,k 1 ,k1m,m 1m!将（13）式写成一般 m阶形式，进而得到与矩阵元素为皿1，k 1,川,m 2,m 1m,04)m, kk 1式中同样形式的耦合系数矩阵Ci ,k(i1)!(k 1)!(i k)!i 2| ,m, k 1,川,i(14)2数值验证及结果分析2.1 一维迁移方程所发展非均匀时间步长显式初始条件为：RK积分方法首先用于求解一维线性问题。控制方程及uu0, x 0,100 txu x,00.(15)自t 0时刻起，左边界x 0处一正弦波u（0,t） sin（ t）开始被激励，并向内场传播。 方程的空间离散采用 4阶精度的DG空间离散2，时间积分采用4

19、阶精度的4阶S SPRK格式。为了验证所发展耦合程序的正确性，计算网格采用由 非均匀有限元网格，整场CFL数一致，因此时间步长非均匀。图3块均匀网格组成的2和图3分别为时间步长比率 t1 : t2: t3为2:1:2和3:2:3的网格和t 100时刻的计算结果，此时第 弦波刚好传出右边界，波布满整个计算区域。计算结果显示，当正弦波分别穿过网格块 交界x 40和x 60时，未见明显抖动。表明所发展耦合方法保证了不同时间步长的 网格单元间数值信息的正确传递，验证了其有效性。1个正E 0 uh(x,t) uh(x 20 ,t)| dx(16)2(,t) dh p 1.2 i 0 1此外，我们对比了

20、t 100时刻第1和第20周期的正弦波数值解（如图 4）,式（16）得到了取近似多项式阶数 P为1至4的DG离散的L2空间数值误差 E;ii 20(P 1)Uh( ,t) Uh根据公式中的图2 t : t2 : t3 = 2：1：2的网格和t= 100时刻数值结果n u43J5nsUHg83D-|Ix图3 t : t2 : t3 = 3:2:3的网格和t= 100时刻数值结果图4 t= 100时刻第1个和第20个周期波示意图图5是误差E与网格单元数 N的对数坐标图，即误差收敛曲线。图中实线是斜率 为2P 1的参考线。图5显示，当网格加密，P为1至4的SSPRK-DG离散的误差收敛 率达到2P

21、1 ;结果表明，所发展线性耦合方法较好的保持了DG方法的超收敛特性。2.2 一维 Burger 方程在验证了所发展耦合方法对线性问题的有效性后，我们进一步探索了其对非线性问题的求解能力。一维 Burger方程是典型的非线性偏微分方程，形式如下：(17)上 _u-2t2u 0XX初始条件选为一双曲线函数u (x, 0)1 tanhxX。2x 0,50其中，为粘性系数。为了能与解析解比较，Xo4,空间离散采用4阶精度的0.2.DRP有限差分离散12，时间离散则采用3阶TVD RK格 式。离散网格则是由两个子网格块组成的非均匀有限差分网格，图6、7中的点划线表示网格交界；时间步长满足t1 : t2X

22、1 ： X2 2:1。为保证两块网格的非均匀时间步长积分顺利进行，时间耦合方法将在重叠区域内插值所涉及各点上进行。t= 0、16、26和36时刻数值解U及1阶空间偏导数与相应解析解的对比结果表明，本文发展的非均 匀时间步长显式 RK线性耦合方法适用于这一非线性算例。此外，此算例也验证了所发 展方法与有限差分离散结合应用的一些关键技术，拓展了方法的应用范围。3结论本文发展了一种具有普适性的非均匀时间步长显式龙格库塔积分线性耦合方法，实现了非均匀网格上采用满足局部稳定性条件的局部时间步长积分。该方法普遍适用于各种显式龙格库塔时间积分格式，且应用于不同积分方法的区别仅在于耦合系数矩阵C。基于已有的对

23、于计算气动声学中高阶显式RK方法的应用，本文进一步给出了在计算流体力学等领域目前常用的几种显式龙格库塔格式的推广，包括经典4阶格式，线性强稳定性保持格式等。数值验证中，所发展方法首先与间断有限元空间离散结合，在两个非 均匀网格上采用非均匀时间步长求解了一维迁移方程。数值结果表明了方法的有效性。L2空间的误差分析显式了明显的超收敛特性。此外，通过求解一维非线性 Burger方程, 该耦合方法与有限差分方法结合应用的一些关键技术也得到初步验证， 为下一步应用研 究奠定了基础。图6数值解u与解析解对比图图7数值解ux与解析解对比图参考文献1 J.C. Butch, The Numerical Ana

24、lysis of Ordinary Differential Equations, Runge-Kutta and Gene ral Linear Methods, Wiley, New York, 1987.2 B. Cockburn and C.-W. Shu, TVB Runge-Kutta Local Projection Discontinuous Galerkin Finite Element Method for Conservation Laws V: Multi-Dimensional Systems, Journal of Computational Physics, 19

25、98, 141: 199-224.3 D.K. Lin, M. Jiang and X.D. Li, An Optimized Multi-Time-Step Interpolation Scheme for Overset Grids, 2009, AIAA Paper No. 2009-3118.4 L. Liu, X.D. Li and F.Q. Hu, Non-Uniform Time-Step Explicit Runge-Kutta Discontinuous Galerkin Method for Computational Aeroacoustics, 2009, AIAA P

26、aper No. 2009-3114.5 T.L. Garrec, X. Gloerfelt and C. Corre, Multi-Size-Mesh, Multi-Time-Step Algorithm for Noise Computation Around An Airfoil in Curvilinear Meshes, 2007, AIAA Paper No. 2007-3504.6 G. Gassner, F. L?rcher, and C.-D. Munz, A Discontinuou Galerkin Scheme Based on A SpaceTime Expansio

27、n. I. Inviscid Compressible Flow in One Space Dimension, Journal of Scientific Computing, 2007, 32: 175-199.7 S. Gottlieb, On High Order Strong Stability Preserving Runge-Kutta and Multi Step Time Discretization, Journal of Scientifics Computation, 2005, 25: 105-128.8 S. Gottlieb and C.-W. Shu, Tota

28、l Variation Diminishing Runge-Kutta Schemes, Mathematics of Computation, 1998, 67(221): 73-85.9 F.Q. Hu, M.Y. Hussaini, and J.L. Manthey, Low-Dissipation and Low-Dispersion Runge-Kutta Schemes for Computational Acoustics, Journal of Computational Physics, 1996, 124: 177-191.10 A. Jameson, W. Schmidt, and E. Turkel, Numerical Solutions of The Euler Equations By The Finite Volume Methods Using Run