固端法：二维有限元先验定量误差估计与控制

袁 驷，袁 全

(清华大学土木工程系，北京 100084)

文献[1]通过对矩阵位移法[2]和有限元法[3−4]的分析比较，得出一个结论：一维有限元的误差主要来源于各个单元的“固端解”项。其后，基于恢复单元固端解这一思想，超收敛计算的单元能量投影(Element Energy Projection，简称EEP)法得以创立并取得长足发展，不仅对一维有限元法 (Finite Element Method, 简称 FEM)[5−8]，对二维有限元线法(Finite Element Method of Lines, 简称FEMOL、线法)[9]、二维乃至三维有限元法[10−11]都建立了EEP 超收敛算法，也得到了数学理论上的证明[12−13]。更有意义的是，用EEP 超收敛解替代精确解来估计常规有限元解的误差，使得基于EEP 技术的自适应有限元求解得以实现，其最突出的特点是可以得到按最大模逐点满足用户给定的误差限的解答，可谓是数值精确解[14−15]。目前，这种自适应有限元方法不仅有效地应用于各种线性问题，也在特征值问题和多种非线性问题中得到了广泛而有效的应用[16−18]，而近期发展的网格局部加密技术为一类刁难奇异问题的自适应求解提供了更高性能的求解方案[19−20]。

纵观各类自适应求解，几乎都有一个共同(弱)点：因为解答事先未知，只能用后验误差方法，按照有限元求解、误差估计、更新网格三步循环迭代求解。这里的关键问题是：缺少先验定量的误差估计。这是因为，目前几乎所有的先验误差估计，都包含了事先不可计算的因素在其中，难以定量，只能是定性的。

本文作者经过对文献[1]的反思和进一步研究，在文献[21]中对于一维Ritz 有限元法提出了先验定量误差估计的“固端法”，从而可以不经有限元计算，便可以一举给出满足精度要求的网格划分。继“固端法”在一维有限元中的初步成功，本文尝试将其拓展到二维有限元。作为初始探索，本文以Poisson 方程为例，采用最常规的4 结点线性元，用EEP 技术预先做出误差估计，能直接给出满足精度要求的网格划分。本文对这一最新进展做一简要介绍，并给出初步的数值结果。

1 一维问题简述

一维有限元的误差主要来自于其丢失的单元“固端解”项[1]，可以分两种情况论述。

1) 精确单元：其形函数是齐次控制微分方程的通解，结点位移是精确的，固端解不折不扣地就是精确单元内部的误差。常截面杆件矩阵位移法的单元即是精确单元，参见文献[1]中的图1：其状态(Ⅱ)为有限元解，而状态(Ⅰ)为固端解，亦是有限元解的误差。

2) 近似单元：其形函数不是齐次控制微分方程的通解，结点位移是近似的，单元内部误差由固端解和非固端的有限元解共同组成。但是，有限元的数学理论已有证明，有限元的端结点位移相比于单元内部位移是超收敛的，而且具有最佳超收敛性[4]；以C0类单元为例(文献[1]的表2)， m次单元在单元内部为 O(hm+1)收敛，而在结点上是O(h2m)收敛的。所以，对于近似单元，特别是高次单元，相对于单元内部位移的误差，结点位移的误差是高阶微量，可以合理地将其略去，亦即近似单元误差的主要来源亦为固端解。

这样，精确单元和近似单元的误差的主要来源便得到了统一，即单元的固端解项。然而，更大的利好是，求固端解是局部单元的问题，并不需要作整体的有限元求解。这就使得不经有限元整体求解而预先对有限元解的误差做出定量估计成为可能。

对于二维问题，半离散的有限元线法可以看作是广义一维有限元法(结点延伸为结线)，其误差的主要来源为线法单元两端结线(单元边结线)固定的解，所以仍可以说是来源于“固端解”。而二维有限元法可以看作用一维有限元求解线法的结线位移，而当结线位移事先固定为零时，线法和有限元法就没有区别了，在有限元网格上也可以定义固端解了。

2 模型问题和有限元解

2.1 Poisson 方程

