谱元法应用分析

谱元法应用分析
汪大伟;李南生
【摘要】首先推导了谱元法在Helmholtz方程和热传导方程中的计算方法,然后以一个规则区域上的Laplace方程和路堤温度场分布的简化模型为算例,对比了有限元法和谱元法的计算效果,显示了谱元法在解决特定实际工程问题中的优越性.【期刊名称】《佳木斯大学学报（自然科学版）》
【年(卷),期】2016(034)003
【总页数】3页(P338-340)
【关键词】有限元法;谱单元;h-型;p-型;高精度;谱收敛性
【作者】汪大伟;李南生
【作者单位】同济大学土木工程学院,上海200092;同济大学土木工程学院,上海200092
【正文语种】中文
【中图分类】TV312
1 Helmholtz方程
Helmholtz方程广泛的应用价值。

该方程表示如下：
2u(x1,x2 )-λ2u(x1,x2 )=f(x1,x2)
其中∂Ω表示区域边界。

由于Helmholtz方程满足自耦性，由这一性质立即可以得到其等价形式，即该问题等价于求其变分形式的极值[1]，但这里直接由该方程的弱形式出发，即：找到
u∈S，使得∀w∈V，有
-∫Ωw·udΩ-λ2∫ΩwudΩ=∫ΩwfdΩ
其中|u(x1,x2 )=0(x1,x2)∈∂Ω}
由于考虑齐次Dirichlet边界条件，所以S,V属于相同的函数空间。

传统的有限元算法，在每个单元上采用低阶多项式来构造分段多项式空间，通过细化网格，缩小单元尺寸，增加基函数个数，扩大搜索空间的方法，来尽可能地逼近原函数，减小逼近误差，而不考虑基函数的类型，显然这种方法的效率是比较低的(当然该方法在某些方面，如复杂边界条件下，有其灵活性)。

而谱(元)方法的有效性在于，其选择的基函数为奇异Sturm-liouville问题的特征函数，所谓奇异Sturm-Liouville问题，即满足一定边界条件的下述微分方程：
(p(x)y′)′+(q(x)+λr(x) )y=0,x∈(a,b)
对于多维的情况，我们的插值函数取相应一维插值函数的张量积形式。

具体构造如下：对于几何坐标插值，有
其中x=(x1,x2 ), 为物理坐标，ξ=(ξ1,ξ2)为为x映射到参数坐标上的对应坐标值，xA和NA分别为节点A出的坐标值和坐标插值函数，而坐标插值函数按一维张量积的形式构造为
其中δAB表示Kronecker delta符号。

类似地，有对未知函数的插值公式
其中uA为节点A处对应的待求函数值，其他定义同上。

注意，这里虽然用了相
同的符号NA和A来分别表示坐标及未知函数的插值函数和节点个数，但是两者可以不相同。

对于试函数w，上面已经说明，采用与u相同的函数空间
下面进行单元节点的布置，对于多维情况，本文采用Gauss-Lobatto-Legendre 节点的Lagrange张量积的形式，在每个参数坐标轴方向，其节点由下式给出：[2] LN′(rk )=0 k=1,2,…,N-1
rN=1,
其中r0～rN表示沿着参数坐标轴方向布置的节点，LN为N阶Legendre多项式，LN′为其导数，即N-1阶Lobatto多项式，上式表示沿每个参数坐标方向，N+1
个节点包括2个端点，以及中间的N-1个节点取N-1阶Lobatto多项式的零点。

为了表示方便，将上面得到的弱形式(3)写成下列形式：
-∫Ωw,iu,idΩ-λ2∫ΩwudΩ=∫ΩwfdΩ
其中i=1,2，w,i和u,i分别表示对u和w的第i个坐标分量的偏导数，这里采用
了爱因斯坦求和约定，以下采用相同的表示方法。

对该式按单元求和
将该积分式转化为参数坐标表示，在每个单元上我们代入式(6)，(7)，(8)并化简可以得到最终需要求解的线性方程组：
其中Λ表示单元集成，单元矩阵分别为
其中I，J=1,2，且
det[J]表示Jacobian行列式，为
而对于上述积分的求和，需要采用Gauss-Lobatto数值积分方法，即选取与前面(11)式给出的函数插值点相同的单元节点作为Gauss积分点，最后，对单元集成，得到Helmholtz问题的总体线性方程组
到此实现了Helmholtz问题的谱元法求解。

