张荣培 王迪 蔚喜军 温学兵†

1) (沈阳师范大学数学与系统科学学院,沈阳 110034)

2) (北京应用物理与计算数学研究所,北京 100088)

波的传播往往在复杂的地质结构中进行,如何有效地求解非均匀介质中的波动方程一直是研究的热点.本文将局部间断 Galekin(local discontinuous Galerkin,LDG)方法引入到数值求解波动方程中.首先引入辅助变量,将二阶波动方程写成一阶偏微分方程组,然后对相应的线性化波动方程和伴随方程构造间断Galerkin格式;为了保证离散格式满足能量守恒,在单元边界上选取广义交替数值通量,理论证明该方法满足能量守恒性.在时间离散上,采用指数积分因子方法,为了提高计算效率,应用Krylov子空间方法近似指数矩阵与向量的乘积.数值实验中给出了带有精确解的算例,验证了LDG方法的数值精度和能量守恒性;此外,也考虑了非均匀介质和复杂计算区域的计算,结果表明LDG方法适合模拟具有复杂结构和多尺度结构介质中的传播.

1 引 言

波的传播是能量传输的一种基本形式,出现在许多科学、工程和工业领域.波动传播问题的研究对地球科学、石油工程、电信和国防工业具有重要意义[1−4].波动方程是描述波传播的数学模型,例如声波方程、地震波方程等.本文考虑如下二维非均匀介质下的二阶波动方程.由于所求解的区域比较复杂,以及传播介质的复杂和不均匀性,大多数波动方程不能精确求解,对于波动方程数值方法的研究就显得非常重要.

以往二维波动方程数值解的计算方法主要是有限差分方法 (finite difference method,FDM)[5,6]、有限体积方法 (finite volume method,FVM)[7]和有限元方法 (finite element method,FEM)[8]等.局部间断 Galerkin (local discontinuous Galerkin,LDG)方法是Cockburn和舒其望[9]在1998年提出的.与传统的间断Galerkin (discontinuous Galerkin,DG)方法[10]比较,LDG 方法通过引入辅助变量将含有高阶导数的微分方程写成只含有一阶导数的偏微分方程组,然后用DG方法进行空间离散.Chung和Engquist[11]在交错三角网格上提出了能量守恒的DG格式,但是交错网格在高维情形格式构造复杂.Chou等[12]发展了一类能量守恒的LDG方法求解二维非均匀介质中的线性二阶波动方程,但是他们的工作只限于矩形网格.本文将在三角网格上构造能量守恒的LDG格式.首先将二阶标量波动方程写成一阶速度-应力形式,通过选取广义交替数值通量,对相应的线性化波动方程和伴随方程构造DG格式.对于高维问题,空间离散后得到的是一组规模非常大的常微分方程组.显式时间离散格式不需要求解代数方程组,但是时间步长有严格的限制;隐式时间离散允许比较大的时间步长,但是需要求解大规模代数方程组.本文采用指数积分因子方法求解常微分方程组,该方法是在2013年由Nie等[13]提出的求解刚性常微分方程组的有效数值方法.隐式积分因子方法的发展分为两个方向:其一是紧致隐积分因子格式[14,15],该格式将解写成矩阵形式并在每个方向计算指数矩阵;其二是Krylov积分因子格式[16],该格式针对指数矩阵与向量的乘积运算,应用Krylov子空间方法近似.

2 局部间断有限元方法

在二维区域Ω求解如下二阶波动方程:

其中,u 是位移;ρ 是物体的密度;A 是波速,本文取常数;f (x,y,t) 是外部源力项.引入两个辅助变量:p=A∇u,w=ut,(1)式可写成如下的等价一阶偏微分方程组:

边界条件考虑齐次 Dirichlet边界:u=0,(x,y)∈∂Ω×[0,T].所对应的初始条件为 :u(x,y,0)=u0(x,y),p (x,y,0)=p0(x,y),∀ (x,y)∈ Ω.当f=0时,系统(2)和(3)式具有能量守恒性:

下面针对系统(2)和(3)式构造LDG格式.首先将二维计算域Ω离散为有限个三角形单元用表示的所有边界的集合,其中是所有内部边缘的集合.定义LDG近似空间为

矢量函数q的平均值和跳跃项可以类似定义为

要注意的是,标量函数v的跳跃项 [v]是与法线平行的向量,向量函数q的跳跃项 [q]是标量.由于讨论的是齐次 Dirichlet边界条件,即在∂Ω上时,u=0.因此对于这种边界条件,将[v]和{q} 设为[v]=vn,{q}=q,其中n是指向外侧的单位法向量.

2.1 LDG方法

在(2)式和(3)式两侧同时乘以试探函数vh,qh,并在每个单元上进行积分,通过分部积分可以得到问题(2)式和(3)式的LDG格式,即找到wh∈Vh,ph∈Σh,使得任意试探函数vh,q,对所有的K∈ Γh满足

2.2 能量守恒性

在构建二维波动方程的数值方法时,能量守恒通常被考虑在内.数值方法是否能够保持能量守恒可以判断该方法是否能更好地模拟长时间的波传播.本节证明LDG方法可以保持能量守恒.考虑到数值通量(10)式和(11)式的定义,在所有单元上对(8)式和(9)式求和可得

对(12)式和(13)式左右两端分别求和,可得

为了得到(12)式和(13)式的能量守恒性,首先证明如下引理.

