时空有限元自适应分析需要在空间和时间两方向进行有限元离散，须满足不同方向的特定要求进行自适应分析。譬如，时间方向的单元应具备无条件稳定性，而空间方向的单元需满足相应的连续性要求。对于自适应分析，则进一步要求空间和时间有限元都具备可靠的误差估计器。因此，时空有限元须应对空间有限元和时间有限元两方面的挑战。

笔者研究团队近年来先后提出了初值问题(时域问题)^{[1 − 2]}和边值问题(空间问题)^{[3 − 4]}降阶单元自适应分析算法，经初具规模的研究表明：将时域问题和空间问题的降阶单元结合为一体的时空降阶单元，可以有效应对时、空两方面的挑战。

降阶单元自适应分析的基本作法非常简单，可统一简述如下：对于$m$次降阶单元，采用$m + 1$次常规单元求解，得到常规有限元解答(称为满阶解)后，略掉其中最高次($m + 1$次)项，得到低一次($m$次)的降阶解，作为降阶单元的解答；转而将略掉的最高次项作为降阶解的误差估计器。图1和图2分别给出了一维和二维降阶线性元作法的示意，即，降阶单元将常规单元的解答分解(${u^h} = u_{\text{R}}^h + \varepsilon _{\text{R}}^h$)为降阶解$u_{\text{R}}^h$和误差项$\varepsilon _{\text{R}}^h$，利用满阶解和降阶解在精度上的“阶差”，很自然地提供了一个可作逐点估计的误差估计器，也自然地发展出以降阶单元作为最终解答的自适应算法。

图  1  一维降阶单元示意图

Figure  1.  Schematic diagram of one-dimensional reduced element

下载: 全尺寸图片 幻灯片

图  2  二维降阶单元示意图

Figure  2.  Schematic diagram of two-dimensional reduced element

下载: 全尺寸图片 幻灯片

降阶单元的概念和作法首先由文献[1 − 2]在初值问题中提出，思路上是对初值问题中凝聚单元的逆构思和逆运用^{[5 − 7]}，现已成为结构动力计算中的一种新型高效、可按最大模自适应步长的时程单元。在此基础上，降阶单元进一步推广到边值问题的自适应有限元分析，对于一维边值问题、二维Poisson方程和弹性力学平面问题等均做了成功的自适应分析求解^{[3 − 4]}。降阶单元具有很多优点，简单概括几点如下：

① 不限单元次数，可统一构造出$m\;( > 0)$次降阶单元；② 不限单元形状，二维问题可以统一构造出四边形和三角形单元；③ 不受网格约束，可用于具有局部加密的非规则网格划分；④ 对于初值问题是无条件稳定的，适用各类初值问题；⑤ 一维问题的结点收敛阶是单元内收敛阶的2倍，可用于长时域问题； ⑥ 结点收敛阶比单元内收敛阶至少高1阶，使误差局限于单元内部，大幅增强局部加密的适用性和可靠性；⑦ 不限问题类型，通用于线性和非线性问题；⑧ 内置了最大模误差估计器，可以给出按最大模满足误差限的自适应有限元解答。

本文以弦振动方程(亦为波动方程，弹性杆轴向振动方程)为模型问题展开讨论，对相应的时空降阶单元及自适应分析算法研究的最新进展做一个综述性介绍。为展示本法的可行性和有效性，文中给出多个典型数值算例，包括初始速度、初始位移、刚度突变、有无阻尼、集中荷载、均布荷载乃至移动荷载等问题，均可获得满意结果。

1   问题描述与有限元求解

1.1   模型问题

考虑如下的拓展的一维弦振动方程(波动方程，wave equation)及其定解条件：

式中，系数$\overline m$、$c$、$p$、$q$均为$x$的函数。式(1)对应着多个物理模型：若$\overline m = p = 1$、$c = q = 0$，它是经典的波动方程；同时它也是经典的(拓展的)弦振动方程，其中各个系数依次代表质量、阻尼、张力刚度和弹性地基刚度；同时它还是弹性杆轴向强迫振动方程，其中各个系数依次为质量、阻尼、轴向刚度和弹性地基刚度。本文取弹性弦振动作为物理模型，为了与单元次数符号$m$区分，质量采用了符号$\overline m$表示。

1.2   时空区域半离散

首先对空间方向用一维$m$次单元进行离散，相当于用有限元线法(finite element method of lines，FEMOL，简称线法)^[8]单元对时空域进行离散，结线取在时间方向，结线位移$\{ d\}$保留为连续的时间t的函数。由于含有阻尼项，不便直接采用变分法求线法解，采用Galerkin法推导出以结线位移为未知量的二阶常微分方程组(二阶运动方程)。