1.2 热传导方程
热传导方程是典型的二阶抛物线形偏微分方程，其具体形式为
其中系数c和k对于在热传导问题中分别表示比热容和热传导系数，同时该方程
满足一定的边界条件和初始条件，同上节所述的Helmholtz方程相比，热传导方
程主要的不同在与其未知函数是随时间变化的，将函数按时间和空间分离进行展开，展开式如下：
(22)
待求未知量uA(t)随时间变化。

这样，按照与上节类似地方法，可得该问题归结为
求解下列一阶常微分方程组
其中表示未知量对时间求导。

对于这一问题可以采用许多不同的数值解法，如Euler法、Adams法、Runge-Kutta法，本文的算例将采用Crank-Nicoleson方法，它属于Euler方法的一种，将(23)转化为下列线性方程组的求解
其中上标n和n+1分别表示第n和n+1时间步时的函数节点函数值。

2.1 Laplace方程
令Helmholtz方程的系数和右边项为零，得到经典的Laplace方程
其中(x,y)∈Ω，为了验证解得精度，采用虚构解方法(Method of Manufactured Solutions)来制造一个精确解[8]，令方程的精确解为
u为无限光滑解，为了简化计算我们取求解域为Ω=[0,1]×[0,1]，并且取齐次Dirichlet边界条件，使满足该精确解，具体如下：
u=sin1·e-y, x=1
u=sinx, y=0
u=e-1sinx, y=1
计算结果列于表1和表2，其中在有限元法中本文采用了常用的二阶单元。

分析表中结果，对于有限元法，我们看到单元数每增加4倍，误差则将为1/8左右，这与前面所述理论分析的结果是一致的，即‖。

而谱元法的收敛速度明显高于传统的有限元法，其误差按指数级减小，达到了所谓的谱收敛性。

2.2 热传导方程
算例为路堤在受到周期性变化的环境温度影响下，其内部的温度场分布情况。

设该区域土的比热容为C=1608J·kg-1·℃-1，热传导系数k=1.42W·m-1·℃-1[3],路堤温度场计算时下界面EH取为一定深度下的水平面，温度边界值为恒值
TG=10.26℃。

路基为南北走向，不考虑其阴阳坡效应，忽略路堤边坡和路肩、路面表面的温度差异，统一取为地表面气温变化的余弦时程函数：
Tu (t)=Tf+Im{Aexp[i(ωt+φ) } (29)
其中Tf=22℃, A=20℃, φ=2.1981, 初始温度为T(x,y,0)=10.26℃，周期为365d，即w=6.342×10-8 π。

谱元法在一般的实际应用中具有精度高，收敛速率快，计算高效的特点，相较于传统的有限元法优势明显，因而有很大的应用价值。

但是，同时也应当看到该方法也有其不足之处。

首先对于几何形状复杂的问题，不如传统的有限元法灵活性高，还有难点有待解决；同时，对于大规模的计算问题，由于该方法所得到的线性方程组的系数矩阵一般都为稠密矩阵，所以给计算带来了很大的困难，常用的直接法和迭代法一般只适用于稀疏矩阵方程组的求解，因而需要找到更加高效的算法来解决这一问题，例如如果原函数是按Fourier展开，使用快速傅里叶变换(FFT)则可以大
大提高计算效率。

【相关文献】
[1] S.C. Mikhlin, Variational Methods in Mathematical Physics[M], Macmillan,New York, 1964.
[2] A.T. Patera. A Spectral Element Method for Fluid Dynamics: Laminar Flow in Channel Expansion[J]. Journal of computational physics,1984,54:468-488.
[3] Li N, Tang B. An Approximate Analytic Method for Calculating the Temperature Field
for a Phase Transition in an Embankment in a Permafrost Region[J]. Heat and Mass Transfer, 2015, 51(8): 1097-1109.。