温彦良,李 彬,唐小微

(1.辽宁科技大学 矿业工程学院，辽宁 鞍山 114051；2.泰山学院 机械与建筑工程学院，山东 泰安 271000；3.大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024)

如何有效地提高计算效率，一直是一件非常重要和有意义的事情，众多学者在这方面的努力一直没有停止过。例如，王林等[2]研究基于改进MP的稀疏表示快速算法，该算法能准确提取滚动轴承故障特征且提高了计算效率；为减少谐波合成法中功率谱矩阵分解的计算量,祝志文等[3]提出了谱解矩阵双轴插值算法和递归插值算法，两种方法均使风场模拟计算效率大幅度提高。

在众多提高计算效率的方法中，时间自适应分析方法受到越来越多的关注。它将误差控制在某一允许范围内，同时尽可能地减少计算时间，这样就在保证计算精度的同时节省了计算时间，最终提高了计算效率。在时间自适应研究领域，国内外很多学者做出了不懈的努力。例如，Zienkiewicz等[4]提出了一种简单的时间自适应方法，并成功应用到动力分析中。Zeng等[5-6]对Zienkiewicz等提出的方法进行了改进。Kavetski等[7]采用了启发式和误差评估的方法进行了时间自适应。Kuo等[8]结合时间元素具有大时间步幅及动量平衡具对不连续载重的平滑化的优点，提出了一套加权式动量时间元素的逐步时间积分方法。王开加等[9]在海上浮基风电平台绕流数值模拟分析中，采用了时间自适应技术。毛卫男等[10]提出了一种适应时间的方法，并应用于土壤冻结的水热耦合模型中。Li等[11]提出了一套时间步长自动优化方法，并验证了该方法在提高计算效率方面的有效性。Naveed等[12]将高阶变分离散时间的自适应时间控制应用到了对流-扩散反应方程中。Vahid等[13]应用时间适应方法，在裂纹扩展的相场公式中进行了大规模的并行处理。正如Marc等[14]提到的，目前时间自适应方法是后验式误差评估时间自适应方法。不同于现有的时间自适应方法，ATAM是先验式研究领域的一次有益尝试。

ATAM在应用时，我们不可能为了一个具体的算例重新编写程序，那样就失去了提高计算效率的意义。目前已经有很多成熟的数值计算程序，将ATAM嵌入到现有程序中，提高原程序的计算效率和计算稳定性是我们开发这个算法的最终目的。

在《与荷载同步变化的时间步自动调整方法》公式的推导过程中，选择了Newmark-β法的解做为近似解。当采用中心差分法的解做为近似解时[15]，会得到与前者相似的结果，但略有不同。本文将比较两者的不同，从而得出ATAM在实际应用中需要注意的问题。

1 自适应时间步长确定方法

结构动力数值分析方法中，动力控制微分方程通常具有如下形式

(1)

由中心差分法，上述动力微分方程tn+1时刻位移近似解ui+1(t)可以表达为

(2)

(3)

(4)

计算相对误差时需要精确解，然而精确解通常是难以获得的，这里采用泰勒级数的展开法获得ti+1时刻位移的解来替代精确解进行误差评估，其表达式为

(5)

(6)

将式(2)、(5)和(6)代入式(4)中，整理后得

(7)

(8)

由于ATAM只调整时间步长，不改变原程序的计算方法，所以自适应后程序的收敛完全依赖于原计算方法。

当采用Newmark-β法的解做为近似解时

(9)

两者关系为

(10)

两者在时间步长确定公式上的细小差别，在实际应用中却代表着两种不同的类型。类型A：原程序在时域离散时采用的方法，与自适应时间步长推导过程中采用近似解的方法完全相同，包括相关的参数设置也相同。例如在《与荷载同步变化的时间步自动调整方法》中，原程序采用Newmark-β法对时域进行了离散，自适应时间步的确定采用了Newmark-β法解做为近似解，且控制参数相同β=0.25。类型B：方法不相同。例如在自适应时间步长的确定中采用了中心差分解为近似解，而原程序采用Newmark-β法。在实际应用中，我们遇到的情况更多是后者。为了方便比较，本文选用了与《与荷载同步变化的时间步自动调整方法》中同一算例的同一结点进行了分析。

2 自适应分析

Bernoulli-Euler梁：一个等截面均质的简支梁，初始位移与初始速度均为0，有一个集中力荷载从梁最左端向右端匀速运动，其数值模型如图1所示。

移动荷载下Bernoulli-Euler梁运动微分方程[16]

(11)

图1 移动荷载下简支梁及节点分布Fig.1 Simply supported beam under moving load and node distribution

(12)

将梁平均分为20个计算单元，节点分布如图2所示，各节点的初始位移和初始速度均为0，有一个集中力荷载由梁最左端向右端匀速运动。计算参数均假设为无量纲参数(其中L、l分别为梁长和单元长)

EI=100,c=0,L=10,l=0.5,

(13)

2.1 以中心差分解为近似解的自适应分析

选取梁的11号节点，采用固定时间步长为0.02的Newmark-β法，在固定荷载从3号节点匀速移动到19号节点的时间段内，给出11号节点挠度随时间变化的图形，并与解析解进行了对比，结果见图2。

相对误差随时间变化的结果，见图3。相对误差在0～0.005 6的范围内变动，取得最大时所对应的计算精度就是Newmark-β法在固定步长为0.02时所达到的计算精度。

对比自适应前后的两个局部放大图(图2(b)与图5(b))，可以看出自适应前，近似解等时间步间距不规则的分布在解析解附近。自适应后，近似解与解析解的相对误差基本保持一致，时间步间距不相同。自适应后的点明显少于自适应前的点，而点的多少代表计算次数的大小，点越少计算次数越小，计算时间就越少。

