Chebyshev谱元法模拟波动的两种集中质量矩阵

Two Kinds of Lumped Mass Matrixes of Simulation of Wave Problem with Chebyshev Spectral Element Method

摘要: 波动模拟是工程结构抗震分析的重要手段之一。Chebyshev谱元法是一种高精度的微分方程求解方法，模拟波动问题时具有高精度和高效率的特点，因此受到广泛关注。Chebyshev谱元法求解波动问题时的质量矩阵是一致质量矩阵，即空间耦合矩阵，当和时域分析方法联用时，为时空耦合格式，每步计算需要联立求解线性方程组，计算效率受到了制约。本文在Chebyshev谱元法一致质量矩阵的基础上导出了两种集中质量矩阵：数学集中质量矩阵和物理集中质量矩阵，给出了两者的数学表达，并采用这两种集中质量矩阵下的Chebyshev谱元法结合时域中心差分法求解一维波动问题，此种模拟方案为时空解耦方法。数值分析表明，采用Chebyshev集中质量矩阵配合时域中心差分法模拟波动的方案具有较高的计算精度，并且每步计算不需要联立求解线性方程组，可以大幅度提高计算效率。其中，数学集中质量矩阵的计算精度要高于物理集中质量矩阵。

Abstract: Wave simulation is an important procedure for seismic analysis of engineering structure. Chebyshev Spectral Element Method is a kind of method to solve differential equations with high accuracy. Chebyshev Spectral Element Method has properties of high accuracy and high efficiency when it is used to simulate wave problem. The mass matrix of Chebyshev Spectral Element Method solving wave problem is consistent mass matrix which is space coupled. When it is used together with normal analysis procedure in time domain, it is time and space coupled. In each time step, it is needed to solve linear equations, and the efficiency is limited. In this paper, two kinds of lumped mass matrixes are derived based on Chebyshev Spectral Element Method consistent mass matrix. The math equations of the two mass matrixes are given. A procedure is given to solve one dimension wave problem. Chebyshev Spectral Element Method with the two matrixes is used in space domain, and central differential method is used in time domain. This procedure is decoupled in time and space domain. Numerical analysis shows that, this procedure has high accuracy, and it is not needed to solve linear equations in each time step, so compute efficiency can be largely enhanced. Within the two mass matrixes, math lumped mass matrix has higher accuracy than physical lumped mass matrix.

1. 前言

工程波动模拟是工程结构抗震分析的重要课题。谱元法是谱方法 [1] 和有限元方法结合的一种高精度微分方程求解方法。上世纪80年代，Patera提出了谱元法 [2]，其具有谱方法高精度的特点，同时又兼具有限元法的分块离散、便于求解复杂介质和复杂边界问题的特点。因此，谱元法在流体力学、计算传热学等领域的复杂问题求解中得到了广泛应用 [3] [4]。Seriani等首先将谱元法用于弹性波方程的求解 [5]，为工程波动模拟拓展了一条新的路径。此后众多研究者对谱元法模拟工程波动进行了深入的研究并取得了长足的进步 [6] [7] [8] [9]。谱元法根据基函数的不同主要分为两种：Legendre谱元法和Chebyshev谱元法。其中Chebyshev谱元法采用Chebyshev多项式作为基函数，积分方案可由公式显式表达，不需要采用数值积分，提高了积分精度并降低了积分过程的计算量，因此得到了广泛应用 [10]。Chebyshev谱元法通过变分原理导出的质量矩阵是一致质量矩阵，矩阵每个位置的元素都不为0，是空间耦合的矩阵。因此在模拟中，当结合传统的时域求解方法时，无论是显式方法还是隐式方法，计算过程都是时空耦合的。每一时步计算都需要联立求解线性方程组，这极大地增加了计算量。同时，在一致质量矩阵中，当某质点产生加速度时，矩阵的所有质点都会瞬间获得加速度，这和波动传播的物理概念并不相符。实际波动传播中，波速是有限的，一个质点获得的加速度在瞬时只能影响相邻质点，而不会使得远处的质点同时获得加速度。相较于一致质量矩阵，集中质量矩阵非对角线位置的元素皆为0，从物理概念上更加符合波速有限的基本原理；同时由于其空间解耦的特性，当结合时域显式差分方法时，计算过程不需要联立求解线性方程组，大幅度地降低了计算量。因此，在Chebyshev谱元法的基础上，提出相应的集中质量矩阵，并结合时域显式差分法，形成时空解耦的模拟方案，便成为了非常实际的考虑。

