降阶单元的新进展：内置了最大模误差估计器的自适应有限元法初探

袁 驷，杨 帅，袁 全，王亦平

(清华大学土木工程系，土木工程安全与耐久教育部重点实验室，北京 100084)

自适应有限元分析可以提升求解质量、精度和效率，是近年来数值计算方法研究的热点[1-10]。自适应求解的前提是，能够对有限元解进行可靠的误差估计，以指导自适应网格划分。一维有限元的误差估计相对简单，对于二维和三维有限元，纵观当前大多数后验误差估计方法，都存在或部分存在如下缺点和弊端：① 不能在单个单元上估计误差，需要结构化网格，需要若干单元联片；② 不能逐点估计误差，难以按最大模控制误差；③ 线性元由于整体上缺失超收敛性，难以构造可靠的误差估计器；④ 误差估计缺乏理论证明，可靠性缺乏理论保障。

对于m(≥1)次常规多项式单元，有限元数学理论有一个最基本的误差估计，即若问题足够光滑，则一维单元到三维单元在单元内部任一点的位移无例外地具有O(hm+1)的收敛性[11-14]。这一估计意味着二分加密网格(h→h/2)或提高单元次数(m→m+1)都能获得精度更高的解，可用来估计原有限元解的误差。如此便很自然地派生出两种略显原始且笨拙的误差估计策略：一是用二分加密网格上的解估计原网格上解的误差，本文称为“双元法”；二是用高一次单元的解估计原单元解的误差，本文称为“双阶法”。显然，为了得到两套不同精度的解答，这两种方法都需要作两次有限元求解，计算代价过大成为其难以让人接受的最大弊端；此外，各自还有其他不利因素，将在后文讨论。

本文提出降阶单元自适应分析方法，对其初步的研究表明，该法可以克服上述所有缺点和弊端。降阶单元的概念和作法首先由文献[15]在初值问题中提出，是对初值问题中凝聚单元的逆构思和逆运用[16-17]，现已成为结构动力计算中的一种新型高效、可按最大模自适应步长的时程单元。本文在此基础上，进一步将其推广到边值问题，其基本作法是：采用m+1次常规单元做常规求解，得到常规有限元解后，略掉其m+1次项，得到m次的降阶解，作为降阶单元的解答，将略掉的m+1次项转而作为降阶解的误差估计器。可见，降阶单元只经一次有限元求解便得到两套不同精度的解答，而两套解在精度上的“阶差”很自然地提供了一个可作逐点估计的误差估计器，也自然地发展出以降阶单元作为最终解答的自适应算法。理论和算例均表明，该算法思路精巧、理论清晰、算法简单、实施灵活，既可靠又通用。

本文暂只限于介绍降阶线性元(m=1)及其自适应算法，文中以常规的二阶非自伴两点边值问题为模型问题进行推导和说明，也给出二维Poisson 方程和弹性力学平面问题的数值算例，以展示本文方法的广泛适用性和有效性。

1 问题描述及有限元解

1.1 模型问题及有限元解

考虑如下二阶非自伴两点边值问题：

式中，p、r、q均为x的函数，且p≥p0＞0，p0为常数。定义与问题(1)相应的双线性型和线性型为：

采用Galerkin 法求解式(1)所示问题的近似解，可归结为求解u∈HE1，使得：

式中：u、v分别为试探函数和检验函数；HE1为满足本质(位移)边界条件直到一阶导数均平方可积的函数空间。

1.2 常规二次元

不失一般性，考虑典型区间x∈[0,h]上的单元，h为单元长度。按照常规做法，将试探函数uh和检验函数vh均取为二次多项式；又为了后期方便，将形函数取为升阶谱形式：

式中：

注意：Nˆ3是一个两端点为0 的二次“泡状”函数；相应的，uˆh3为广义结点位移，并不代单元中点位移。则Galerkin 有限元解归结为求解，使得：

2 降阶线性元及其特性

2.1 降阶线性元的结构

在得到上述二次单元的解答后，可将各个单元的解答写为：

式中：

图1 一维降阶单元示意图Fig.1 Schematic diagram of one-dimensional reduced element

图2 二维降阶单元示意图Fig.2 Schematic diagram of two-dimensional reduced element

2.2 降阶线性元的特性

降阶线性元解并不等价于常规线性元解。根据一维有限元的数学理论[11-14]，常规线性元解、二次元解及降阶线性元解在单元内部和端结点处的收敛阶，归纳于表1，由表1 可见并可推断出如下特性：

表1 一维常规单元和降阶单元位移解的收敛阶Table 1 Convergence orders of displacement solutions for one-dimensional conventional and reduced elements

1) 降阶线性元的解是线性元和二次元混合解：内部是线性解，结点处是二次单元的解。