定义与式(1)相应的双线性型和线性型为：

采用Galerkin法求解式(1)，则其近似解可归结为求解$u \in H_E^1$，使得：

式中：$u$、$v$分别为试探函数和检验函数；$H_E^1$为满足本质边界条件($u(0,t) = 0,\;\;\;v(0,t) = 0$)直到一阶导数均平方可积的函数空间。

时空单元上的位移$u_{}^h(\xi ,\;t)$的半离散采用FEMOL插值形式，其试探函数和检验函数采用Hierarchical(升阶谱)形函数，形式如下：

其中：

式中，${\hat N_i}(\xi )$为单元内部Hierarchical形函数。将以上位移模式代入到式(2)中，可以推导出FEMOL单元系数矩阵和向量如下：

所有空间单元代入到式(3)中，集成后可得整体系数矩阵和向量(去掉式(6)中的上标e)。再由${\{ v\} ^h}$的任意性，有FEMOL方程：

1.3   一阶时程单元

对于线法导出的式(7)中的二阶常微分方程组(二阶运动方程)，为获得无条件稳定算法^{[1, 7]}，将其等效地转化为一阶常微分方程组(一阶运动方程)：

式(8)可表示为如下标准的一阶初值问题：

即，结线上的位移$\{ d\}$和速度$\{ v\}$整合为广义结线位移$\{ {\boldsymbol{u}}\}$。对于时间方向，结线位移$\{ {\boldsymbol{u}}\,(\;t)\}$采用相同的$m$次Hierarchical(升阶谱)形函数离散：

其中：

且$[{\boldsymbol{I}}]$为$(m + 1) \times (m + 1)$阶单位矩阵。至此，双m次降阶时空单元得以构造。

采用Galerkin法可推导出单步法的时程单元的递推矩阵方程组，可逐步长(逐时程单元)自适应积分求解^[1]。

1.4   降阶解的提取

为了区分，降阶解用下标${\text{R}}$(代表Reduced)加以标识。由于采用了升阶谱形函数，降阶解的提取十分方便，只需将最高次($m$次)的结点(结线)位移设成0即可，即令$u_m^h = 0$和$\{ {\boldsymbol{u}}\} _m^h = \{ 0\}$。编程时，可直接对相应矩阵和向量中相应分量赋0值即可：

2   误差估计及其他相关算法

2.1   误差估计

记空间方向单元长度为${h_e}$，时间方向单元长度(时间步长)为${h_t}$。降阶单元用于自适应求解时，取降阶解$u_{\text{R}}^h$为最终的解，将满阶解$u_{}^h$看作“精确解”，则可方便地获得时程单元${h_t}$上结线位移和二维时空单元(${h_e} \times {h_t}$)上逐点的误差估计：

2.2   初值的处理

对于给定的初始位移${\overline u_{ 0}}(x)$和初始速度${\overline v_{ 0}}(x)$，采用$m$次分片多项式自适应插值获得其逼近解$\overline u_{\,0}^h(x)$、$\overline v_{\,0}^h(x)$，使其按最大模度量满足用户给定的误差限：

插值完成后，可得到初始空间网格${\pi _0}$。

2.3   悬挂单元的处理

对于每个时间步长(单元)${h_t}$，空间方向上都有相应的单元网格划分，本文将其称为“单元层”。由于各个单元层的单元网格划分允许不同，因此会出现悬挂单元。参见图3，第2层上单元$({\tilde e_{ 1}})$、$({\tilde e_{ 2}})$为悬挂单元，悬挂在第1层的承担单元$(e)$上。

图  3  悬挂单元示例

Figure  3.  Example of hanging elements

下载: 全尺寸图片 幻灯片

暂将${u^h}$简记为$u$。又记搭接边上三个单元的位移为：

此时，悬挂单元结点解${\boldsymbol{d}}_{{S_1}}^{\,({{\tilde e}_{\,1}})}$和${\boldsymbol{d}}_{{S_1}}^{\,({{\tilde e}_{\,2}})}$并非完全独立，应该用承担单元的结点解${\boldsymbol{\overline d}}_{{S_2}}^{\,(e)}$的线性组合表示。为处理悬挂结点，引入2个转换矩阵(${\boldsymbol{z}}$矩阵)：左半边的${{\boldsymbol{z}}_1}\,$，右半边的${{\boldsymbol{z}}_2}\,$。则通过插值条件可得转换关系：

