自适应步长时程分析的新进展：从凝聚单元到降阶单元

2023-01-04袁驷，袁全

工程力学 2023年1期

袁驷，袁全

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

结构动力响应问题是工程计算的重要课题，时程积分方法是最有效的方法之一[1−2]，而自适应步长求解也成为近几年研究的热点之一。就自适应求解而言，多数研究都是对离散解进行误差估计并建立自适应算法[3−4]，但在时域上逐点按最大模控制误差的自适应求解更为困难，也极为少见，而这恰是本文的目标。

采用有限元法进行时程分析，其解是连续的，为最大模控制误差提供了前提条件；为了逐点控制误差，需要有更高精度的超收敛的解。依据有限元数学原理[5]，袁驷等[6−8]提出的单元能量投影(Element Energy Projection，简称EEP)法是一种非常有效的后处理的超收敛算法，可用于进行逐点误差估计和精度修正。该EEP 算法已广泛应用于一维[9−10]、二维[11]，乃至三维有限元中[12]。EEP技术在自适应步长时程积分的有限元分析中也起到了关键的作用[13−15]：基于EEP 技术，可对常规单元构建单步法递推公式[13]，可对有限元解逐点进行误差估计[13−14]，可创建自适应步长算法公式[15]，进而可再次运用EEP 技术对结点位移修正[16]来实现更高效的自适应求解。

为了能够成功实现高效可靠的自适应步长求解，经本文作者课题组持续研究，先后构造了常规单元[13−14]、修正单元[13−14]、凝聚单元[17−18]，乃至本文提出的降阶单元，不断提升时程单元的性能。文本将对这一系列单元做一回顾、分析和评价，从而引出最新的降阶单元。降阶单元是一个采用了非常规算法的常规单元，简单实用，高效可靠，是至今最高性能的单元之一。

为了有助于整体理解，这里对文中内容做几点统一说明：

2) 文中所说结点，如无特殊说明，均指单元端结点，亦即有限元网格结点；

4) 文中以常系数单自由度体系为例展开讨论，但方法完全适用于变系数多自由度问题；

5) 文中所称自适应步长时程单元，是指该单元配备了误差估计器从而可以自适应调整步长的单元；

6) 由于可以逐单元积分求解，不失一般性，仅考虑两端结点坐标为(t1,t2)=(0,h)的典型单元e，其中h为单元长度，亦为时间步长；

7) 表1 给出了文中所述6 种单元的性态比较，标×的和精度比< 2 的部分是该单元的缺点或不足，正因为这些不足，才导致我们构建了高性能的凝聚单元和降阶单元。

表1 各类单元性态比较(次)Table 1 Comparison of various elements of degree

注： m 次降阶单元在单元上的整体精度是m+2阶的，但是在估计误差和控制误差时被降阶为m+1阶的。

方程类别单元编号单元类型无条件稳定 EEP计算结点修正凝聚计算结点精度(收敛阶) 单元精度(收敛阶) 精度比二阶运动方程1常规单元[13]×要无−2m m+1<2 2修正单元[13]×要要−2m+2 m+12 3常规单元[14]√要无−2m m+1<2 4修正单元[14]×要要−2m+2一阶运动方程m+12 5凝聚单元[17]√要无要2m+2 m+12 6降阶单元√无无无2m+2 m+1 m+2,2

8) 总体结论是，作为自适应步长时程单元，表1 中的单元可分为三类：两种常规单元(单元1和单元3)为基本不可用单元；两种修正单元(单元2 和单元4)为基本可用单元；而凝聚单元和降阶单元(单元5 和单元6)为高性能单元。

1 二阶运动方程及其有限元求解

本节介绍求解二阶运动方程的两类有限元：常规单元和修正单元。简言之，虽然修正单元远优于常规单元，是基本可用单元，但由于这两种单元都是有条件稳定的，均不属于高性能单元。

1.1 二阶运动方程

运动方程最常见的表达是如下的二阶常微分方程初值问题：

式中：m、c和k分别为质量、阻尼和刚度；P为外荷载；u为动位移；u0、v0分别为初始位移和初始速度；为时域的上界。

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

则相应的Galerkin 法归结为求解u∈使得：

1.2 常规单元

本文将以上按常规算法构造的单元称为二阶运动方程的“常规单元”[13]，该单元简单易行，但有以下几点缺点(参见表1 的单元1)：

1) 其解uh在结点上是 2阶收敛的，与单元精度比小于2，长时间域问题难以控制住误差；

2) 对于线性元，精度比为1，EEP 解不具有超收敛性，是失效的；

3) 是有条件稳定的，这是一个“硬伤”。

1.3 修正单元

2 一阶运动方程及其有限元求解

