[有限元方法基础理论] 质量集中有限元(谱元)

[有限元方法基础理论] 质量集中有限元(谱元)

集中质量有限元,也称谱元,要求具有对角质量矩阵,且对角线上的值等于行和。可以通过用数值积分公式以及其他方法来构造谱元。

不妨考虑一个初边值问题,
u t − Δ u = f  in  Ω , t > 0 u = 0  on  ∂ Ω , t > 0 ,  with  u ( ⋅ , 0 ) = v  in  Ω \begin{aligned} &u_{t}-\Delta u=f \quad \text { in } \Omega, \quad t>0\\ &u=0 \quad \text { on } \partial \Omega, t>0, \quad \text { with } u(\cdot, 0)=v \quad \text { in } \Omega \end{aligned} utΔu=f in Ω,t>0u=0 on Ω,t>0, with u(,0)=v in Ω
它的半离散方法可以写为:寻找 u h : [ 0 , ∞ ) → S h u_{h}:[0, \infty) \rightarrow S_h uh:[0,)Sh,使得
( u h , t , χ ) + ( ∇ u h , ∇ χ ) = ( f , χ ) , ∀ χ ∈ S h , t > 0 , u h ( 0 ) = v h \left(u_{h, t}, \chi\right)+\left(\nabla u_{h}, \nabla \chi\right)=(f, \chi), \quad \forall \chi \in S_{h}, t>0, \quad u_{h}(0)=v_{h} (uh,t,χ)+(uh,χ)=(f,χ),χSh,t>0,uh(0)=vh
写成矩阵的形式,就是
B α ′ ( t ) + A α ( t ) = f ~ ( t ) ,  for  t > 0 ,  with  α ( 0 ) = γ \mathcal{B} \alpha^{\prime}(t)+\mathcal{A} \alpha(t)=\widetilde{f}(t), \quad \text { for } t>0, \quad \text { with } \alpha(0)=\gamma Bα(t)+Aα(t)=f (t), for t>0, with α(0)=γ
这里的 B \mathcal{B} B是质量矩阵, A \mathcal{A} A是刚度矩阵。以上就是标准的有限元离散到代数方程组,用到的一些符号的具体的含义不言自明,就不细说了。

所谓的质量集中有限元,我们就是希望,这里得到的质量矩阵 B \mathcal B B,满足条件:它是个**对角矩阵,对角线上的值刚好等于行和。**不妨用 B ‾ \overline{\mathcal{B}} B来表示。

那么问题来了,通过什么样的方式,上述的要求才能被满足呢?我们这里提两种方法,一个是通过数值积分的方法构造,另外一个是通过构造分片常数有限元空间的方法来构造。实现的时候,总是方便的,但是在收敛性分析的时候,我们希望有它的一种构造方式作为桥梁。

采用数值积分公式计算单元质量矩阵时,当求积点与插值点一致时,就可以得到谱元。积分公式要求:达到指定阶,且权重全正。

对于分片 k k k 次多项式有限元,精确计算质量矩阵的积分 公式需安达到 2 k 2 k 2k 阶。对常见二阶 PDE,积分公式只需要达到 2 k − 2 2 k-2 2k2 阶便可保 证数值解精度。

考虑最简单的线性元 P 1 P_{1} P1插值点为三个顶点 { v 0 , v 1 , v 2 } \left\{v_{0}, v_{1}, v_{2}\right\} {v0,v1,v2}, 将其作为求积点,积分公式可达到 1 阶且权重 ( 1 3 ) \left(\frac{1}{3}\right) (31) 为正,满足要求。下面我们来验证一下,它满不满足质量集中的条件。

