熟知，运动方程在数学上归结为二阶常微分方程(组)初值问题，常用中心差分法、Newmark-β法、Wilson-θ法等逐步长积分方法来求解，而数值解答的收敛性、精度、数值稳定性、计算效率等都强烈依赖于时间步长的选择和优化[1－2]。本文知难而进，把自适应步长求解作为目标。
有限元法(FEM)是一种通用有效的数值计算方法，在工程计算中得到极为广泛的应用[3]。有限元中最简单、灵活、方便的单元是线性单元。但是，线性元只有2阶收敛、位移精度不高、导数精度更差、缺乏超收敛性态等不利因素，使得线性元常常被排除在高精度高效能计算之外。本文返璞归真，简之又简，采用Galerkin型线性有限元。
要实现上述两点，即对运动方程实现线性有限元的自适应步长求解，需要有一个核心关键技术作为前提，即：用来估计误差、调整步长、修正精度的超收敛计算技术。袁驷等[4－6]基于数学和力学概念，提出了一种有限元后处理超收敛算法——单元能量投影(Element Energy Projection，简称EEP)法，且已经建立了一套完整、统一的 FEM 自适应求解算法，已对各类一维、二维乃至三维线性、非线性边值问题取得了一致成功[7－9]。最近，采用与边值问题类似的高次(>1)单元和整域网格细分的方法，还成功地应用于运动方程的自适应求解[10－11]。
本文深入发掘并用巧用足 EEP超收敛计算技术，构建了Galerkin型线性元的单步法递推公式，开发了单元步长自动优化和结点位移精度修正两项关键技术，可在整个时域上得到误差分布均匀且逐点按最大模度量满足给定的误差限的解答——堪称数值解析解。文中给出了单自由度和多自由度的数值算例以验证本法的有效性。
1 问题描述和Galerkin有限元法
1.1 问题描述
为了表述方便，本文仅以单自由度系统为例展开讨论；所述方法同样适用多自由度系统，后文将给出其具体的数值算例。
单自由度系统的运动方程可以归结为如下的二阶常微分方程初值问题[1－2]：
其中：u、
分别为位移、速度和加速度；m、c和k分别为质量、阻尼和刚度(均为常数)；P为外荷载且是时间t的函数；L是式中定义的微分算子。
为构造以上问题的Galerkin法弱形式，定义双线性型和线性型：
其中，
为所有满足u(0)=u0且直到一阶导数均平方可积的函数空间，而
为所有满足v(T)=0且直到一阶导数均平方可积的函数空间。可以验证，式(3)等效于原问题式(1)。
1.2 Galerkin有限元法——线性单元
先按整体求解的思路，将定义域(0,T]划分为Ne个单元，整体结点编码从定义域最左端为结点0，依序递增，最右端为结点Ne。考虑典型单元e，两端点坐标为t1、t2，长度为h。单元上的试探函数和检验(权)函数均由线性插值得到：
则问题式(1)的有限元解归结为求解
使得：
其中，
为通常意义下的线性有限元试探空间和检验空间。记u、uh和vh为广义精确解、有限元解和任一检验(权)函数，则有限元解的误差eh=u-uh满足投影定理[3]：
本文后边将把投影定理近似地用在单个单元上，导出EEP超收敛计算公式。
基于式(5)，可以求得单元刚度矩阵的显式表达式：
按照有限元常规的做法进行整体集成，并引入初始条件后，可以得到整体刚度方程。由初值问题的特性所决定，整体刚度矩阵是一个下三角矩阵，并不需要整体联立求解，可采用逐步长递推求解，其递推公式为：
其中，上标(i)为单元码。以上是本法逐步长求解的基本递推公式。
1.3 与其他方法比较
为了加深对本法的理解和认识，这里将其与两个常用的中心差分法和 Newmark-β线性加速度法做一简单比较。
中心差分法的基本公式为[2]：
为了与其比较，将刚度系数的表达式(7)代入式(11)，得到本法相应的计算公式：
可见，本法同中心差分法相比只在含k的项和右端荷载项处理不同。换言之，本法公式可以看作是利用Galerkin法构造的另一种差分格式。
Newmark-β线性加速度法是一个自起步的单步法。用该法计算第 1个结点的位移，并利用本法式(9)的符号，可以得到该法的计算公式[2]：
可见，除了荷载项处理不同外，该法与本法的差别只是多了一些二阶微量(含h2)的项。由于该法也只有二阶收敛，因此所有含h2的项基本上是冗余计算，将其略掉并不降低其收敛阶。换言之，本法可以看作是简化的、线性化的Newmark-β线性加速度法。
此外，经计算，本法的稳定性条件优于中心差分法(是其
倍)，与 Newmark-β线性加速度法相同。总之，本文方法在简单方便性上与上述两种方法大致相当，但又具有下文将要介绍的一系列独特优势。
1.4 单步法递推格式
从本法第1个结点的递推公式(9)可以看出，本法是显式格式的、自起步的单步法：即给定初始位移u0和速度
，便可得到下一结点(时刻)的位移。但从后续结点的递推公式(10)来看，每个步长(单元)的求解，都需要两个单元3个结点的位移值，表面上看不再是单步法了。其实不然。若将式(10)中所有单元(i-1)的项都放在一起，即
，则发现它恰恰代表单元
的第2个端点处速度
的近似解，而且是下文将要介绍的EEP超收敛解，特将其记为
。后文将证明
与结点位移
具有同阶精度。如此对式(10)做简化，可得到更加明确的单步法递推公式：
其中，在第一个单元
。此公式非常简捷明了：每计算完一个单元步长，就把该单元终点的位移
和速度
同时计算出来，作为下一个单元起点的初始条件，而无须保留当前单元的任何其他(刚度系数和荷载系数)信息。
2 自适应步长的确定
2.1 EEP超收敛解公式
EEP(Element Energy Projection，单元能量投影)超收敛计算是本文的核心技术[4－6]。在得到常规有限元解答
后，假设投影定理式(6)在单个单元e上仍近似成立，利用分部积分并作一些略舍即可推导出EEP超收敛位移和导数的计算公式如下[5]：
其中：ta为单元[t1,t2]上的任一点，()a表示在点ta处取值；上标*表示超收敛解。可以看出，EEP解
的计算只涉及一些积分运算，对于线性元，甚至可以推出显式表达式；这也是采用线性元的另一个优点。
数学上已经证明，对于
次单元，用式(15)和式(16)计算的
均可达到
阶收敛，比有限元解
阶至少高一阶。但是，线性元是个例外，其结点位移
和单元内部位移
阶收敛，用式(15)和式(16)计算的
也只有
阶收敛；亦即，线性元的超收敛解
并没有“超过”有限元解uh的收敛阶，这也是线性元一直被排除在EEP超收敛计算之外的原因之一。
但是，线性元的u?在单元内部的数值精度仍有很大的改善，仍然可以用于检验有限元解在单元内部的误差。为了具体展示超收敛解的效果，这里给出一个数值算例；采用的数据是下文第4节的例1，只计算一个单元，步长h=0.2。各种位移的误差如图 1所示，可见用 EEP解计算的有限元解的误差u*-uh和真实的误差u-uh极为接近。同时还可以看出，有限元解与超收敛解在单元端点的值相同，因此，EEP超收敛解无法用来估计有限元端点位移的误差；换言之，用EEP解估计单元内的误差，其可靠性是以高精度的有限元结点位移为前提的。超收敛导数见图2。可以看到，EEP超收敛导数对有限元解有数量级的改善，而且精确满足左端点给定的初始速度条件。
图1 有限元误差与超收敛误差对比(位移)
Fig.1 Errors in FE and EEP solutions (Displacement)
图2 EEP超收敛导数图
Fig.2 EEP derivatives
至此，可以证明式(14)中的
恰恰是该结点处的EEP超收敛导数。在式(16)中，将ta取在单元的两个端点，注意到
，再利用分部积分，不难得到：
式(17)中的第2个方程即为式(14)中
的表达式。
2.2 自适应步长的计算公式
本文的自适应求解的终极目标为，在精确解u未知的情况下，事先给定误差限Tl，寻求一个优化的单元网格分布，使得所有单元上有限元解答uh按最大模度量满足给定的误差限Tl，即逐单元满足
。实施中，由于精确解u未知，使用EEP超收敛解u*代替u来估计误差，即要求有限元解uh逐单元满足：
现有的自适应算法中，大都是仿照边值问题的做法，对整体网格进行细分，已经插入的结点不再改变[10―11]。这种做法的好处是初值和边值问题通用，但是需要联立求解整体刚度方程并对整体网格细分，有时会产生过大的冗余精度和冗余单元数目。事实上，对于初值问题，本可以遵循其特性，从初始时刻起对每一单元步长都作自适应调整优化，逐步长地递推求解。
这里，一个关键技术是如何调整单元步长。基于有限元数学理论，本文给出一个简单、有效且实用的计算公式(具体的数学分析和推导从略)。设已给定步长
并已求得其上的解uh和u*，而新的单元长度h应尽量使得
，则其计算公式为：
其中，α≤ 1是一个安全系数。数值试验也验证了2/5的分数阶次是最优的。
如此，本文的逐单元调整步长的自适应算法似可初步归结为“三步走”策略：1) 有限元解；2) 超收敛解；3) 调整步长。这一算法，对于线性元来讲，并不完善。实际计算时，在一定步数之前都是有效的；之后，结点和单元内点位移的误差都开始超限。这是由于线性元的结点与单元内点位移精度都是2阶收敛的；而对于初值问题(如无阻尼的自由振动)，即使步长满足稳定性条件，结点位移的误差往往仍不可避免地在相当长的时间区域上逐步积累增大，随着结点位移误差的超限，单元内部EEP超收敛解的精度随之失去保障，误差估计失效。这里，一个严峻的挑战是结点位移的误差控制，这也是各类数值方法的共同难点。下面将引入基于EEP超收敛计算的结点位移精度修正技术来破解这一难点。
3 结点位移精度的修正
3.1 结点位移误差的估算
记超收敛解误差e*=u-u*。由于超收敛解u*与有限元解uh在结点上相等
因此在结点上误差相同
如果能够求解出超收敛解误差e*，则可得到有限元结点位移
的误差，对
进行修正，就得到更高精度的结点位移。下面仍然用线性元求解e*的有限元近似解εh[10]。
将e*=u-u*代回原问题式(1)，可以得到e*的控制微分方程和初值条件：
其中，残余荷载P*=P-Lu*。注意，由于u*满足位移和速度的初始条件，因此式(20)中初始条件均为齐次条件(右边为 0)。显而易见，求e*的式(20)与求u的式(1)除了荷载和初值不同外，其他所有定解条件都相同。这就意味着，在原有限元网格上求解e*的近似解εh，其整体刚度矩阵及其逆矩阵均无改变，只是荷载向量需要重新生成，再经一个回代过程即可。换言之，结点位移误差的计算仅相当于多计算一种荷载工况。
对于本文的单步长方法，具体实施步骤如下：
1) 在当前单元(步长)，通过递推公式(14)计算有限元解uh；
2) 用式(15)计算EEP超收敛解u*；
3) 如需要自适应调整步长，用式(19)调整步长，直至当前单元自适应求解完毕；
4) 将u*代入P*=P-Lu*，进行荷载线性型积
分，得到单元荷载向量
；
5) 用递推公式(14)(用
代替Pe)得到结点位移误差
和导数误差
；
6) 将
分别迭加到
上，得到修正的结点解
大量数值计算表明，如此修正后的结点位移由原来的h2阶收敛提升为h4阶收敛，可谓取得了难得的“翻倍”超收敛的效果！值得一提的是，由于本法是逐步长修正(用单元左端的修正解，求右端有限元解，再修正右端解)，它比整体修正(在整体区域上对有限元各结点解同时修正)的精度又有翻倍提高的效果。这是逐步长求解法的又一优势。
3.2 算法小结
综合以上的零散算法，可以构成一个用于求解运动方程的简捷、高效、可靠的自适应步长线性有限元分析方法。其中，EEP超收敛解起到了关键的作用，用它代替精确解，估计单元上有限元解的误差，估计最佳步长，还可以对结点位移的精度进行大幅修正。
为方便，本文引入“误差比”：即最大误差与误差限之比，用最大误差加上划线来表示。记
，则误差比记为
其值≤1即满足了误差限。本算法的突出特点是自适应调整步长。在实际计算中，最佳步长虽然可以用公式迭代算出，但考虑到超收敛解毕竟也是近似解，要留有一定的安全空间，所以误差限Tl用稍加保守的误差限代替(α·Tl)，其中安全系数一般取为α=0.8即可；另外，对步长的精度也不宜过分苛求，只要误差比在一个可接受区间内即可，本文取
为可接受区间。大量数值算例表明，很多单元都无需调整步长，而需要调整时，最多迭代2次即可调整到位。
综合起来，本文的完整的自适应算法归纳为：
1) 给定误差限Tl和初始步长h0；
2) 在当前单元上用递推公式(14)求有限元解uh；
3) 用式(15)计算 EEP超收敛解u*，计算单元内最大误差比
4) 如果
，则用式(19)调整步长，返回到2)(一般调整一两次即可)；
5) 如果
，对结点位移进行修正：
a) 将u*代入P*=P-Lu*，进行荷载线性型积分，得到
；
b) 用
代替Pe，用递推公式(14)得到结点位移误差
和导数误差
；
c) 将
分别迭加到
上，得到修正后的结点解
。
6) 返回2)；直至到达或超过终点时间T。
4 数值算例
根据以上算法，本文用Fortran 90编制了相应的程序代码。限于篇幅，本节给出两个算例，分别为单自由度(图3)和多自由度算例(图4)，且均采用无量纲计算。每个算例都给出用精确解计算的各个结点位移和单元上位移的误差比
的分布图，可见所有结点和单元都满足精度要求，自适应求解是成功的。
例1.有阻尼单自由度体系的简谐荷载反应
求解数据为m=k=1，c=0.04，
，误差限Tl=1/1000，精确解如图 3所示。首先用定步长h=0.2计算，并与中心差分法比较，如图5所示。可看到，有限元法起始一段的误差小于差分法，在之后的则与中心差分法基本相同。还可以看出，结点位移的误差比很快(几秒左右)就超过 1，大约在50 s前后达到峰值，然后开始衰减。这个解答的误差远远超出了误差限Tl=1/1000。若加以结点位移修正，结点位移误差比分布如图6所示，最大的误差比
由无修正的 24.4降到修正后的0.0955，大约是256倍。
图3 例1的解析解u
Fig.3 Exact solutionufor Example 1
图4 例2的解析解u1,u2,u3
Fig.4 Exact solutionsu1,u2andu3for Example 2
图5 有限元与
差分法结点位移误差比
对比图
Fig.5 Error ratios for nodal displacements of FE and FD
图6 修正后
的有限元结点位移误差比
分布
Fig.6 Error ratios for recovered FE nodal displacements
再用完整的自适应步长算法求解，初始步长h0=1，共用了1003个单元，共调整了472次步长，不到50%的单元做了步长调整，最大单元、最小单元长分别为hmax=1.967和hmin=0.093。结点位移和单元位移误差比的分布示于图7和图8，可见满足精度要求，结点位移保持高精度，且各单元误差分布比较均匀。
例2.多自由度体系有阻尼受迫振动
本例意在展示本文方法对多自由度系统的适用性。物理模型来源于3层剪切型框架结构，计算数据为：
图7
例1的结点位移误差比
Fig.7 Error ratios for nodal displacements of Example 1
图8 例1
的单元位移误差比
Fig.8 Error ratios for element displacements of Example 1
简谐荷载向量为F=(sin(10t)sin(10t)sin(10t))T，初始位移和初始速度均取零，T取10 s，解析解如图4所示。计算时初始步长
，误差限T=10-4。共用了 328l个单元，共调整了45次步长，每个步长最多调整2次，不到50%的单元做了步长调整，最大单元、最小单元长分别为hmax=0.1104和hmin=0.0195。3个位移分量的结点位移和单元位移误差比的分布示于图 9～图 12，可见圆满满足精度要求，结点位移保持高精度，且各单元误差分布均匀。
图9 例2中
u1的结点位移的误差比
Fig.9 Error ratios for nodal displacements ofu1of Example 2
图10 例2中u2的结点位移的误差比
Fig.10 Error ratios
for nodal displacements ofu2of Example 2
图11 例2中u3的结点位移的误差比
Fig.11 Error ratios
for nodal displacements ofu3of Example 2
图12 例2
的单元位移误差比
Fig.12 Error ratios for element displacements of Example 2
5 结论
本文是一个初步的尝试，也是一个成功的进展，算法本身有进一步优化的空间，方法具有广泛的应用背景和发展前景。
参考文献：
[1]Clough R W. Dynamics of structures [M]. New York:McGraw-Hill, 1995: 68.
[2]刘晶波, 杜修力. 结构动力学[M]. 北京: 机械工业出版社, 2005: 132.Liu Jingbo, Du Xiuli. Structural dynamics [M]. Beijing:China Machine Press, 2005: 132. (in Chinese)
[3]Strang G, Fix G J. An analysis of the finite element method [M]. Englewood Cliffs, N J: Prentice-Hall, 1973:1.
[4]袁驷, 王枚. 一维有限元后处理超收敛解答计算的EEP法[J]. 工程力学, 2004, 21(2): 1―9.Yuan Si, Wang Mei. An element-energy-projection method for post-computation of super-convergent solutions in one-dimensional FEM [J]. Engineering Mechanics, 2004, 21(2): 1―9. (in Chinese)
[5]袁驷, 林永静. 二阶非自伴两点边值问题 Galerkin有限元后处理超收敛解答计算的EEP法[J]. 计算力学学报, 2007, 24(2): 142―147.Yuan Si, Lin Yongjing. An EEP method for post-computation of super-convergent solutions in one-dimensional Galerkin FEM for second order non-self-adjoint boundary-value problem [J]. Chinese Journal of Computational Mechanics, 2007, 24(2): 142―147. (in Chinese)
[6]Wang Mei, Yuan Si. Computation of super-convergent nodal stresses of Timoshenko beam elements by EEP method [J]. Applied Mathematics and Mechanics(English Version), 2004, 25(11): 1228―1240.
[7]Yuan Si, He Xuefeng. A self-adaptive strategy for one-dimensional FEM based on EEP method [J]. Applied Mathematics and Mechanics (English Version), 2006,27(11): 1461―1474.
[8]Yuan Si, Xing Qinyan, Wang Xu, et al. Self-adaptive strategy for one-dimensional finite element method based on EEP method with optimal super-convergence order[J]. Applied Mathematics and Mechanics (English Version), 2008, 29(5): 591―602.
[9]袁驷, 徐俊杰, 叶康生, 等. 二维自适应技术新进展:从有限元线法到有限元法[J]. 工程力学, 2011, 28(增刊II): 1―10.Yuan Si, Xu Junjie, Ye Kangsheng, et al. New progress in self-adaptive analysis of 2D problems: from FEMOL to FEM [J]. Engineering Mechanics, 2011, 28(Suppl II):1―10. (in Chinese)
[10]邢沁妍. 基于EEP法的一维Galerkin有限元自适应分析[D]. 北京: 清华大学, 2008.Xing Qinyan. Adaptive analysis of 1D Galerkin FEM based on EEP super-convergent method [D]. Beijing:Tsinghua University, 2008. (in Chinese)
[11]邢沁妍, 杨杏, 袁驷. 离散系统运动方程的Galerkin有限元 EEP法自适应求解[J]. 应用数学和力学, 2017,38(2): 133―143.Xing Qinyan, Yang Xing, Yuan Si. An EEP adaptive strategy of Galerkin FEM for dynamic equations of discrete systems [J]. Applied Mathematics and Mechanics (Chinese Version), 2017, 38(2): 133―143. (in Chinese)