本文以一维波动为例，在Chebyshev谱元法的基础上导出了两种集中质量矩阵：数学集中质量矩阵和物理集中质量矩阵，给出了相应的数学表达。采用这两种集中质量矩阵，结合时域中心差分法，提出了时空解耦的Chebyshev谱元法模拟方案。数值实验结果表明，此种模拟方案具有很高的计算精度，并且在每步计算中不需要联立求解线性方程组，从而大幅度地提高了计算效率。

2. 波动模拟的Chebyshev谱元模型

为便于说明问题，在这里以一维波动问题为例。考虑如下的一维波动方程：

$\rho \frac{{\partial }^{2}u}{\partial t^2}=\mu \frac{{\partial }^{2}u}{\partial x^2}$ (1)

其中 $\rho$ 为传播介质密度， $\mu$ 为传播介质的抗剪模量。

对于式(1)，导出其等效弱积分形式：

${\displaystyle \underset{\Omega }{\int }\rho v\ddot{u}\text{d}x}+{\displaystyle \underset{\Omega }{\int }\nabla v\mu \nabla u\text{d}x}=0$ (2)

将一维空间区域划分为若干空间单元，单元的大小设定为 $\Delta x$。 $x_{j+1}-x_j=\Delta x$，单元内部的等效弱积分形式可以写为：

${\displaystyle \underset{\Omega }{\int }{\rho }^{e}v{\ddot{u}}^e\text{d}x}+{\displaystyle \underset{\Omega }{\int }\nabla v{\mu }^e\nabla u^e\text{d}x}=0$ (3)

将目标函数 $u^e$ 用一组空间正交基向量展开

$u^e\left(x\right)={\displaystyle \underset{i=0}{\overset{N-1}{\sum }}{a}_i{A}_i\left(x\right)}$ (4)

其中， $A\left(x\right)$ 为Chebyshev多项式。

在单元内部设立一组空间域控制点 $x_i$， $i=0,1,\cdots ,N-1$，取Chebyshev-Gauss-Lobatto的节点形式，具体表达式如下：

${\xi }_i=\cos \left(\frac{\pi i}{m}\right),i=0,1,\cdots ,m$ (5)

对于式(4)可以变换为如下以空间控制点的目标函数值表达的形式：

$u^e\left(x\right)={\displaystyle \underset{i=0}{\overset{N-1}{\sum }}{u}^e\left(x_i\right)l_i\left(x\right)}$ (6)

将式(6)代入式(3)并选取检验函数 $v=T\left(x\right)$， $T\left(x\right)$ 取Chebyshev多项式，可以得到单元内的半离散格式：

$M^e{\ddot{u}}^e\left(t\right)+K^e{u}^e\left(t\right)=0$ (7)

其中

$M_{ij}^e=\rho J\frac{4}{N^2{\bar{c}}_i{\bar{c}}_j}{\displaystyle \underset{k,l=0}{\overset{N}{\sum }}\left(\frac{1}{{\bar{c}}_k{\bar{c}}_l}{T}_k\left({\epsilon }_i\right)T_l\left({\epsilon }_j\right){\displaystyle \underset{-1}{\overset{1}{\int }}{T}_k}\left(\epsilon \right)T_l\left(\epsilon \right)\text{d}\epsilon \right)}$ (8)

$K_{ij}^e=J{\left(J^{-1}\right)}^2\frac{4}{N^2{\bar{c}}_i{\bar{c}}_j}{\displaystyle \underset{k,l=0}{\overset{N}{\sum }}\left(\frac{1}{{\bar{c}}_k{\bar{c}}_l}{T}_k\left({\epsilon }_i\right)T_l\left({\epsilon }_j\right){\displaystyle \underset{-1}{\overset{1}{\int }}{T^{\prime }}_k}\left(\epsilon \right){T^{\prime }}_l\left(\epsilon \right)\text{d}\epsilon \right)}$ (9)

其中J为空间单元变换至标准单元[0,1]的雅克比矩阵， $J^{-1}$ 为雅克比矩阵的逆矩阵。对于一维空间单元

$J=\frac{\Delta x}{2}$， $J^{-1}=\frac{2}{\Delta x}$。

将单元质量与刚度矩阵按照空间对应位置集成为整体质量和整体刚度矩阵，则得到半离散的整体求解方程：

$M\ddot{u}\left(t\right)+Ku\left(t\right)=0$ (10)

求解半离散方程(10)便可以得到整个波动模拟的全过程时空解答。

3. 集中质量矩阵模型