在单元 τ \tau τ上,它的以顶点为积分点的数值积分公式为,
Q τ , h ( f ) = 1 3  area  ( τ ) ∑ j = 1 3 f ( P τ , j ) ≈ ∫ τ f d x Q_{\tau, h}(f)=\frac{1}{3} \text { area }(\tau) \sum_{j=1}^{3} f\left(P_{\tau, j}\right) \approx \int_{\tau} f d x Qτ,h(f)=31 area (τ)j=13f(Pτ,j)τfdx
刚好是三个顶点的平均值再乘以面积。为了方便,我们接下来用如下符号来表示数值积分意义下的内积,
( ψ , χ ) h = ∑ τ ∈ T h Q τ , h ( ψ χ ) (\psi, \chi)_{h}=\sum_{\tau \in \mathcal{T}_{h}} Q_{\tau, h}(\psi \chi) (ψ,χ)h=τThQτ,h(ψχ)
下面只要证明,对于
( u h , t , χ ) h + ( ∇ u h , ∇ χ ) = ( f , χ ) , ∀ χ ∈ S h , t > 0 , u h ( 0 ) = v h \left(u_{h, t}, \chi\right)_{h}+\left(\nabla u_{h}, \nabla \chi\right)=(f, \chi), \quad \forall \chi \in S_{h}, t>0, \quad u_{h}(0)=v_{h} (uh,t,χ)h+(uh,χ)=(f,χ),χSh,t>0,uh(0)=vh
代数化得到的质量矩阵,满足集中质量的要求即可。事实上,若令 u h ( t ) = ∑ j = 1 N h α j ( t ) Φ j u_{h}(t)=\sum_{j=1}^{N_{h}} \alpha_{j}(t) \Phi_{j} uh(t)=j=1Nhαj(t)Φj,则代数方程组可以写为:
∑ j = 1 N h α j ′ ( t ) ( Φ j , Φ k ) h + ∑ j = 1 N h α j ( t ) ( ∇ Φ j , ∇ Φ k ) = ( f , Φ k ) , k = 1 , … , N h \sum_{j=1}^{N_{h}} \alpha_{j}^{\prime}(t)\left(\Phi_{j}, \Phi_{k}\right)_{h}+\sum_{j=1}^{N_{h}} \alpha_{j}(t)\left(\nabla \Phi_{j}, \nabla \Phi_{k}\right)=\left(f, \Phi_{k}\right), \quad k=1, \ldots, N_{h} j=1Nhαj(t)(Φj,Φk)h+j=1Nhαj(t)(Φj,Φk)=(f,Φk),k=1,,Nh
只要证明两点:

  • ( Φ j , Φ k ) h = 0 \left(\Phi_{j}, \Phi_{k}\right)_{h}=0 (Φj,Φk)h=0 for j ≠ k j \neq k j=k,既然是0狄利克雷边界。
  • ∥ Φ j ∥ h 2 = ( Φ j , Φ j ) h = ∑ k = 1 N h ( Φ j , Φ k ) \left\|\Phi_{j}\right\|_{h}^{2}=\left(\Phi_{j}, \Phi_{j}\right)_{h}=\sum_{k=1}^{N_{h}}\left(\Phi_{j}, \Phi_{k}\right) Φjh2=(Φj,Φj)h=k=1Nh(Φj,Φk)

第一点是自然的,因为我们可以看到,不同的形函数的乘积在顶点处的值都等于 0,那么数值积分值就等于 0。下面来说明第二点。

首先,需要说明的是,在同一个单元上,两个不同的形函数的之间的乘积积分等于 1/12 的三角形面积,而形函数自身和自身的乘积的积分等于 1/6 的单元面积,即
∫ τ Φ j Φ k d x = 1 12  area  ( τ )  和 ∫ τ Φ j 2 d x = 1 6  area  ( τ ) \int_{\tau} \Phi_{j} \Phi_{k} d x=\frac{1}{12} \text { area }(\tau) \text { 和} \int_{\tau} \Phi_{j}^{2} d x=\frac{1}{6} \text { area }(\tau) τΦjΦkdx=121 area (τ) τΦj2dx=61 area (τ)
这个可以通过标准单元变换和简单计算得到。可以看到,若用 D j D_j Dj表示和顶点 P j P_j Pj有关的三角形区域的并集,那么原来的那个质量矩阵的行和就表达为
∑ k = 1 N h ( Φ j , Φ k ) = 1 3  area  ( D j ) \sum_{k=1}^{N_{h}}\left(\Phi_{j}, \Phi_{k}\right)=\frac{1}{3} \text { area }\left(D_{j}\right) k=1Nh(Φj,Φk)=31 area (Dj)
这个不清楚的话,画一下图,就可以看明白。一个基函数在相关的某个三角形单元上,自身和自身作用了一次,和同一三角形上其他顶点对应的基函数各作用了一次。