本文的分析以二维Poisson 方程为模型问题，其形式如下：

2.2 二维有限元解

用二维有限元求解此问题时，本文采用双向m 次四边形单元，单元试探函数为：

2.3 拟线法网格

求解域 Ω不必是规则区域，但是为适合二维EEP 超收敛计算，网格划分须是“拟线法网格”，即先利用FEMOL 的离散方式用一组结线对求解区域进行半离散，然后再沿结线维度进一步离散，得到的网格即是求解FEMOL 常微分方程组(Ordinary Differential Equations，简称 ODEs)的一维有限元网格，也是原问题的二维有限元网格，如图1所示。

图 1 二维问题的逐维离散Fig. 1 Dimension-by-dimension discretization of 2D problems

2.4 有限元线法

以下简单介绍FEMOL[22]的概念。图2 为一个典型FEMOL 二次单元的几何映射。

图 2 FEMOL 单元映射Fig. 2 Mapping of FEMOL element

2.5 线法的 EEP 解

3 固端法

3.1 固端解

3.2 线性元公式

3.3 先验误差控制

4 数值算例

表 1 例1 弹性扭转问题的结果Table 1 Results of elastic torsion problem in Example 1

图 3 方形区域弹性扭转问题的误差Fig. 3 Error of elastic torsion problem on a square region

例2. 非规则四边形求解域问题

本例考虑一个Poisson 方程定义在如图4 所示的非规则四边形求解域上，四周为齐次本质边界条件，问题描述如下：

f 由u 反求得到，形式复杂，在此略去。

本例的区域为非规则区域，荷载也是一个函数，是一个很典型的例题。本例将采用两种方法估算误差：一是基于式(8)或式(10)的“直估法”；二是基于式(6)的“先验法”。

图 4 非规则四边形求解域Fig. 4 An irregular quadrilateral solution domain

为了实施“直估法”，即用式(8)或式(10)直接估算单元大小，需对本例荷载和区域作一些近似处理：本例的荷载 f是函数，大约在区域中心，可找到其最大值 f0≃0.9；本例的区域非矩形，将其近似为边长为5 的方形区域，以便估算沿边长所用的单元数。以上处理，相当于在 5×5的区域上有均布荷载 f0作用的Poisson 方程。这样就可以用式(8)或式(10)直接估算单元大小(模仿例1)，并直接确定网格。“先验法”则用式(6)逐单元估算固端解误差，需投入少量的计算，但比“直估法”更加精准一些。

问题(Ⅰ)求解结果列于表2。其中直估法的计算极为简单，例如对于2×2 网格(参见图5)，由式(8)可直估误差为 f0h2/8=0.9·(5/2)2/8=0.7031。可以看出，简简单单的直估法虽然在网格比较稀疏时比先验法估计得要偏高一些，但随着网格加密，二者完全趋于一致。这是因为，本例的最大误差发生在区域中心，而当网格很密、单元很小时，区域中心单元上的荷载趋于常数 f0，因此与直估法趋于一致。这也辅证了直估法的合理性。还可以看到，无论是直估法还是先验法，所估算的误差都略大于真实误差；作为误差估计器，是偏于安全和可靠的。

表 2 例2 非规则区域问题(I)的结果Table 2 Case (I) results of irregular domain problem in Example 2

图 5 二维有限元网格Fig. 5 A 2D FEM mesh

问题(Ⅱ)求解结果列于表3。最后一列给出最大误差比，可以看出，其值均小于1，亦即预估网格的结果都满足误差限。还可看出，真实误差均略小于预估误差，安全而可靠。

表 3 例2 非规则区域问题(II)的结果Table 3 Case (II) results of irregular domain problem in Example 2

5 结论

本文提出二维有限元先验定量误差估计的“固端法”，可以根据给定的误差限，直接确定允许的单元大小，一举得到允许的网格划分，极大地简化了二维有限元误差估计的计算，大量减少自适应有限元求解的迭代步骤，并大幅提升自适应有限元求解的效率。该法的更为深入、广泛、系统的研究成果将另文介绍。