结构动力响应问题是近年来的热门课题，常用的时程积分方法更是常用的重要方法，但其计算效率、收敛精度、数值稳定性等均与步长的选择和优化相关^[1-2]，因此自适应步长求解也成为近几年研究的难点和热点之一。就自适应求解而言，有对离散解进行误差估计并建立自适应算法^[3-4]的研究，但在时域上逐点按最大模控制误差的高效自适应求解比较少见，而这恰是本文追求的目标。

袁驷教授领导的课题组依据有限元原理^[5]提出的单元能量投影^[6-8](element energy projection，EEP)法是一种非常有效的有限元后处理的超收敛算法，可用于逐点误差估计和精度修正。该EEP算法已应用于一维^[9-10]、二维有限元^[11]，乃至三维有限元^[12]，并且提出一类基于EEP技术的线性有限元时程积分的自适应方法^[13-14]。该法采用简单的线性元，构建单步法递推公式，然后运用EEP法进行超收敛计算，并对结点位移修正来提高精度并进行后验自适应求解。此外，该课题组还对边值问题提出了一种先验的自适应方法，称为固端法^[15-16]，该法可以不经有限元计算，通过EEP解的超收敛性进行直接误差估计，并一举给出满足精度要求的网格划分。

本文基于边值问题的固端法思想，提出一类不经有限元求解、便能得到先验EEP解的自适应步长算法。本文算法相比于文献[13]的后验法，能够在不计算有限元解的情况下，得到满足精度要求的单元步长。文中以线性元为例，给出了代表性算例，用以验证该法的有效性和高效性。算例表明，本文方法在大幅提高计算效率的同时，仍严格保障了解答精度，为其后相关的深入广泛研究提供了有效的可行性方案。

1   问题描述与有限元解

文献[13]采用线性Galerkin有限元，对运动方程进行自适应步长求解提出了整套方案，本节对此做一简要概述。

1.1   问题描述

本文考虑单自由度常系数运动方程：

式中：$L$为式中定义的微分算子；$m$、$c$和$k$分别为质量、阻尼和刚度；$f$为外荷载函数；$u$为动位移；${u_0}$、${v_0}$分别为初始位移和速度；$\bar T$为时域的上界。为了兼顾一般问题，时间自变量采用$x$。

为构造Galerkin弱形式，定义双线性型和线性型：

记$H_u^1$为所有满足$u(0) = {u_0}$且直到一阶导数均平方可积的函数空间，而$H_v^1$为所有满足$v(\bar T) = 0$且直到一阶导数均平方可积的函数空间，则问题(1)(见式(1))的Galerkin法归结为求解$u \in H_u^1$使得：

1.2   Galerkin有限元解

由于可以逐单元步长积分求解，不失一般性，仅考虑两端结点坐标为$(0,\;h)$的典型单元$e$。采用线性单元，将试探函数$u_{}^h$与检验函数$v_{}^h$均取为线性函数：

则Galerkin有限元解归结为求解${u^h} \in S_u^h \subset H_u^1$使得：

其中，$S_u^h$和$S_v^h$为通常意义下的有限元试探空间和检验空间。基于式(5)，可求得单元刚度矩阵和单元等效结点荷载向量：

则Galerkin有限元求解归结为如下递推公式^[13]：

其中，在第一个单元，$u_1^h$和$v_1^*$分别取${u_0}$和${v_0}$。注意，${m^{ - 1}}$与步长$h$无关，可以一次性求出，而$k_{12}^{ - 1}$则对不同步长$h$需要重新计算，对于多自由度问题更是主要的耗时计算的环节。通常情况，线性元的结点位移和单元内位移同为$O({h^2})$收敛。

1.3   EEP超收敛解

EEP超收敛计算是本文误差估计的核心技术^[6-8]。在得到常规有限元解答${u^h}$后，可推导出EEP超收敛位移的计算公式^[13]如下：

其中：上标$*$表示超收敛解。对于线性元，有限元解${u^h}$与EEP解$u_{}^ *$具有相同收敛阶，同为$O({h^2})$；但研究表明，$u_{}^ *$仍可用来代替精确解$u$对${u^h}$进行误差估计。