对于通常的 Chebyshev谱元法，通过变分原理导出的质量矩阵是一致质量矩阵，从矩阵形式上看是满阵，矩阵的每一个位置上的元素都不为零，这在物理上为空间耦连。在这里，分别导出两种适用于Chebyshev谱元法的集中质量矩阵。

3.1. 数学集中质量矩阵

首先给出网格参数取值 $\Delta x=1$， $N=\text{5}$ 时的一个Chebyshev一致质量矩阵，如下：

$M=\left[\begin{array}{ccccc}0.0444& -0.0071& -0.0158& 0.0198& -0.0079\\ -0.0071& 0.2539& 0.0507& -0.0507& 0.0198\\ -0.0158& 0.0507& 0.3301& 0.0507& -0.0158\\ 0.0198& -0.0507& 0.0507& 0.2539& -0.0071\\ -0.0079& 0.0198& -0.0158& -0.0071& 0.0444\end{array}\right]$

对所列的Chebyshev一致质量矩阵进行观察，远离对角线的位置上的元素数值相较于对角线位置上的元素数值是比较小的，这在物理上是可以得到解释的。因为质点之间越是互相远离，其相互产生扰动的程度会越小。根据一致质量矩阵的此种特性，从数学角度提出一种集中质量矩阵。

将一致质量矩阵的每一行的元素全部相加，放置在这一行与矩阵对角线相交的位置上，这一行的其他元素置0。公式表达如下：

一致质量矩阵表达为M，数学集中质量矩阵表达为 $M_s$，M和 $M_s$ 的大小均为 $N\times N$。

$M_s\left(i,i\right)={\displaystyle \underset{j=1}{\overset{N}{\sum }}M\left(i,j\right)},i=1,2,\cdots ,N;j=1,2,\cdots ,N$ (11.a)

$M_s\left(i,j\right)=0,i=1,2,\cdots ,N;j=1,2,\cdots ,N;i\ne j$ (11.b)

3.2. 物理集中质量矩阵

廖振鹏 [11] 及其研究组针对有限元法模拟波动问题，提出了有限元法的集中质量矩阵。从物理概念出发，将有限元网格点所占据的质量代表值赋予质量矩阵对角线的相应位置，并将此集中质量矩阵和中心差分法相结合提出了时空解耦的波动模拟方法，取得了非常好的效果。在此给出Chebyshev谱元模拟的物理集中质量矩阵：空间节点和其相邻节点平均分配他们的共有质量，非对角线位置上的元素皆置0。公式表达如下：

单元空间布置点的坐标为 $x_i$， $i=1,2,\cdots ,N$，物理质量矩阵表达为 $M_w$

$M_w\left(i,i\right)=\frac{x_{i+1}-x_{i-1}}{2}\rho ,i=2,\cdots ,N-1$ (12.a)

$M_w\left(i,i\right)=\frac{x_2-x_1}{2}\rho ,i=1$ (12.b)

$M_w\left(i,i\right)=\frac{x_N-x_{N-1}}{2}\rho ,i=N$ (12.c)

$M_w\left(i,j\right)=0,i=1,2,\cdots ,N;j=1,2,\cdots ,N;i\ne j$ (12.d)

将式(10)中的M替换为式(11)和式(12)中的 $M_s$ 和 $M_w$，则式(10)可以改写为：

$M_s\ddot{u}\left(t\right)+Ku\left(t\right)=0$ (13)

或者

$M_w\ddot{u}\left(t\right)+Ku\left(t\right)=0$ (14)

式(13)和式(14)即为经过Chebyshev谱元离散后，两种集中质量矩阵下对应的半离散方程，求解式(13)和式(14)即可得到波动的时空解答。

4. 时空解耦模拟方案

对于式(13)和式(14)，为了形成时空解耦的计算格式，采用中心差分法进行求解。以式(13)为例，其具有如下的过程：

将式(13)写成离散表达式

$M_s{\ddot{u}}_i\left(t\right)+K{u}_i\left(t\right)=0$ (15)

式中 ${\ddot{u}}_i,u_i$ 分别代表 $t_i$ 时刻的加速度反应和位移反应。

中心差分法基于有限差分代替位移对时间的求导，假定时间域[0,T]内采用等时间步长， $\Delta t_i=\Delta t$，则加速度的中心差分近似表达为

${\ddot{u}}_i=\frac{u_{i+1}-2u_i+u_{i-1}}{\Delta t^2}$ (16)

将式(16)代入式(15)中去，可以得到

$M_s\frac{u_{i+1}-2u_i+u_{i-1}}{\Delta t^2}+K{u}_i=0$ (17)