2) 降阶线性元的解，在单元内部的精度为h2阶，在结点处的精度为h4阶，可称之为超线性解。

3) 结点处h4阶的超收敛精度，可以非常有效地优先大幅减小结点位移的误差，将误差集中在单元内部，使得误差模式得到统一。故在做误差估计时，只估计单元内部误差即可，基本上可以排除结点误差的影响。

4) 二次单元内部解uh具有h3阶精度，比降阶线性解高一阶，用来逐点估计uhR的误差具有理论上的合理性和有效性，其误差恰为，堪称为预先内置了最大模估计器。

综上，降阶单元的核心思想是：常规二次单元的解，包含了一个降阶线性元解和一个误差估计项。用降阶单元的解作为最终解，用内置的误差估计项估计误差，指导网格细分，即得到一个十分简单方便的自适应有限元求解策略。

3 自适应算法

3.1 自适应求解目标及算法

本文的自适应求解的终极目标为：由用户事先给定误差限tol，寻找一个优化的有限元网格，使得该网格上的最终有限元解uh按照最大模度量满足误差限，即逐单元满足：

注意，按最大模自适应得到的是一个逐点满足误差限的解答，可看作在误差限tol意义下的数值精确解。

则，实际计算时，停机准则为逐单元满足：

可见，降阶单元的误差估计既不用逐点搜索，也不用逐点计算，简单、精确、高效。

至此，降阶单元的自适应求解算法可归纳为：

1) 给定初始网格和误差限tol。

2) 求二次有限元解uh，提取降阶线性元解。

3) 逐单元检验式(11)是否成立?

4) 对不满足式(11)的单元，将其二等分(一维单元)或四等分(二维单元)，返回到2)；

5) 直至所有单元都满足式(11)，则得到满足tol的降阶线性元解。

本小节将对本文方法与引言中提及的双阶法和双元法分别作一简要分析和比较。

双阶法：在同一网格下，分别用m次和m+1次的常规单元进行求解，取m次单元的解作为最终解，用m+1次单元的解来估计其误差，指导自适应过程。与降阶法相比，此法存在两点主要弊端：一是需要在同一网格上进行两次求解，计算代价过大(降阶法只需一次求解)；二是最终有限元解的端结点精度没有提高(降阶单元至少提高一阶)，导致误差模式纷杂，对某些问题自适应过程不稳定、最终单元数会无序增加(见第4 节的例1)。

双元法：给定一个初始网格，对其所有单元进行二分加密后得到第二个网格，对两重网格都采用m次常规单元进行求解，将稀疏网格的解作为最终解，用密集网格的解来估计其误差，指导自适应过程。与降阶法相比，此法存在3 点主要弊端：一是需要进行两次求解，计算代价过大；二是高次元二分后增加的自由度几乎翻倍(一个三次元二分后增加3 个自由度)，进一步增加了计算代价；三是检验解和最终解是同阶收敛，简单直接的误差估计并不可靠。

相比之下，无论是计算效率还是计算精度，降阶法比双阶法和双元法都具有明显优势，是更佳选择。降阶法巧妙地将两次求解合二为一，在提升了降阶解结点精度的同时，为各个单元内置了各自的误差估计器，自给自足、灵活方便。

4 数值算例

本文采用Fortran90 所编写的程序代码计算一维和二维例题，以验证并展示本法的有效性和可靠性。为方便起见，本文引入“误差比”，即误差与误差限之比，记作又记Ne为单元数，Nadp为自适应步数。例1.二阶非自伴两点边值问题

问题描述如下：

其精确解如图3 所示，表达式较复杂，不再给出。为了展示降阶法优于双阶法，本例中也给出双阶法的部分结果，其误差比记作。给定误差限tol=1×10-4，以均匀分布的16 个单元作为初始网格，计算过程和结果如图4～图9 所示。

图3 例1 的位移精确解Fig.3 Exact displacement solution for example 1

图4 降阶法误差比( Ne=16,Nadp=0)Fig.4 Error ratio of reduced order method(Ne=16,Nadp=0)

由所示结果可见，双阶法出现单元端结点误差大于内部误差的情况(如图5、图7 所示)，误差并非集中在单元内部，使得网格的调整次数增多，最终经过8 步自适应调整，共划分349 个单元作为最终网格(图9)。而降阶法在自适应过程中始终保持着如图4、图6 所示的稳定且统一的误差模式，即单元端结点误差远小于内部误差，最终仅经过5 步自适应调整，共划分123 个单元作为最终网格(图8)。可见两法虽然都给出了满足误差限要求的解答，但降阶法相比于双阶法更加简单、高效。例2.Poisson 方程——四边形区域问题