1.4   自适应步长公式

本文的自适应目标与文献[13]相同，即在精确解$u$未知的情况下，事先给定误差限$tol$，寻求一个优化的单元步长$h$，使得有限元解答${u^h}$按最大模度量，逐单元满足：

其中，${u^h}$和$u_{}^*$为当前单元步长${h_0}$上的解答。若当前解不满足式(9)，则采用文献[13]的步长公式估算新的步长：

其中，$\alpha \leqslant 1$是一个安全系数(通常取0.8)。步长公式的数学证明见文献[17]。

1.5   结点位移修正

由于线性元结点收敛阶与有限元解和EEP解收敛阶均相同，而对于初值问题，结点位移的误差是发生数值耗散的最主要和关键的因素，为此文献[13]引入了基于EEP超收敛解的结点位移修正技术^[18]。

将${e^*} = u - {u^*}$代回式(1)，限于第一个单元，可以得到${e^*}$的控制微分方程和初值条件：

其中，残余荷载${f^*} = f - L{u^*}$。对此问题，在现有的有限元步长上，再次求解${e^*}$的有限元解${\varepsilon ^h}$，并将结点2的位移修正为$u_2^h + \varepsilon _2^h$。这一修正技术的突出特点是，只需用${f^*}$代替$f$形成一个新荷载向量，无需求逆计算，只经回代过程即可；且如此修正后的结点位移的收敛阶由原来的$O({h^{\text{2}}})$提升为$O({h^{\text{4}}})$，使得结点位移精度翻倍，极大提升了自适应的性能和效率。

1.6   自适应步长求解

综合以上可得一个自适应步长求解的完整算法^[13]，简要归纳为：① 用初始步长计算有限元解；② 计算EEP解并进行误差比较；③ 如果不满足要求则调整步长，并返回①；④ 如果满足要求则对右端点进行结点位移修正，并进行下一个单元计算。

2   先验法

1.6节算法是一套完整、可靠、有效的算法，其中为了自适应调整步长，采用了后验误差估计$u_{}^ * - {u^h}$。为了计算$u_{}^ * - {u^h}$，需要事先计算有限元解$u_{}^h$，亦即需要多次计算式(7)中的$k_{12}^{ - 1}$。以下探讨一种不经有限元求逆计算$k_{12}^{ - 1}$的先验误差估计并确定步长的方法，以期进一步提升求解效率。

2.1   误差估计的简化

误差估计的常规做法是采用EEP解${u^*}$代替精确解$u$，亦即用${e^{*h}} = {u^*} - {u^h}$代替${e^h} = u - {u^h}$来估计有限元解${u^h}$的误差，且要求${e^{{\text{*}}h}}$满足误差限来确定临界步长。由于${e^{*h}}$在单元两端点均为0，故单元上${e^{{\text{*}}h}}$的最大值必取在单元内的极值点上。这样，确定临界步长归结为如下一组非线性代数方程组的求解，即求解$x、\;h$，使得：

为了能得到式(12)的先验计算结果，本文采取以下两项简化举措。

简化一，将单元上的荷载近似为常数${f_0}$。本文取单元内中点值，即${f_0} = f({h \mathord{\left/ {\vphantom {h 2}} \right. } 2})$，则${e^{*h}}$可简化为显式表达的三次多项式：

其中，$H = 6h{\kern 1pt} {k_{12}} = 6m + 3ch + k{h^2}$。

${\hat e^{{\text{*}}h}}$虽然是x的三次多项式，但是对于h是有理分式，其中${1 \mathord{\left/ {\vphantom {1 H}} \right. } H}$涉及${({k_{12}})^{ - 1}}$的计算，随着自由度数值的增加，计算效率随之大幅降低。事实上，若完成了${({k_{12}})^{ - 1}}$的计算，就相当于完成了有限元求解，便不能称为先验法了。

简化二，将${H^{ - 1}}$进行Taylor级数展开。记为：

