《全矩阵无积分间断有限元方法实现.doc》由会员分享，可在线阅读，更多相关《全矩阵无积分间断有限元方法实现.doc（11页珍藏版）》请在三一办公上搜索。
1、精品论文全矩阵无积分间断有限元方法实现谭勤学，任静，蒋洪德（清华大学热能工程系）5摘要：本文主要针对间断有限元方法计算量及存储量开销巨大的问题，提出一种节约计算量及存储量的实现方法。使用该方法，在相同自由度条件下，间断有限元方法计算速度及存储 量与普通有限体积方法相当。另外结合有限体积隐式求解方法，实现间断有限元方法精度多 重网格，进一步有效加快间断有限元方法收敛速度，改善收敛特性。 关键词：间断有限元，精度多重网格，无积分, 矩阵运算10中图分类号：V231.1Quadrature-free implementation of Discontinuous GalerkinMethod Usi
2、ng Matrix Product FormulaTan Qinxue, Ren Jing, Jiang Hongde15(Tsinghua University, Beijing 100084)Abstract: In this paper, a discontinuous galerkin only include matrix operation formulation that avoids the use of discrete quadrature formulas is described. Using this method, the memory requirements a
3、nd computing speed are almost equal to the finite volume method When degrees of freedom of the same. A p-multigrid discontinuous Galerkin method that combine explicit DGM and implicit FVM is20presented for the solution of thecompressible Euler equations on unstructured grids. And the computation can
4、 converge much better and faster with this p-multigrid method.Key words: Discontinuous Galerkin Method, p-multigrid, Quadrature Free, Matrix Operation0引言25间断有限元方法(DGM, Discontinuous Galerkin Method)最早是由 Reed 和 Hill 在 1973 年 1提出的，用于中子输运方程的求解。该方法结合了有限体积方法(FVM)和有限元方法(FEM) 的基本思想，特别易于处理复杂边界及边值问题，同时具有灵活处理
5、间断的能力，且可以通 过适当选取基函数，提高单元插值多项式的次数来实现任意高阶精度，另外具有容易实现自 适应、并行计算等优势2。由于这些优异性质，在欧洲探讨用于未来工业设计高阶 CFD 程30序的 ADIGMA 计划中3，间断有限元方法被认为是最有前途的算法之一。近年来 DGM 得到了非常大的发展，其中以 Cockburn&Shu 提出的龙格库塔间断有限元 方法尤为突出，得到了广泛的应用。Cockburn 和 Shu3-7引入龙格库塔 DGM 求解 Euler 方程。 Bassi&Rebay8将 DGM 用于求解可压 NS 方程，在求解过程中将混合有限元的思想引入 DGM 中。Cockburn
6、&Shu9将 Bassi 和 Rebay 处理梯度的方法和龙格库塔 DGM 结合提出了35LDG(Local Discontinuous Galerkin)方法，用于求解可压缩 NS 方程。在文献10，Cockburn 和Shu 综述了使用 LDG 求解不可压问题。同时，现阶段 DGM 也有很多问题需要解决。其中准确有效，鲁棒性好的限制器是 DGM面临的重大挑战之一：Cockburn3-7中提出了修正 MinMod 限制器的 TVB 限制器，Luo11使基金项目：高校新教师基金(20090002120033)作者简介：谭勤学(1986-),男，博士研究生，燃气轮机二次空气系统及高精度算法通信联
7、系人：任静,女，教授，燃气轮机透平传热. E-mail: renj- 11 -用了 Barth&Jesperson 的限制器并结合使用 Krivodonova12的间断探测器来判断光滑区。另外40DGM 较大的计算量及存储量开销，是影响 DGM 实用化的重要限制因素之一：DGM 计算中 使用 k 阶多项式进行插值，三维情况下插值多项式的数量为(k+1)(k+2)(k+3)/6，此时体积分 要求 2k+1 阶精度，面积分至少要求 2k 阶精度，在使用龙格库塔求解时要求 CFL1.0/(2k+1)， 这些因素致使 DGM 的计算量和存储量都非常庞大。现阶段，p 多重网格方法11，无矩阵隐 式方法1
8、3，无积分方法14，h 多重网格方法15等方法应用于 DGM 中来加快收敛速度和减少45存储量。e本文主要针对 DGM 计算量及存储量开销巨大的问题，提出一种节约计算量及存储量的 实现方法，以期实现在相同自由度条件下，DGM 计算量及存储量与普通有限体积方法相当。 另外结合有限体积方法，实现间断有限元方法 p 型多重网格，能有效地加快计算收敛速度。1控制方程50守恒形式的三维 Euler 方程：tU + F (U ) = 0(1.1) u v w u1 uu + p vu wuU = u2 , Fx = uv , Fy = vv + p , Fz = wv (1.2) u3 uw vw ww
9、+ p E Hu Hu Hu此时为封闭方程还需要引入状态方程等关系式：p = (P,T ), C p = C p (P,T ), H = C p dT , E = H 552间断有限元方法(1.3)对方程(1.1)，将方程两端与测试函数 相乘，在任意单元上分部积分，可得到控制方程的变分形式为：U d +F (U ) r dS F (U )d = 0(2.1)ee t其中：e 表示单元区域，e 表示单元边界，使用有限元空间近似函数U h , h 代替函p60数U , ，对于间断有限元取数值解空间和测试函数空间相同，且都取为 p 阶多项式函数空间 P p ：U h =i = Ni =1i = NU
10、i ( x, t)i ,i Pn(2.2)hiii = ( x, t) , P pi =1DGM 离散过程中可直接取 h = i(2.3)652.1数值积分方程(2.1)中的面积分和体积分常用高斯数值积分计算：如图 1 将积分单元映射到标准单元，然后在标准单元上利用标准单元上的高斯积分点进行加权求和，其积分公式如下：nqpe d = refJ d ref= i ( , , ) Ji ( , , ) Wi ( , , )i =1nnn(2.4)J ( , ) = ( x, y, z) , x =( , , )i =1xii , y =i =1yii , z =i =1zii(2.5)70N 表示
11、积分点数量，Wi 分别数值积分点权重， Ji插值形函数，同有限元的等参单元。表示雅克比行列式，i 表示坐标的752.2形函数的选取图 1.六面体单元与标准六面体单元之间的映射形函数i 可以与具体物理单元相关，如下面所示的形函数：1 = 1,2 = ( x x0 ),3 = ( y y0 ),4 = ( z z0 ) = ( x x )2 ,= ( y y )2 ,= ( z z )2 ,.(2.6)506070形函数i 也可以为参考单元坐标的函数，其只于单元类型如四面体，六面体等相关，而80与具体物理单元相关，如下面所示的形函数1 = 1,2 = ,3 = ,4 = = 2 ,= 2 ,= 2
12、 ,= ,= ,= , .(2.7)567567对于上述两种形式的形函数在使用之前均可以使用式(2.8)或式(2.9)对其进行正交化，以 保证U1 代表单元平均值%ei = i i d , i = 2, 3,.(2.8)ii85%= refi d ref , i = 2, 3,.(2.9)3全矩阵运算无积分 DGM 实现方法Shu14中提到一种节约内存及加快计算速度的实现方法，但其做了雅克比矩阵在单元内 各积分点为常量的假设，并对数值通量进行重构。本文推导一种全矩阵运算无积分的精品论文DGM 实现方法，且不再需要雅克比矩阵为常数的假设，需要注意：在下文推导中将使用基90于参考单元的形函数。3.
13、1时间项离散时间项导数单元的积分，结合式(2.2)及式(2.4)，式(2.1)中时间项如下计算：Uj = Nq i = N U ( x, t )k i , jj j ,kje th d = ji =1i W Jt(3.1)同式(2.4)，对于标准单元上的积分有式(2.8)Urefj = Nq i = N U ( x, t) i , jj j ,k95refh d =tji =1i W t(3.2)由于其为虚拟时间积分，解的精度仅依赖与空间项的离散，因此对上述积分可近似假设单元各积分点的雅克比行列式为一近常数 J j J ，如用单元平均雅克比行列式等，而不会影响到计算的精度，结合式(3.1)，式
14、(3.2)，此时上述积分可写成如下所示的矩阵形式： U1 UU t d J d ref= J M . , Mi , j =i j d ref(3.3)hhetreftref =U nsf t100其中 Nq 表示数值积分点数量，而 N 表示形函数数量， J 表示近似平均雅克比行列式，其中矩阵 M 仅依赖参考单元，与具体物理单元无关。3.2单元体积分结合式(2.2)及式(2.4)，可得单元体积分过程： F (U )d =( ) J 1 J F (U )d erefref F F F x ,1 y ,1 z ,1 (3.4) A . + B . + C . Fx ,nqp Fy ,nqp z ,n
15、qp A=i , j W , F F= JF+ F+ F i , j jx xxyyzz 105B=i , j W , F = JF+ F+ F (3.5)i , jjy xxyyzz C=i , j W , F = JF+ F+ F i , jjz xxyyzz 其中i = 1,., N , j = 1, ., Nq , Nq 表示体数值点数量。3.3单元面积分结合(2.2)及式(2.4)，可得单元面积分过程： Fn,1 r e F (U ) ndS D . , Di , j = i , jW j(3.6)110 Fn,nsp 其中i = 1,., N , j = 1,., Ns , Ns
16、表示面数值积分点数量， Fn, j 表示第 j 个面数值积分点的法向通量。3.4数值离散过程小结总结上述离散过程，间断有限元方法离散后方程可写为下述所示半离散形式： U1 F F F F tNn,1 i x ,1 y ,1 z,1 J . +M 1 D . = M 1 A . + M 1B . + M 1C . (3.7) U i =1 F F F F = nsf n,nsp ix ,nqpy ,nqpz,nqpt115120分析比较上述离散过程对计算效率和内存开销方面的影响：1.计算效率方面，问中所述实现方法将积分过程以矩阵与向量的乘积，相关矩阵可以预 先计算好，程序编写简单，计算量减少，并
17、可以利用一些高效矩阵运算库，计算效率能有较 大提高。2.内存代价方面，由于单元类型非常有限，矩阵 M，A，B，C，D 内存开销可以忽略。 另外此实现方法仅需要存储单元内各积分点的雅克比矩阵，以及面积分点的雅克比矩阵行列 式，而不需要其它开销，内存开销和同自由度有限体积方法一致。而若使用式(2.6)类型的具 体物理单元相关的形函数，一般需要存储形函数，形函数的导数，而这部分内存开销如式(3.2) 计算，开销较大125130135Nface4 N Nq + N Nsi =1另外若使用式(2.6)类型形函数，若仍然使用上述推导所示，将积分过程以矩阵与向量的乘积替代，此时每个单元的 M，A，B，C，D
18、 都各不相同，其内存开销将非常庞大。3.经本文程序测试，使用该实现方法的高精度间断有限元方法运算速度约为 42 万自由 度/秒，而本文所编写有限体积自由度约为 60 万自由度/秒。4精度多重网格方法显示龙格库塔间断有限元方法，由于稳定性要求其对 CFL 数的限制相当苛刻，而高精 度隐式算法对内存需求非常庞大。为进一步加快高精度 DGM 的收敛速度，并兼顾内存开销， 本文结合有限体积方法和间断有限元方法，在 p 型多重网格的计算过程中，低阶上求解采用 有限体积方法。其原因在于：在于 DGM 一阶精度相同自由度条件下，有限体积方法能达到 二阶精度，并且有限体积方法能非常方便的使用隐式解法。在本文中
19、所用 p 型多重网格为两 层的 V 型多重网格，低阶有限体积使用 LU-SGS 方法求解，其具体操作过程如下所示：Lp1)在高阶精度上，用显示龙格库塔间断有限元方法推进一个时间步获得初始值U 0p2)计算高阶精度残差 R(U 0 )0pL3)将U 0， R p (U p) 限制到低阶上得到 R0，U0，具体限制过程如下：pLLppLLpp1pU0 = I p U 0= (U0) ， R0= I%pR(U 0) = R(U0 )1140其中(U 0) 代表单元平均值的第一项，而 R(U 0)1 =F 1d ， = 1.0 ，即为残e11p差 R(U 0) 的第一项LLL4)计算低阶计算源项： F
20、= R(U 0) R05)使用有限体积方法离散，隐式 LU-SGS 方法求解： R(U L ) = F LpppLLp6)修正高阶结果：U= U0+ I (U U 0) ，具体沿拓过程 I 如下：()( 0 )0145U p= U p+ U L U L150111(U p ) 为代表单元平均值的第一项5算例与分析本文中测试算例：扰 NACA0012 翼型无粘流动，远场流 Ma 数 0.63，进口攻角 2 度 ，所 用网格如图 2 所示，其中网格数 20700，全三棱柱非结构网格：图 3 分别给出了是否使用精度多重网格收敛过程，图 4、图 5 分别给出了计算时间与迭 代步数与计算残差的情况。对比
21、是否使用 p 多重网格的最终收敛结果，最终结果一致。另外 对比收敛情况，如图 4 及图 5 所示，使用 p 多重网格的龙格库塔间断有限元方法，其收敛特 性有所改善，且收敛速度有明显的提高。155图 2 计算所用网格(右：局部放大图)(a)三步龙格库塔求解160(b)三步龙格库塔+p 型多重网格图 3 二阶精度 HLL 格式计算所得到的 Ma 数分布云图图 4 收敛步数对比1656结论图 5 收敛时间对比本文主要针对间断有限元方法计算量及存储量开销巨大的问题，提出一种全矩阵无积分 DGM 实现方法。该方法使用基于参考单元的形函数，通过矩阵运算代替积分过程，可以有 效节约计算量，提高计算效率。经本
22、文程序测试，在相同自由度条件下，其内存消耗及有限170175180185190195200205体积方法差别不大，计算速度约为有限体积方法的 2/3。另外结合有限体积隐式求解方法，实现间断有限元方法 p 多重网格：高阶间断有限元求 解计算采用 3 步显示龙格库塔方法，低阶精度的计算采用隐式 LU-SGS 有限体积求解，经算 例验证能有效加快间断有限元方法收敛速度，改善收敛特性。参考文献 (References)1 Reed, W.H. and T.R. Hill, TRIANGULAR MESH METHODS FOR THE NEUTRON TRANSPORT EQUATION.2 刘儒勋与
23、舒其望著, 计算流体力学的若干新方法. 2003, 北京: 科学出版社. 244.3 COCKBURN, B. and C.W. SHU, The Runge-Kutta local projection P1-discontinuous Galerkin finite element method for scalar conservation laws. 1988.4 COCKBURN, B. and C.W. SHU, TVB RUNGE-KUTTA LOCAL PROJECTION DISCONTINUOUS GALERKIN FINITE-ELEMENT METHOD FOR CON
24、SERVATION-LAWS.2.GENERAL FRAMEWORK. MATHEMATICS OF COMPUTATION, 1989. 52(186): p. 411-435.5 COCKBURN, B., S.Y. LIN and C.W. SHU, TVB RUNGE-KUTTA LOCAL PROJECTION DISCONTINUOUSGALERKIN FINITE-ELEMENT METHOD FOR CONSERVATION-LAWS .3. ONE-DIMENSIONAL SYSTEMS. JOURNAL OF COMPUTATIONAL PHYSICS, 1989. 84(
25、1): p. 90-113.6 COCKBURN, B., S.Y. LIN and C.W. SHU, THE RUNGE-KUTTAL OCAL PROJECTION DISCONTINUOUS GALERKIN FINITE ELEMENT METHOD FOR CONSERVATION LAWS IV: THE MULTIDIMENSIONAL CASE. mathematics of computation, 1990. 54(190): p. 545-581.7 Cockburn, B. and C. Shu, The Runge-Kutta Discontinuous Galer
26、kin Method for Conservation Laws V: Multidimensional Systems. JOURNAL OF COMPUTATIONAL PHYSICS, 1998(141): p. 199-224.8 Bassi, F. and S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations. JOURNAL OF COMPUTATIONAL P
27、HYSICS, 1997. 131: p.267-279.9 Cockburn, B. and C.W. Shu, The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM JOURNAL ON NUMERICAL ANALYSIS, 1998. 35(6): p. 2440-2463.10 Cockburn, B. and C.W. Shu, Runge-Kutta discontinuous Galerkin methods for convection-dom
28、inated problems. Journal of Scientific Computing, 2001. 16(3): p. 173-261.11 Luo, H., J.D. Baum and R. Lhner, Fast p-Multigrid Discontinuous Galerkin Method for Compressible Flows at All Speeds. AIAA Journal, 2008. 46(3): p. 635-652.12 Burbeau, A., et al., A Problem-Independent Limiter for High-Orde
29、r Runge-Kutta Discontinuous GalerkinMethods. Journal of Computational Physics. 169(1): p. 111-150.13 Crivellini, A., et al., An implicit matrix-free Discontinuous Galerkin solver for viscous and turbulent aerodynamic simulations. COMPUTERS & FLUIDS. 50(1): p. 81-9314 Atkins, H.L., H.L. Atkins and C.
30、 Shu, Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations. AIAA Journal. 36: p. 775-78215 Klaij, C.M., et al., h-Multigrid for space-time discontinuous Galerkin discretizations of the compressibleNavierStokes equations. Journal of Computational Physics. 227(2): p. 1024-1045.