式中，${{\boldsymbol{z}}_1}\,$和${{\boldsymbol{z}}_2}\,$矩阵可有显式表达式。

对于悬挂单元，本文算法不允许出现多于2个单元悬挂在同一个大单元上。若发现此种情况，则实施回溯加密：即退回到前一步长的单元层，将“一对多”的大单元在空间方向二分，在新网格上继续求解。

2.4   移动荷载的处理

移动集中荷载的强迫振动反应是一个典型的问题。设移动荷载${P_0}$以速度$\beta$从左至右经过弦，时刻t时的作用点为${x_c} = \beta \,t$。图4所示为移动荷载经过某单元的示例。由于集中荷载在数学上为Dirac函数$\delta (x - {x_c})$，其积分为单位阶跃函数$H(x - {x_c})$，经分析推导，式(6)中的等效结线荷载向量可用单位阶跃函数表示为：

其中，${t_{c1}} = {{{x_1}} \mathord{\left/ {\vphantom {{{x_1}} \beta }} \right. } \beta }$，${t_{c2}} = {{{x_2}} \mathord{\left/ {\vphantom {{{x_2}} \beta }} \right. } \beta }$。

图  4  移动荷载经过单元示意

Figure  4.  Example of moving load on an element

下载: 全尺寸图片 幻灯片

3   自适应算法

3.1   自适应求解目标及算法

如前所述，本文的降阶单元自适应分析算法，将降阶解$u_{\text{R}}^h$取为最终解，用满阶解$u_{}^h$代替精确解来估计$u_{\text{R}}^h$的误差。本文的自适应求解的终极目标为：由用户事先给定误差限$tol$，寻找一个优化的有限元网格，使得该网格上的降阶单元解$u_{\text{R}}^h$按照最大模度量满足误差限，即逐单元满足：

由于精确解$u$一般事先未知，在实施中，采用以下停机准则：

为了满足式(19)，对于时程单元的各个结线解需要事先满足：

由于速度$\dot d_i^h$常会出现高梯度乃至跳跃情况，加之$u_{}^h$的计算式(4)中并没用到结线上的速度，通常可放松对$\dot d_{{\text{R,}}i}^h$的误差控制。大量数值算例表明，对于绝大多数所求解的问题，随着误差限的缩小，满足式(19)和式(20)的解基本上都可以满足式(18)。

3.2   自适应求解算法

至此，时空降阶单元自适应求解算法可简略地归纳为：

1) 给定误差限$tol$、降阶单元次数${m_{\text{R}}}$，对${\overline u_0}$和${\overline v_0}$做自适应插值，得到初始空间网格${\pi _0}$；

2) 生成FEMOL的矩阵和向量$[A],\;[B],\;[C],\;\{ F\}$ (参见式(6))，形成一阶运动方程式(9)；

3) 时程降阶单元自适应步长求解式(9)^[1]直到各个结线上的解满足式(20)，得步长${h_t}$，及结线上的满阶解${\{ {\boldsymbol{u}}\} ^h}$和降阶解$\{ {\boldsymbol{u}}\} _{\text{R}}^h$；

4) 对该单元层各个单元的降阶解检验误差是否满足式(19)；

5) 若满足，则空间网格不动，回到步骤3)，执行下一步时间步长；

6) 若未满足，则将该单元在x方向二分，调整设定初始解解后，回到步骤2)；

7) 直至步长到达时域终点。

4   数值算例

为方便起见，引入“误差比”，记作$\overline e_{\text{R}}^h = {(u - u_{\text{R}}^h)} / {tol}$，最大模的误差比记为$\overline e_{\max }^h = {{\max | {u - u_{\text{R}}^h} |} / {tol}}$，而估计的误差比记作$\overline e_{\max }^{*h} = {{\max | {{u^h} - u_{\text{R}}^h} |}/ {tol}}$。又记${N_{{\text{adp}}}}$为自适应步数，${N_h}$为时间步长(单元)数，${N_e}$为时空区域总体单元数。借鉴通常的提法^{[9 − 10]}，误差估计可有20%的容差，折算为误差比，则真实误差比$\overline e_{\max }^h {\leqslant} 1.25$，均算成功。对于本文方法的算例，随着误差限减小，网格加密，多数都能严格满足式(18)，即$\overline e_{\max }^h {\leqslant} 1$。