实际计算时，$\bar H$的展开式中可取2项或者3项。

将$\bar H$展开式(14)代入式(13)，可得到先验误差估计${\hat e^{{\text{*}}h}}$的进一步简化的表达式。

2.2   先验误差估计

在计算步长的时候，首先需要求得最大误差。为了求${\hat e^{{\text{*}}h}}$的极值，给定初始步长${h_0}$，则由$({\hat e^{{\text{*}}h}})' = 0$可得一个$x$的二次方程，意味着在单元上最多有2个极值点(记为${x_1}$、${x_2}$，且${x_1} \leqslant {x_2}$)，均可以用公式直接求出。

记${x_{\rm m}}$为取得最大误差值$\max | {{{\hat e}^{{\text{*}}h}}} |$的坐标点。若单元内只有一个极值点，则该点即为${x_{\rm m}}$。若有两个极值点，则取使绝对误差较大的${x_{\rm m}}$。

2.3   先验步长计算

有了当前步长${h_0}$上的${x_{\rm m}}$，则可令$| {{{\hat e}^{{\text{*}}h}}({x_{\rm m}})} | = tol$来估算一个临界步长$h$。但由于临界步长的迭代计算涉及到虚根计算，仍显复杂，本文采用效率更高的步长公式计算新的步长：

其中，${\alpha _0}$为先验算法的安全系数，可以比后验公式(10)中的$\alpha$更小一些，以保证先验法的有效性和安全性。与式(10)类似，先验法自适应目标为$\max | {{{\hat e}^{{\text{*}}h}}} | \leqslant tol$；但是与后验法不同的是，先验法无需计算有限元解。

2.4   双验法求解策略

有了先验误差估计和步长计算的手段，可以对第1.6小节中的求解方案作出改进，即：用先验法估计步长，用后验法检验步长。这样，先验法和后验法各自发挥自己所长，既提高了求解效率，又不失精度的保障。这一“双验法”自适应步长求解策略可归纳为：

a) 给定初始步长${h_0}$；

b) 用先验法检验误差，必要时调整步长$h$；

c) 计算有限元解；

d) 计算EEP解并用后验法检验误差；

e) 如不满足要求则调整步长，并返回到c)，满足则继续下一步；

f) 对右端点进行结点位移修正，并进行下一个单元计算。

3   数值算例

本文采用符号计算软件Maple计算有阻尼简谐运动和无阻尼自由振动的算例，以验证本法的有效性。为方便起见，本文引入误差比，即误差与误差限之比，记作${\bar e^h}$。

例1. 有阻尼简谐运动。

计算数据为：$m = 1$，$\;k = 1$，$c = 0.04$，$f = \sin (0.2x)$，${u_0} = 0$，${v_0} = 1$，$\bar T = 256{\text{ s}}$，初始步长${h_0} = 0.5$。表1和表2分别是$tol = 0.001$和0.01的相关计算结果。需要注意的是，无论是先验法还是后验法，理论上都是基于$h$足够小时才足够有效，且$h$越小越可靠。表1中的误差限比较小，步长$h$足够小，先验法几乎可以取代后验法，而且所有单元逐点满足误差限，精度完全保证，效率大幅提高。表2中的误差限相对较大，使得$h$变大，后验法矫正次数增加，但直接后验法失控率反而更高；从误差比大于1的单元数和最大误差值来看，双验法比直接后验法更加安全可靠，可以更好地满足自适应要求，这也是先验法的一种功效。在图1和图2是取Taylor级数2项，$tol = 0.001$的先验法数据，可见该方法单元步长分布合理，误差得到有效控制。

表  1  有阻尼简谐运动计算数据 (256 s, $tol = 0.001$)

Table  1.  Results of damped harmonic motion (256 s, $tol = 0.001$)