鉴于二阶运动方程难以构造无条件稳定的单元，本节介绍求解一阶运动方程的有限元：常规单元和修正单元。简言之，常规单元是无条件稳定的，但是精度比低于2；而修正单元精度比达到期望值2，可惜又变为有条件稳定的。由于这两种单元各有不利和不足，也难属于高性能单元。

2.1 一阶运动方程

现将二阶运动方程等效地转换成一阶方程组(称为一阶运动方程)[14]：

本文将本节中按常规算法构造的单元称为一阶运动方程的“常规单元”[14]，该单元除了可以同时给出位移和速度，也能同时对位移和速度进行误差控制这些优点之外，其最大的优点是：各个次数的常规单元均是无条件稳定的[19]。尽管如此，该单元仍有以下缺憾(参见表1 的单元3)：

1) 其解uh在结点上是2 阶收敛的，与单元精度比小于2，长时间域问题难以控制住误差；

2) 对于线性元，精度比为1，EEP 解不具有超收敛性，是失效的。

2.3 修正单元

记EEP 超收敛解的误差e∗=u−u∗。将其代回式(11)的原问题，可以得到e∗的控制微分方程和初值条件：

本文将按以上做法对结点位移进行修正后的单元称为一阶运动方程的“修正单元”[14]。修正单元大幅提高了结点位移的精度，使其精度比对各次单元均可达到2，特别是也解救了线性元，成为基本可用单元。本以为修正单元在无条件稳定的常规单元的基础上作精度修正，仍可保留其无条件稳定性，但令人遗憾地发现，修正单元改进了常规单元的前两个缺点，但将原本的优点反转为一个缺点(参见表1 的单元4)：不是无条件稳定的。

3 凝聚单元和降阶单元

一阶运动方程的常规单元是无条件稳定的，所欠缺的是精度比不够高；而修正单元提升了精度比，但却丧失了无条件稳定性。本节讨论的两类单元均是在保留常规单元性态的基础上，进一步提高精度比而构造的。首先介绍凝聚单元，然后推出本文提出的降阶单元。简言之，这两类单元都是无条件稳定的，无须额外的结点修正技术即可将精度比直接提升为2，而降阶单元还自带超收敛解，无须EEP 超收敛计算。这两类单元均属高性能单元。

3.1 凝聚单元

3.2 降阶单元

降阶单元不涉及新的公式，仅需用语言表述如下：

图1 例1 的一个单元解(h=3.5)Fig. 1 One element solution of example 1 (h=3.5)

以上一系列算法和操作使得降阶单元成为一个十分简单、简洁、简化的“三无”单元：无凝聚操作，无EEP 计算，无结点修正。凝聚单元和降阶单元的无条件稳定性、各方面的收敛阶等都有严格的数学证明，表现也很类似，但在实施上降阶单元具有很大优势，是一个简单朴实、功能齐全、性能卓越的单元(参见表1 的单元6)。

4 数值算例

表2 有阻尼简谐振动结果 (256 s，三次元)Table 2 Results of damped harmonic motion (256 s,Cubic)

图2 例1 步长分布图Fig. 2 Step-size distribution for example 1

图3 例1 单元位移误差比图Fig. 3 Error ratios of element displacements for example 1

例2. 多自由度简谐振动

物理模型来源于3 层剪切型框架结构，计算数据为：

简谐荷载向量为F=10(sin(10t) sin(10t)sin(10t))T，初始位移和初始速度均取0，取10 s，初始步长取为h0=0.5。表3 为三次元的相关计算结果，由其可知，对多自由度本法同样有效，其性质与单自由度相似，步长分布均匀且合理。图4、图5 为三次元、tol=10−3的数据结果。由表3 和图4 可见，由于其多自由度特性，尽管自适应次数和迭代次数仍然较少，但是多于单自由度问题。由图5 可知，误差分布多在误差上下限区间内。

图5 例2 单元误差比图Fig. 5 Error ratios of elements for example 2

表3 多自由度简谐运动计算结果 (10 s,三次元)Table 3 Results of multi-degree-of-freedom system motion (10 s,Cubic)

图4 例2 步长分布图Fig. 4 Step-size distribution for example 2

例3. 与凝聚单元比较算例

为了对降阶单元和凝聚单元的特点作一比较，用凝聚单元计算了上述两个算例，将结果比较列于表4 和表5。可见，两种单元均能够圆满完成自适应求解。总体上，相同次数的降阶单元所需要的单元数和变单元长度次数均比凝聚单元少。二者的差别主要是采用了不同的误差估计器，凝聚单元采用EEP 解估计误差u∗−uh，而降阶单元采用(t)估计误差。整体上看，降阶单元的表现稍优于凝聚单元。