《龙格库塔间断有限元方法在计算爆轰问题中的应用方案》由会员分享，可在线阅读，更多相关《龙格库塔间断有限元方法在计算爆轰问题中的应用方案（13页珍藏版）》请在金锄头文库上搜索。
1、龙格库塔间断有限元方法在计算爆轰问题中的应用摘要 构造了求解带源项守恒律方程组的龙格库塔间断有限元（RKDG）方法，并分别结合源项的Strang分裂法和无分裂法数值求解模型守恒律方程和反应欧拉方程。为了和有限体积型WENO方法进行比较，设计了计算源项的WENO重构格式。对一维带源项守恒律的计算表明，对于非刚性问题，RKDG方法比有限体积型FVWENO方法的误差更小，而对于刚性问题，RKDG方法对于间断面位置的捕捉更为精确。对于一二维爆轰波问题的计算结果表明，RKDG方法对于爆轰波结构的分辨和爆轰波位置的捕捉能力更强。关键词 龙格库塔间断有限元方法；爆轰波；反应Euler方程；刚性源项引言在非平
2、衡气体动力学中，爆轰问题通常被描述成无粘性的化学反应流动问题，其物理模型是一个非齐次双曲守恒律方程组，通常被称作反应欧拉方程组（reactive Euler equations），其中非齐次项（源项）通常被解释为由于化学反应引起的混合组分的质量变化率1。最简单的反应欧拉方程组假定混合气体仅仅由两种组分构成：已燃气体（burnt gas）和未燃气体（unburnt gas）2。当混合气体达到点火温度时，未燃气体通过一个不可逆的化学反应转化为已燃气体，因此，混合气体状态可以用一个标量，即未燃气体的质量分数Y来表示，进一步假设混合气体各组分具有相同的比热比和比气体常数R，则二维反应欧拉方程组可以写成
3、如下形式 , （1）其中，各个变量和函数的具体形式为式中是混合气体的密度，分别是沿方向的速度分量，分别是压力，总能量，温度。称为化学反应率，通常用Arrhenius模型 （2）其中是化学反应率因子，是活化能。此外混合气体的状态方程和温度分别为 （3）其中为单位质量未燃气体发生化学反应所释放的热量。在实际爆轰波问题中，化学反应的时间尺度远小于流体流动的时间尺度，所以非齐次方程组具有很强的刚性，给数值求解带来很大困难。无论对源项用不用算子分裂方法，采用显格式还是隐格式，在处理间断解问题时，都可能得到非物理解。这是因为当源项不能用足够的空间和时间分辨率进行求解时，激波捕捉方法会得到错误的爆轰波传播速
4、度3,4,5。文5发现，当空间分辨率不够时，一个同时含有已燃气体和未燃气体的网格中会发生虚假反应，数值爆轰波会以每个时间步一个空间步长的速度传播。文6也指出计算的数值误差极可能促成温度敏感化学反应的提早发生。因此，要获得正确的爆轰波位置，常用的解决方案是采用合适的点火温度模型来抑制虚假反应3，或者采用网格自适应方法7, 8, 9来保证反应区域内有足够多的网格点，或者采用高精度高分辨率方法来保证对反应区有足够高的分辨率。由于后两种途径有更好的普适性而受到更多的关注。由Cockburn和Shu等人发展的龙格库塔间断有限元方法（Runge-Kutta Discontinuous Galerkin M
5、ethod，RKDG）是一类具有高精度和高分辨率的数值方法10，在解决含有间断现象的问题中发挥着越来越重要的作用，它被广泛地发展和应用于水动力学，气动力学，波传播，半导体中的电荷传输等问题。该方法既保持了一般有限元方法和有限体积方法的优点，又克服了各自的不足。该方法可采用局部高阶插值的方法构造基函数，具有灵活处理间断和边界条件以及可显式求解的能力，克服了一般有限元方法不适于间断问题的缺点，以及一般有限体积方法必须通过扩大模板进行重构来提高精度的不足。RKDG方法所具有的单元上连续分布的高精度高分辨率逼近特性有望能更好地模拟爆轰波问题。本文应用RKDG方法，结合化学源项的无分裂方法及算子分裂方法
6、，对爆轰波问题进行数值模拟，并和有限体积型WENO（FVWENO）方法的计算结果进行比较。由于FVWENO方法需要利用流动变量单元平均值来计算源项的单元积分，本文采用Simpson积分公式，其积分点处的流动变量采用WENO重构以提高计算精度。对于带源项的一维守恒律的数值测试结果表明，对于非刚性问题，RKDG方法计算结果的误差更小；而对于刚性问题，RKDG方法对于间断波位置的捕捉能力更强。对典型的一二维爆轰波算例，包括二维不稳定爆轰波和爆轰波绕过90拐角的衍射问题的数值模拟结果显示RKDG法仍有一定优势。1 数值方法1.1 空间离散本文的空间离散采用RKDG方法10。设是区域的一个有限剖分。单元
7、表示多边形单元的一条边界，表示单元边界的外法向。是单元上的局部有限元空间，取作次多项式集合。，在间断有限元空间中寻找近似解，其中。首先在单元上用试探函数乘以方程（1）的两端，并用近似解代替方程（1）的精确解，用代替试探函数，将（1）写成变分形式，并由分部积分和Green公式，可得 （4）对于（4）式的后三项，采用数值积分计算 （5） （6） （7）流通量采用与之相容的数值流通量代替，数值流通量定义为为某种形式的Riemann解算器。分别代表边积分和单元积分的Gauss积分点个数，和分别代表边积分和单元积分的Gauss积分点位置10。这样得到半离散的数值格式 （8）为了计算方便，在单元中取正交基
8、函数（如勒让德多项式），则质量矩阵成为对角矩阵，故有限元解表示为 （9）在式（8）中，取试探函数为所有基函数，得到半离散的常微分方程组 （10）其中，为单元的质量矩阵，和分别为 （11）和 （12）高阶RKDG方法会产生数值振荡，需要在每一个时间步后对数值解应用斜率限制器，本文采用文献10的TVB限制器。当基函数的阶数增加时，流通量的选择对于数值结果的影响较小10，所以本文采用局部的Lax-Friedrichs流通量。1.2 无分裂方法和分裂方法对于空间离散后得到的常微分方程组（10），如果源项的刚性不太强，可以直接采用常微分方程组的求解方法进行时间推进。本文的无分裂方法是采用三阶强稳定保持的
9、龙格库塔（SSPRK）方法11直接求解方程组（10），该方法在空间离散取阶基函数时，可以保证光滑解的三阶时间和空间精度。但当源项的刚性比较强时，通常要对流动项和化学源项分裂求解，即对于方程组（10）分别求解流动项的半离散方程 （13）和化学源项的半离散方程 （14）本文采用的分裂法是具有二阶时间精度的Strang分裂方法12 （15）其中和分别表示以和为时间步长求解（14）和（13）的算子。在分裂方法中，守恒律的半离散方程组（13）仍采用显式的三阶的SSPRK方法求解，而源项的方程组（14）采用常微分方程求解器VODE求解13。1.3 有限体积型WENO方法本文将比较RKDG方法和高精度的FV
10、WENO方法14。在1.1节介绍的空间离散中，如果基函数只取为1，则可得到守恒律方程组（1）的有限体积方法。在二维矩形网格上，半离散有限体积法的方程为 （16）上式中对流项的离散采用维数分裂的基于特征变量的五阶WENO格式14。在应用无分裂方法计算时，由于FVWENO方法的待求变量是流动变量的网格单元平均值，如果直接用它计算（16）式中源项的积分，则只有二阶空间精度。为尽量提高精度，本文对一维问题采用四阶的Simpson公式来计算积分 （17）其中积分点的值和采用基于特征变量的五阶WENO重构得到。对于的重构会产生负权，故需要对重构结果进行处理15。对于二维问题，可以采用（17）式的张量积形式
11、，这时积分具有四阶精度，但涉及重构积分点处值，计算公式相对复杂。为计算方便，本文采用一种具有三阶精度的二维Simpson积分 （18）其中由方向的单元平均值重构得到，由方向的单元平均值重构得到，采用两个方向重构值的平均值。2 数值算例2.1 精度测试第一个例子是一个带有源项的Burgers方程 （19）这个方程的解是连续的，其解析解为 （20）式中控制方程的刚性，即为初值。为了测试算法的精度，取，此时问题的刚性不强。计算区间取为-15,25，RKDG方法中TVB限制器的常数取为5，计算到0，CFL数取为0.18。表1是数值解的误差及其数值精度阶。表1 带有源项的Burgers方程Table 1
12、 Burgers equation with source termNStrang split methodUnsplit methodDG2RWENO5RDG2RWENO5RWENOMR327.6610-54.7910-37.5510-54.7810-35.2510-3645.2910-63.869.6210-42.314.7710-63.989.6410-42.315.0410-43.381284.5810-73.532.4210-41.991.9110-74.642.4310-41.993.1910-53.982568.3710-82.455.7610-52.071.3610-83.81
13、5.7610-52.072.4110-63.725121.8410-82.181.4310-52.011.4210-93.261.4310-52.012.8010-73.1110244.3310-92.093.5610-62.001.6910-103.073.5610-62.003.5310-82.9920481.0510-92.048.9010-72.002.0910-113.028.9010-72.004.4610-92.99注：DG2表示P2基函数的RKDG方法，WENO5表示五阶FVWENO方法，WENOM表示五阶FVWENO方法+源项四阶精度积分公式（17）从表1中可以看出，用Str
14、ang分裂方法，由于在时间上只有二阶精度，在时间步长与空间步长同阶的情况下，即使采用高精度的基函数DG方法也只能达到二阶精度。而无分裂的元DG方法和修正FVWENOM方法可以达到和SSPRK方法同阶的三阶精度，但没有采用高精度源项积分的FVWENO只达到二阶精度，无论是分裂方法还是无分裂方法，RKDG方法在网格规模相同的情况下，其数值解的误差明显要比FVWENO方法的小得多。2.2 间断计算能力测试第二个例子是一个带有源项的对流方程 （21） 取初值为 （22）是间断的初始位置，这个例子的解析解是一个以速度1向右传播的间断解。文5指出，当网格数不够多时，计算的间断波速度变慢甚至是不动的。为了衡量方法的间断捕捉能力，定义数值平均