引理对于所有的试探函数

是恒成立的.

证明对(15)式左端项进行分部积分可得

利用平均值和跳跃项的定义,把各个单元的总和改写为边界形式,通过简单的计算可得到

将数值通量(10)式和(11)式代入(15)式右端项得到

即(15)式成立.

利用引理可以得到(8)式和(9)式的能量守恒性.

定理 (能量守恒性)在齐次Dirichlet边界条件和数值通量(10)式和(11)式的定义下,系统(8)式和(9)式是能量守恒的,即当 f=0 时,

证明选取试探函数vh=wh代入(8)式中可得到

在 (9)式中,选取试探函数qh=ph,能够得到

对(19)式和(20)式左右两边分别求和,并考虑到(15)式,可得到

3 时间离散

二维波动方程通过LDG方法进行空间离散后,与一维二阶波动方相类似,也得到一组非线性常微分方程组,为了节约这种复杂代数方程组求解的计算成本,本文利用指数积分因子方法来进行时间离散.

首先将LDG格式(8)式和(9)式对于所有单元联立,可得到全局非线性常微分方程组的矩阵方程形式,即

其中,W=(w1,w2,···,wJ)T和P=(p1,p2,···pJ)是所有单元上的自由度,wj和pj表示在第j个单元上的自由度,M1是由 3×3 矩阵组成的对角分块矩阵,M2为6×6矩阵组成的对角分块矩阵,其逆矩阵容易求解.将(22)式和(23)式进一步写成如下形式:

在方程(24)式两端同乘积分因子e-Bt,关于时间 tn到tn+1进行积分,可得

采用梯形积分公式计算(25)式中的积分,得到二阶指数积分因子格式:

对于高维情况,矩阵指数的计算将遇到巨大的困难,因为指数矩阵是大而稠密的.对于这种困难,可以通过使用Krylov子空间方法近似指数矩阵和向量的乘积来解决[16].

4 数值算例

首先给出带有精确解的波动方程测试方法的精度和能量守恒性,然后应用LDG方法求解波的传播问题.第二个算例是求解非均匀介质介质中波的传播问题,第三个算例是L形区域中波的传播问题.在下面的计算中,考虑线性元,边界条件为齐次 Dirichlet 边界,A=I 是单位矩阵,f(x,y,t)=0.

算例1考虑具有常系数的二维波动方程utt=∇2u,(x,y)∈[0,2]×[0,2],初始条件为u(x,y,0)=sin(πx)sin(πy),ut(x,y,0)=0,这个问题的精确解为时间步长取Δt=0.001,表1列出了网格数分别为16×16,32×32,64×648和1 28×12 时数值解 wh和p 的误差和收敛阶数,通过表格可以发现w的精度接近2,p的精度接近1,通常观察到这种数值精度的振荡行为就可以说明所设计的LDG方法是能量守恒的.

算例2考虑波动方程(1)在单位正方形Ω=[0,1]2上的传播,将其剖分为64×64个均匀大小的正方形,然后使用一个对角线将每个正方形细分为两个三角形,如图1(a)所示.本算例考虑非均匀介质,密度函数ρ定义为

表1 数值解 wh和p 的误差和收敛阶数Table 1.Error and convergence order of numerical solution w h and p.

图1 (a) 算例 2 的网格剖分和数值解 w h 在不同时刻 (b) t=0.2,(c) t=0.3,(d) t=0.4 时的波传播Fig.1.(a) The triangulation mesh of Example 2;contour plot of solution w h at different time:(b) t=0.2 ;(c) t=0.3 ;(d) t=0.4.

图2 算例 2 的能量随时间的演化Fig.2.Energy evolution with time for Example 2.

初始条件w0(x,y)=2exp[-500(x-0.5)2+(y-0.5)2],p0(x,y)=0.图2 给出了数值解 wh在时间 t=0.2,0.3,0.4 时的等高线,可以看出,波大约在t=0.3时触及界面 x=0.65,在那之后,波以更快的速度通过界面.这说明传播系数不同会导致波的传播速度不同.同时也给出离散能量的时间演变图(图2).从图2中可以发现,LDG方法可以保证能量守恒.

算例3考虑圆形波在L形区域Ω=[0,1]2[0.7,1]2上的传播,初始条件同算例2.采用非一致三角剖分将其剖分为13578个三角单元,如图3(a)所示.数值解 wh在 t=0.1,0.3和0.45 的波形传播图在图3(b)—(d)给出.从图3可以看出,波在 t=0.3 时刻到达了拐角点.在此之后,波反射得到另外一个圆形波.通过与文献[11]比较,发现本文的数值结果与其数值结果是吻合的.

图3 (a) 算例 3 的网格剖分和数值解 w h 在不同时刻 (b) t=0.1,(c) t=0.3,(d) t=0.45 时的波传播Fig.3.(a) The triangulation mesh of Example 3;contour plot of solution w h at different time:(b) t=0.1 ;(c) t=0.3 ;(d) t=0.45.

5 结 论

研究了二阶波动方程在二维计算区域的传播问题,通过选取广义交替数值通量,建立了LDG方法并分析了格式的能量守恒性.时间离散上利用指数积分因子方法.数值实验验证了能量守恒性质和最优误差收敛阶.下一步的工作考虑将守恒的LDG格式应用到求解非线性波动方程,例如Sine-Gordon等方程上.