我们再看利用数值积分公式计算内积得到的质量矩阵对角线上的值,
∥ Φ j ∥ h 2 = ∑ Q τ , h ( Φ j 2 ) = 1 3  area  ( D j ) \left\|\Phi_{j}\right\|_{h}^{2}=\sum Q_{\tau, h}\left(\Phi_{j}^{2}\right)=\frac{1}{3} \text { area }\left(D_{j}\right) Φjh2=Qτ,h(Φj2)=31 area (Dj)

这个刚好就等于原始质量矩阵的行和,这就说明了,用这样一个数值积分格式得到的质量矩阵刚好满足质量集中的要求。

二次元 P 2 P_{2} P2插值点为三个顶点 { v 0 , v 1 , v 2 } \left\{v_{0}, v_{1}, v_{2}\right\} {v0,v1,v2} 和三条边中点 { e 0 , e 1 , e 2 } \left\{e_{0}, e_{1}, e_{2}\right\} {e0,e1,e2} ,相应的积分公式可 以达到 2 阶精度。但是,三个顶点的权重为 0 0 0。一般采取增加自由度或者扩展多项式空间的方式构造谱元。
三次元 P 3 P_{3} P3可以找到由三个顶点、每条边上两个对称的内点以及三 角形重心构成的 4 阶积分公式。但是,三个顶点的权重是负的!

下面我们看用另外一种方法来构造谱元。

对于一个三角形单元 τ \tau τ,我们把每个顶点和它对边的中点用直线段连接,这些直线段就交于三角形的重心,并且把三角形分成了等面积的 6 块。我们用 B j , τ B_{j, \tau} Bj,τ来表示六块中和顶点 P j P_j Pj相邻的两块,它的面积刚好是三角形面积的1/3。我们用 B j B_j Bj来表示顶点 P j P_j Pj周围的 B j , τ B_{j, \tau} Bj,τ的并区域。这个自己画个图就可以很容易理解。

那么,我们就可以以 B j B_j Bj为分割单元,构造一个分片常数有限元空间 S ˉ h \bar{S}_{h} Sˉh,里面的元素可以表示为,
χ ˉ ( x ) = ∑ j = 1 N h χ ˉ ( P j ) Φ ˉ j ( x ) \bar{\chi}(x)=\sum_{j=1}^{N_{h}} \bar{\chi}\left(P_{j}\right) \bar{\Phi}_{j}(x) χˉ(x)=j=1Nhχˉ(Pj)Φˉj(x)
这里的 Φ ˉ j = 1 \bar{\Phi}_{j}=1 Φˉj=1是关于区域 B j B_j Bj的特征函数。

那么只针对原来的质量矩阵这一块进行处理,
( u ˉ h , t , χ ˉ ) + ( ∇ u h , ∇ χ ) = ( f , χ ) , ∀ χ ∈ S h , t > 0 , u h ( 0 ) = v h \left(\bar{u}_{h, t}, \bar{\chi}\right)+\left(\nabla u_{h}, \nabla \chi\right)=(f, \chi), \forall \chi \in S_{h}, t>0, \quad u_{h}(0)=v_{h} (uˉh,t,χˉ)+(uh,χ)=(f,χ),χSh,t>0,uh(0)=vh
很容易看到,这样离散得到的质量矩阵,是满足集中质量的要求的。事实上, ( Φ ˉ j , Φ ˉ k ) = 0 \left(\bar{\Phi}_{j}, \bar{\Phi}_{k}\right)=0 (Φˉj,Φˉk)=0 j ≠ k j \neq k j=k。并且 ∥ Φ ˉ j ∥ 2 = \left\|\bar{\Phi}_{j}\right\|^{2}= Φˉj2= area ( B j ) = \left(B_{j}\right)= (Bj)= area ( D j ) / 3 \left(D_{j}\right) / 3 (Dj)/3

关于其误差分析,可以利用方法一提到数值积分的构造,参考Thomée V. Galerkin finite element methods for parabolic problems[M]. Berlin: Springer-Verlag, 1984.。

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: Age of Ai 设计师:meimeiellie 返回首页