图5 双阶法误差比( Ne=16,Nadp=0)Fig.5 Error ratio of double order method ( Ne=16,Nadp=0)

图6 降阶法误差比( Ne=32,Nadp=1)Fig.6 Error ratio of reduced order method(Ne=32,Nadp=1)

图7 双阶法误差比( Ne=32,Nadp=1)Fig.7 Error ratio of double order method ( Ne=32,Nadp=1)

图8 降阶法误差比( Ne=123,Nadp=5)Fig.8 Error ratio of reduced order method(Ne=123,Nadp=5)

图9 双阶法误差比( Ne=349,Nadp=8)Fig.9 Error ratio of double order method(Ne=349,Nadp=8)

问题描述如下：

给定问题精确解为：

式中：

荷载f由原方程导出。四边形区域及初始网格如图10 所示，给定误差限tol=1×10-3。

图10 四边形区域示意图及初始网格Fig.10 Schematic diagram of quadrilateral area and initial mesh

本例有意给定图10 所示的 4×4非结构化网格作为初始网格，用以检验本法对非结构化网格的适用性。计算中经过5 步自适应调整，得到如图11所示的共2383 个单元的最终网格。图12 给出初始网格下的误差比分布，可见在单元角结点处误差比最小，即结点位移具有更高精度，符合降阶单元的预期误差模式，也是优势之一。图13 给出最终网格下降阶线性元的误差比分布，可以看出在 ±1以内。本例显示本法适用于非规则区域、非结构化网格，这是本法的一个突出特色，为实际应用提供了极大的方便灵活性。例3.Poisson 方程——SPR 算例

图11 例2 最终网格划分Fig.11 Final mesh of example 2

图12 例2 初始网格域内误差比的分布Fig.12 Error ratio on the initial mesh of example 2

图13 例2 最终网格域内误差比的分布Fig.13 Error ratio on the final mesh of example 2

此问题是SPR 自适应方法用于检验算法的算例[18]，问题描述如下：

给定问题精确解为

荷载f由原方程导出。此问题沿求解域对角线存在应力集中现象。初始网格取1 个单元，给定误差限tol=1×10-3。

计算中经过6 步自适应调整，得到如图14 所示的共763 个单元的自适应最终网格，从中可见，网格逐渐向对角线附近加密，显示出问题的难点所在。图15 给出最终网格下降阶线性元解的误差比分布，可以看出误差比都在 ±1以内。本法对存在一定高梯度的问题可以给出满足误差限要求的解答。

图14 例3 最终网格划分Fig.14 Final mesh of example 3

图15 例3域内误差比?的分布Fig.15 Error ratio of example 3

例4.弹性力学平面问题——Cook 梁

Cook 梁问题是由Cook 首先提出的一个经典考题[19-20]，它是如图16 所示的一端受剪力作用的变截面梯形悬臂深梁，本例的特点是A点的应力奇异。按平面应力问题进行纯数值求解，给定弹性模量E=1 ，Poisson 比ν=1/3 ，厚度t=1，剪力P=1 ， 误差限tol=5×10-2。此外，右端剪力P等效地转化为沿截面为二次抛物线分布，以避免在上、下角点出现与理论不自洽。

图16 Cook 梁示意图及初始网格Fig.16 Schematic diagram of Cook beam and initial mesh

计算中以如图16 所示的2 个单元作为初始网格，经过6 步自适应调整，得到如图17 所示的共251 个单元的最终网格，从中可见，网格逐渐向A点附近加密，显示出问题的难点所在。图18、图19 分别给出最终网格下两个位移分量的误差比分布，可以看出都在 ±1以内。本法对弹性力学平面问题可以给出满足误差限要求的解答。

图17 例4 最终网格划分Fig.17 Final mesh of example 4

图18 例4域内误差比(u)的分布Fig.18 Error ratio (u)ofexample 4

图19 例4 域内误差比(v)的分布Fig.19 Error ratio(v) of example 4

5 结论

理论分析和大量的数值试验表明，降阶法相较于其他自适应方法有如下几点优势：① 只需求解一个二次单元的解答，其包含所需的最终解(降阶单元解)和一个误差估计器；② 可以有效减少自适应所需自由度数；③ 有限元解的误差模式得到统一，即结点误差相较单元内部误差为高阶微量，且各单元误差单调减小，极少出现单元反叛情况；④ 算法简单，便于实施，可靠通用。

本文以二阶非自伴两点边值问题为例，在初值问题降阶单元的基础上，进一步提出边值问题的降阶线性元及其自适应算法。降阶单元具有高精度的结点解，同时内置了有效可靠的误差估计器，从单元构造、数值实施来看，也更加简单、便捷、灵活、通用，具有明显的优势。