熟知，误差估计的准确与否依赖于单元的尺度$h$和真解高阶导数的模^{[11 − 12]}，而弦振动方程的特性是几乎所有问题都存在不同阶数的奇异性：二阶导数都是不连续的，甚至$u$也可以是不连续的(移动荷载问题)。尽管这一特性会导致奇异区域附近的误差估计失准，文中绝大多数算例，即使在很宽松的误差限情况下，也能将误差控制在$\overline e_{\max }^h < 2.5$限度内。

本节的算例均可以构造出三角级数解^{[13 − 15]}，为了检验算法，部分算例以10项($n = 10$)级数解为求解目标，荷载$f$用级数解由原方程导出；而由于级数解引入了振荡性、弱化了奇异性，往往并不能真实反映原问题的特性，部分算例也采用原问题定解条件直接求解。所有算例都仅对位移u控制误差。

例1. 初始常速度问题(级数解为目标解)

考虑一个两端固定的弦，受到一个初始速度$\overline v_0^{}(x) = 1$的激励，问题定义如下：

由于弦在两端固定，速度$v$在两端始终为0，即：$v(0,t) = v(1,t) = 0$。因此，初始速度$\overline v_0^{}$在两端点有跳跃。

取10项级数解作为自适应求解目标，$\overline v_0^{}$如图5(c)所示。取不同单元次数m、不同误差限tol计算，时长$T = 1\;{\text{s}}$，结果见表1。图6给出了$tol = 0.0001$时不同次数降阶单元最终的网格划分。可见，网格疏密分布合理；尽管速度$v$有突变，精度仍得到很好保证，多数算例均满足误差限，高次元的表现更为出色。

图  5  10项级数解的位移u、速度v和初始速度$\overline v_0$

Figure  5.  Displacement u, velocity v, and initial velocity $\overline v_0$ from a 10 termed series solution

下载: 全尺寸图片 幻灯片

表  1  初始常速度问题 (级数解为目标解)

Table  1.  Problem of constant initial velocity (the series solution is the target solution)

$tol = 0.001$							$tol = 0.0001$
降阶单元 次数${m_{\text{R} }}$	自适应 迭代次数${N_{ {\text{adap} }} }$	空间 单元数${N_h}$	时间 单元数${N_e}$	估计 误差比$\overline e_{\max }^{*h}$	真实 误差比$\overline e_{\max }^h$		自适应迭代 次数${N_{ {\text{adap} }} }$	空间 单元数${N_h}$	时间 单元数${N_e}$	估计 误差比$\overline e_{\max }^{*h}$	真实 误差比$\overline e_{\max }^h$
2	32	18	606	0.976	1.791		61	46	2672	0.904	0.952
3	17	14	240	0.978	1.155		37	26	816	0.986	1.012
4	13	10	132	0.998	1.067		24	15	344	0.977	1.019

下载: 导出CSV 

| 显示表格

图  6  不同次数降阶单元求解的最终网格

Figure  6.  Final meshes of the adaptive analysis with reduced elements of deferent degrees

下载: 全尺寸图片 幻灯片

例2. 初始常速度问题(折线逼近常速度)

问题同例1，放弃级数解，荷载$f = 0$。对初始速度采用几乎为常数1的折线逼近(图7(c))：

图  7  降阶2次元求解的位移u、速度v和初始速度$\overline v_0$

Figure  7.  Displacement u, velocity v, and initial velocity $\overline v_0$ from reduced elements of degree 2

下载: 全尺寸图片 幻灯片

设误差限$tol = 0.0001$，仍取不同单元次数m计算。图7给出了降阶2次元自适应求解后的位移$u_{\mathrm{R}}^h$、速度$v_{\mathrm{R}}^h$和采用的初始速度$\overline v_0$。图8给出了不同次数降阶单元最终的网格划分。可见，网格疏密分布合理，特别是位移$u$的斜坡面和速度$v$的平台面的解答都更加接近平面，比级数解更符合原题的真解。

图  8  例2的不同次数降阶单元求解的最终网格

Figure  8.  Final meshes of the adaptive analysis with reduced elements of deferent degrees of Example 2

下载: 全尺寸图片 幻灯片

例3. 移动集中荷载(级数解为目标解)

考虑一个两端固定的弦，移动荷载${P_0}$以速度$\beta$从左至右经过弦。问题的数学定义如下：

简单取${P_0} = 1$，$\beta = 1$。图9给出了$\beta = 1$时位移$u$的精确解，可见位移$u$沿$x = t$线两侧是不连续的，导数是奇异的。取10项级数解作为自适应求解目标，由方程导出荷载$f$，如图9(c)所示。取不同单元次数m、不同误差限tol计算，时长$T = 1\;{\text{s}}$，结果如表2所示。可见，随着误差限的缩小，误差得到更好的控制，多数算例均严格满足误差限，高次元的表现更为出色。图10给出了$tol = 0.01$时不同次数降阶单元最终的网格划分。可见，网格疏密分布合理，算法具有较好的自适应能力。