相比自适应前，自适应后的近似解明显偏离解析解，但是相对误差基本保持一致。这是因为相对误差与自适应前整个计算过程中的最大相对误差一致造成的，因此，自适应前后的计算精度保持一致。由此可见，自适应过程提高原程序计算效率的效果是明显的。

(a)时间步为0.02

(b)局部放大

图3 Newmark-β法的相对误差Fig.3 Relative error of Newmark-β method

图4 自适应后的相对误差Fig.4 Relative error after time adaptive

(a)自适应时间步

(b)局部放大

2.2 两种类型的比较与分析

2.2.1 相同点

在提高计算效率方面，两者都可以显著的提高原程序的计算效率；在荷载与自适应时间步长的关系上，两者也是一致的，可以解决在冲击荷载、振动荷载、爆炸荷载或地震荷载的动力数值计算中，因时间步长的设置不当导致计算结果发散甚至计算中断的问题，在一定程度上提高了原程序的计算稳定性，相关内容可参考《与荷载同步变化的时间步自动调整方法》。下面就两种类型在应用时的另一个共同优点做一个简单介绍。

针对一个具体的数值计算方法，在实现时间自适应后，计算效率具体提高了多少？目前的后验式时间自适应方法无法解决这个问题，因为在没有解析解的情况下，无法获得准确的相对误差，计算精度就无法准确比较，从而计算效率的提高程度也就无法获得。ATAM可以解决这个问题。

在图4的自适应过程中，最大自适应步长为0.394 9，最小自适应步长0.020 9(与0.02基本一致)。以最小步长做为自适应前原程序的固定时间步长，相对误差变化如图6所示，最大相对误差0.006 1。

图6 Newmark-β法的相对误差Fig.6 Relative error of Newmark-β method

根据前面得出的结果，可从两方面验证该方法的可行性。①最小自适应时间步长为0.020 9与固定时间步长0.02基本一致；②图6的最大相对误差0.006 1与图4的0.005 5基本一致。这两点都可说明自适应前后的计算精度基本保持一致，在此基础上比较两者的计算时间，就可实现自适应前后计算精度提高程度的比较。

自适应的过程中，最大时间步长与最小步长之间跨度非常大。在这里最大是最小的18.89倍，如果是冲击荷载、振动荷载、爆炸荷载或是地震荷载，相信这个跨度会大很多。鉴于此，在程序的实现过程中，应对最大自适应步长加以限制，由于ATAM的推导过程中对时间步长要求Δt<<1，所以建议最大时间步长不超过0.4，所有超过0.4的时间步长都按0.4计算。此建议不影响自适应前后计算效率的比较。

2.2.2 不同点

首先，两者代表应用时的两种类型A、B，这在第一部分的自适应时间步长确定方法中已经提到了。由于中心差分法自身的缺陷，很少有数值计算程序选用它来进行时域离散，所以在实际应用中建议采用Newmark-β法的解做为近似解的式(9)来确定自适应时间步长，当遇到原程序采用Newmark-β法对时域进行离散时按A类处理，其他情况按B类处理。

图7 自适应后的相对误差Fig.7 Relative error after time adaptive

在本文与《与荷载同步变化的时间步自动调整方法》中，类型A与B之间存在如下的对应关系

(14)

即：类型B设定的相对误差限是自适应后实际得到的1.333倍。如果ATAM选取Newmark-β法的解做为近似解的式(9)来确定自适应时间步长，而原程序中时域离散的方法是其他方法，这样的对应关系就难以确定，具体关系需要具体问题具体分析，所以这样的对应关系只限于此。

3 荷载与自适应时间步长的关系

关于两者的关系，在《与荷载同步变化的时间步自动调整方法》中已经得出了相关结论，这里补充下由两者的关系得出的另一个适用条件。

根据《与荷载同步变化的时间步自动调整方法》中的结论，可以得到自适应时间步长随时间的变化，如图8所示。

将图8中的时域扩展，得到结果如图9所示。

如图9所示，在位移较小的时候，荷载对时间步长的影响占主导地位，随着位移的增大，位移对时间步长的影响比重会逐渐增大。位移对时间步长有影响是因为在公式的推导中选取了相对误差做为衡量标准，在相对误差保持不变的情况下，位移越大则近似解偏离解析解越远，则时间步长越大。

图8 时间步长随时间的变化Fig.8 Time step change over time

图9 时间步长随时间的变化Fig.9 Time step change over time

由于ATAM的推导过程中要求时间步长Δt<<1，根据前面的结论，位移越大则时间步长越大，所以ATAM在处理大变形问题时，会受到自适应时间步长最大不超过0.4的限制，使得提高计算效率的程度受到影响。

4 结 论

在实际应用中，ATAM分为两类情况：类型A与类型B。由于中心差分法自身的缺陷，很少有数值计算程序选用它来进行时域离散，所以在实际应用中建议采用Newmark-β法的解做为近似解的式(9)来确定自适应时间步长，当遇到原程序采用Newmark-β法对时域进行离散时按A类处理，其他情况按B类处理。

两种类型在计算效率的提高、荷载与自适应时间步长的关系上得出的结论是一致的。两种类型都可以对自适应前后计算效率的提高程度进行准确比较：以自适应过程中的最小时间步长做为自适应前原程序的固定时间步长，这样就保证了自适应前后计算精度保持一致，在此基础上比较两者的计算时间就可以获得计算效率的提高程度。

类型A可以设定预期获得的计算精度，且计算结果与预期相符；类型B不可以。本文中，由于类型A、B存在具体的数值对应关系，所以类型B设定的相对误差限是自适应后实际得的1.333倍，但这样的对应关系只限于此。