由式(17)可得

$\frac{1}{\Delta t^2}{M}_s{u}_{i+1}=-\left(K{u}_i-\frac{2}{\Delta t^2}{M}_s\right)u_i-\frac{2}{\Delta t^2}{M}_s{u}_{i-1}$ (18)

式(18)即为半离散方程(15)的中心差分求解式。

从式(18)可以看出，i + 1时刻点的反应由i时刻和i − 1时刻的反应决定， $u_{i+1}$ 可以由 $u_i$ 和 $u_{i-1}$ 显式表达，计算中不需要联立求解线性方程组，计算过程是时空解耦的，和一致质量矩阵时空耦合的计算格式相比，计算效率上有了大幅度的提高。

5. 数值实验

一维空间计算区域取100 m，在坐标50 m处，施加一个初始位移函数，位移函数为三次B样条函数，具体函数表达式见式(19)、(20)。

$P\left(\tau \right)=16\left[G\left(\tau \right)-4G\left(\tau -\frac{1}{4}\right)+6G\left(\tau -\frac{1}{2}\right)-4G\left(\tau -\frac{3}{4}\right)+G\left(\tau -1\right)\right]$ (19)

$G\left(\tau \right)={\tau }^{3}H\left(\tau \right)$ (20)

式中 $\tau =\frac{x-50}{X}$，X代表初始位移函数施加区域的大小。

初始位移函数的施加区域为(50 m, 70 m)，即 $X=20$。此问题是一个波动的柯西初值问题。初始速度取为0，波速取为1 m/s。

采用本文提出的两种集中质量矩阵下Chebyshev谱元法进行空间离散，而后在时间域采用中心差分法进行求解。把计算结果和解析解进行对比分析。在本例中，谱元法的计算参数是：空间单元大小 $\Delta x=2\text{m}$， $N=5$。为有效地考察两种集中质量矩阵的计算效果，需排除时间域计算精度对计算结果的影响，因此将中心差分法的时步取为0.02 s，可以确保时间域的计算精度能够处在一个较高的水准。

图1和图2给出了在给定网格下的两种集中质量矩阵的计算效果和解析解的对比。图1列出了10 s，20 s，30 s三个时刻的空间波形。图2列出了50 m，60 m，70 m三个空间位置点上的位移时程反应。从图中可以看出，两种集中质量矩阵都很好地模拟了波动的传播过程，无论是空间波形还是时程都与解析解吻合得非常好，计算误差在图示比例尺中可以忽略。

Figure 1. Spatial waveform comparison

图1. 空间波形对比

Figure 2. Comparison of displacement time history

图2. 位移时程对比

图3给出了两种集中质量矩阵和解析解关于空间波形和时程的峰值的相对误差的对比。首先提取出10 s，20 s，30 s和40 s的空间波形，然后分别取这四条波形的峰值位移。将两种集中致量矩阵的计算结果和解析解进行对比，结果见图3(a)，纵坐标表示两种集中质量矩阵计算结果与解析解的相对误差。同样的方式，提取出50 m，60 m，70 m，80 m处的波动时程，然后取四条时程的峰值位移，再将两种集中质量矩阵结果和解析解做比较，对比结果见图3(b)。从图中可以看出，两种集中质量矩阵的计算相对误差都处于较低的水平。在算例所取的网格参数下，最大的相对误差值为0.279%，最小的相对误差值为0.0257%。在两种集中质量矩阵之间进行对比，数学集中质量矩阵的计算相对误差值大约为物理集中质量矩阵的1/5~1/3。

Figure 3. Comparison of peak relative error

图3. 峰值相对误差对比

6. 结论

在传统Chebyshev谱元法模拟波动的一致质量矩阵的基础上，提出了两种集中质量矩阵。以一维波动为例进行了阐述，给出了具体数学表达。将采用两种集中质量矩阵的Chebyshev谱元法和时域中心差分法相结合形成了时空解耦的波动模拟方案。进行了数值实验，得到结论如下：

1) 两种集中质量矩阵都可以很好地模拟波动的传播过程，且具有较高的计算精度。

2) 数学集中质量矩阵的计算精度高于物理集中质量矩阵。在本文给定参数下，数学集中质量矩阵的相对误差约为物理集中质量矩阵的1/5~1/3。

3) 采用本文提出的两种集中质量矩阵的谱元法结合时间域中心差分法可以使得计算过程成为一个时空解耦的过程，不需每时步联立求解线性方程组，可以大幅度提高计算效率。