图  9  移动荷载问题位移的精确解u、10项级数解及级数解导出的荷载

Figure  9.  Exact solution u, series solution of ten terms and its derived load term

下载: 全尺寸图片 幻灯片

表  2  移动荷载 (级数解为目标解)

Table  2.  Problem of moving load (the series solution is the target solution)

$tol = 0.01$							$tol = 0.001$
降阶单元 次数${m_{\text{R} }}$	自适应迭代 次数${N_{ {\text{adap} }} }$	空间 单元数${N_h}$	时间 单元数${N_e}$	估计 误差比$\overline e_{\max }^{*h}$	真实 误差比$\overline e_{\max }^h$		自适应迭代 次数${N_{ {\text{adap} }} }$	空间 单元数${N_h}$	时间 单元数${N_e}$	估计 误差比$\overline e_{\max }^{*h}$	真实 误差比$\overline e_{\max }^h$
2	45	9	136	0.981	2.302		58	21	695	0.999	2.057
3	31	7	89	0.989	1.319		43	13	297	0.997	0.983
4	9	5	28	0.881	0.951		22	8	100	0.905	0.938

下载: 导出CSV 

| 显示表格

图  10  例2的不同次数降阶单元求解的最终网格

Figure  10.  Final meshes of the adaptive analysis with reduced elements of deferent degrees of Example 2

下载: 全尺寸图片 幻灯片

例4. 移动集中荷载(用集中荷载求解)

问题同例3，放弃级数解，直接采用集中荷载，按照式(17)计算荷载向量。误差限$tol = 0.01$，图11给出了降阶2次元自适应求解后的位移、误差分布和最终的单元网格。可见，网格疏密分布合理，特别是位移$u$的斜平面和水平面的解答，比波纹式的级数解更符合原题的真解，误差主要集中在不连续的对角线上。

图  11  降阶2次元求解的位移、误差和最终网格

Figure  11.  Displacement, error, and final mesh from the solution of reduced elements of degree 2

下载: 全尺寸图片 幻灯片

例5. 刚性问题

此例为文献[16]中用其时间、空间混合有限元计算的问题，其数学定义如下：

若将其物理模型看作两端固定的弦，右半跨的刚度是左半跨的400倍，几乎是刚性约束。文献[16]采用了一种混合单元计算，空间方向取了20个均匀划分的单元，时间步长在左半跨取为0.05、右半跨取为0.002；简单而言，其网格是“右密左疏”。用本文方法计算此题，采用降阶2次元(${m_{\text{R}}} = 2$)，$tol = 0.025$，$T = 1\;{\text{s}}$。图12给出了降阶单元自适应后的解答和网格划分，可见按照自适应分析，时间步长大致为0.004，而网格应为“左密右疏”，而且左半侧从下三角稀疏网格随时间逐步加密也符合受力特性。本例从一个侧面体现了自适应求解的必要性，同时也展示了本文方法的良好自适应性能力。

图  12  降阶2次元求解刚性问题的位移和网格划分

Figure  12.  Solution of the stiff problem and the final mesh from reduced elements of degree 2

下载: 全尺寸图片 幻灯片

例6. 有阻尼突加均布荷载问题

考虑一个两端固定的弦，无初位移和初速度，赋予阻尼$c = 1$，突加均布荷载$f = 1$，问题定义如下：

用降阶2次元计算，时长$T = 15\;{\text{s}}$，$tol = 0.001$。图13给出了自适应求解的位移、速度。作为检验，在时域终点$t = 15\;{\text{s}}$，弦基本处于稳态，其计算的跨中位移为0.125 029，与稳态的理论解${1 /8}$十分相符。

图  13  二次降阶单元解的位移、速度和稳态位移

Figure  13.  Displacement, velocity and steady state displacement from reduced elements of degree 2

下载: 全尺寸图片 幻灯片

5   结论

本文集初值问题和边值问题降阶单元的优势于一体，提出时空降阶单元，以弹性弦振动方程(波动方程)为例，构造了自适应有限元分析算法，可以给出按最大模度量满足误差限的降阶单元的解答。本文给出若干强迫振动以及移动荷载反应的典型算例以展示本法的可行性和有效性，以及广阔的发展前景。

注：该文在第33届结构工程学术会议 (2024 阜新) 应邀作特邀报告