误差估计方法	单元数${N_{\rm{e}}}$	自适应次数${N_{{\rm{adp}}}}$	先验迭代次数${n_{{\rm{iter}}\,0}}$	后验迭代次数${n_{{\rm{iter}}\,}}$	误差比大于1的单元数${N_{{\rm{fail}}}}$	最小单元长度${h_{\,\min }}$	最大单元长度${h_{\,\max }}$
直接后验	1003	373		472	0	0.09255	1.21070
求逆先验	1047	286	315	0	0	0.09050	0.75368
3项先验	1044	284	310	4	0	0.09050	0.76208
2项先验	1081	294	329	23	0	0.09048	0.63940

下载: 导出CSV 

| 显示表格

表  2  有阻尼简谐运动计算数据 ( 256s, $tol = 0.01$)

Table  2.  Results of damped harmonic motion (256s, $tol = 0.01$)

误差估计方法	单元数${N_{\rm{e}}}$	自适应次数${N_{{\rm{adp}}}}$	先验迭代次数${n_{{\rm{iter}}\,0}}$	后验迭代次数${n_{{\rm{iter}}\,}}$	误差比大于1的单元数(最大超限值)${N_{{\rm{fail}}}}$	最小单元长度${h_{\,\min }}$	最大单元长度${h_{\,\max }}$
直接后验	342	99		168	16 (1.384)	0.29528	3.18623
求逆先验	381	70	87	4	1 (1.005)	0.28098	1.30528
3项先验	364	85	128	20	2 (1.097)	0.28103	1.35135
2项先验	422	81	116	8	0	0.27923	0.99091

下载: 导出CSV 

| 显示表格

图  1  例1步长分布图

Figure  1.  Step-size distribution for Example 1

下载: 全尺寸图片 幻灯片

图  2  例1单元位移误差比图

Figure  2.  Displacement errors of elements for Example 1

下载: 全尺寸图片 幻灯片

例2. 无阻尼自由振动。

计算数据为：$m = 1$，$\;k = 1$，$c = 0$, $f = 0$，${u_0} = 0$，${v_0} = 1$，$\bar T = 256{\text{ s}}$, $tol = 0.001$, 初始步长$h = 0.5$。由表3可知，采用直接求逆的先验法和后验法是完全等效的；这是由于本算例中$f = 0$，则荷载近似为常数以及逆矩阵的Taylor级数展开带来的误差，都未出现。此外，由于本例步长较小，故Taylor渐近展开也足够精确，可以完全替代后验法。图3和图4为Taylor级数取2项的先验法计算结果；可见，本例步长、误差分布均匀且有周期性规律。

表  3  无阻尼自由振动计算数据 ( 256 s, $tol = 0.001$)

Table  3.  Results of undamped free vibration (256 s, $tol = 0.001$)

误差估 计方法	单元数${N_{\rm{e}}}$	自适应次数${N_{{\rm{adp}}}}$	先验迭代次数${n_{{\rm{iter}}\,0}}$	后验迭代次数${n_{{\rm{iter}}\,}}$	误差比大于1的单元数${N_{{\rm{fail}}}}$	最小单元长度${h_{\,\min }}$	最大单元长度${h_{\,\max }}$
直接后验	2790	677		677	0	0.07749	0.25854
求逆先验	2790	677	677	0	0	0.07749	0.25845
3项先验	2790	677	677	0	0	0.07749	0.25848
2项先验	2790	677	677	0	0	0.07749	0.26791

下载: 导出CSV 

| 显示表格

图  3  例2步长分布图

Figure  3.  Step-size distribution for Example 2

下载: 全尺寸图片 幻灯片

图  4  例2单元误差比图

Figure  4.  Displacement errors of elements for Example 2

下载: 全尺寸图片 幻灯片

4   结论

本文以单自由度运动方程线性元为例，提出了一种先验的自适应步长算法，构建了“双验”法自适应求解策略：

(1) 该法对于给定的求解误差限，可以直接确定步长大小，从而得到合理的有限元解；

(2) 在后验法保障精度的前提下，简明高效的先验法大幅提升了自适应步长求解的效率；

(3) 典型的数值算例验证了本法的有效性，展示了本法在计算效率上的提升。

本研究为未来多自由度运动方程、高次单元等更为一般情况的自适应步长求解提供了良好的基础理论和成功的案例。