结构动力学：理论及其在地震工程中的应用

1、5章 动力反应的数值计算 如果激励作用力或地面加速度是随时间任意变化的，或者体系是非线性的，那么对单自由度体系的运动方程进行解析求解通常是不可能的。这类问题可以通过数值时间步进法对微分方程进行积分来处理。在应用力学广阔的学科领域中，有关各种类型微分方程数值求解方法的文献（包括几部著作中的主要章节）浩如烟海，这些文献包括这些方法的数学进展以及它们的精度、收敛性、稳定性和计算机实现等问题。 然而，本章仅对在单自由度体系动力反应分析中特别有用的很少几种方法进行简要介绍，这些介绍仅提供这些方法的基本概念和计算算法。尽管这些对许多实际问题和应用研究已经足够了，但是读者应该明白，有关这个主题存在大量的知识

2、。 5.1 时间步进法 对于一个非弹性体系，欲采用数值求解的运动方程为 或者 (5.1.1)初始条件 假定体系具有线性粘滞阻尼，不过，也可以考虑其他形式的阻尼（包括非线性阻尼），后面会明显看到这一点。然而由于缺乏阻尼信息因此很少这样做，特别是在大振幅运动时。作用力由一系列离散值给出: , 到N。时间间隔 （5.1.2）图5.1.1 时间步进法的记号通常取为常数，尽管这不是必需的。在离散时刻（表示为时刻）确定反应，单自由度体系的位移、速度和加速度分别为、和。假定这些值是已知的，它们在时刻满足方程 （5.1.3）式中，是时刻的抗力，对于线弹性体系，但是如果体系是非弹性的，那么它会依赖于时刻以前的位

3、移时程和速度。将要介绍的数值方法将使我们能够确定+1时刻满足方程(5.1.1)的反应、和，即在+1时刻 （5.1.4）对于=0，1，2，3，连续使用时间步进法，即可给出=0，l，2，3，所有瞬时所需的反应。已知的初始条件和提供了起动该方法的必要信息。 从时刻到+1时刻的步进一般不是精确的方法，许多在数值上可以实现的近似方法是可能的。对于数值方法，有三个重要的要求：(1)收敛性一随着时间步长的减少，数值解应逼近精确解；(2)稳定性一在存在数值舍入误差的情况下，数值解应是稳定的；(3)精度一数值方法应提供与精确解足够接近的结果。这些重要的问题在本书中均作简要的讨论，全面的论述可在着重微分方程数值解

4、法的书中找到。 本章介绍三种类型的时间步进法：(1)基于激励函数插值的方法；(2)基于速度和加速度有限差分表达的方法；(3)基于假设加速度变化的方法。前两类中各只介绍一种方法，第三类中介绍两种方法。 5.2 基于激励插值的方法对于线性体系，通过在每个时间间隔里对激励进行插值，并利用第4章的方法进行精确求解，能推导出一种非常有效的数值方法。如果时间间隔较短，则线性插值是令人满意的。图5. 2.1所示的时间间隔，激励函数为 （5.2.1a）其中 （5.2.1b）时间变量从0到变化。为数学上简单起见，我们首先考虑无阻尼体系，后面再将该方法扩展到有阻尼体系。待求解的方程为 （5.2.2）在时间间隔内，

5、反应为三部分之和：(1) =0时刻的初位移和初速度引起的自由振动；(2)零初始条件下对阶跃力的反应；(3)零初始条件下对斜坡力的反应。对这三种情况分别采用来自2.1、4.3和4.4中已有的解答，得 （5.2.3a） （5.2.3b）计算时的这些等式，得+1时刻的位移和速度： (5.2.4a) (5.2.4b)将式（5.2.1b）代入后，可将这些等式重写为如下的递推公式： (5.2.5a) (5.2.5b)对于欠临界阻尼体系（即），重复上面的推导，表明式(5.2.5)也适用于有阻尼体系，系数A，B，,的表达式由表5.2.1给出，这些系数取决于体系的参数、k和以及时间间隔。 表5.2.1 递推公式

6、中的系数（）因为递推公式是从运动方程的精确解推导出的，因此对时问步长大小的唯一限制条件是，允许它对于激励函数有一个接近的逼近，并以较密的时间间隔提供反应结果，以使反应峰值不会被漏掉。这种数值方法对于激励由紧密的时间间隔定义的情况（例如对于地震地面加速度的情况）特别有用，从而使得线性捕值即可得到较完美的结果。如果时间步长是常数，则系数A，B，仅需计算一次。这种数值方法所要求的运动方程精确解答仪对线性体系是可行的。如上所述，这种方法用于单白由度体系较便利，但是对于多自由度体系则是不切实际的，除非它们的反应由振型反应的叠加（第12章和第13章）来获得。5.3 中心差分法这种方法是基于对位移时间导数（

7、即速度和加速度）的有限差分近似进行的。步长，则时刻的速度和加速度的中心差分表达式为 (5.3.1)将速度和加速度的这些近似表达式代人方程(5.1.3)中，对线弹性体系，得 (5.3.2)在这个方程中，和假定是已知的（来自于前面时间步内方法的执行）。将这些已知量移到右侧，导得 (5.3.3)或写成 (5.3.4)其中 （5.3.5） （5.3.6）则未知的由下式给出 （5.3.7）+l时刻的解答是根据时刻的平衡条件即方程(5.1.3)确定的，而不是以时刻+1的平衡条件式(5.1.4)确定的，这种方法称为显式方法。观察（5.3.6），为了计算，需要已知的位移因而，为了确定，需要和。特定的初始位移是

8、已知的，为了确定，我们将式（5.3.1）专门用于=0的情况，得 （5.3.8）从第一个式子解出,然后代入第二个式子,给出 （5.3.9）其中初位移和初速度是已知的,由0时刻()的运动方程 可以得到0时刻的加速度为: （5.3.10）表5.3.1总结了可在计算机上执行的上述方法。 表5.3.1 中心差分法如果时间步长选取得不够短，那么由于数字舍入误差的存在，中心差分法将会“放大”，而给出无意义的结果。为了稳定性，特别要求 （5.3.11）对于单自由体系，上式永远不会是一个约束，因为为了获得准确的结果，选择的时间步长将是非常小的。为了充分地定义反应，通常选择；在大多地震反应分析中甚至选择更短时间步

9、长，为了准确地定义地面加速度，通常选取 =0.01到0.02秒。5.4 Newmark法5.4.1 基本方法1959年，NM. Newmark发展了一类时间步进法，它们基于下面的公式： （5.4.1a） （5.4.1b）参数和定义了时间步内加速度的变化，并决定方法的稳定性与精度特征。对于为1/2和1/6l/4的典型选择，从包括精度的所有观点来看都是令人满意的。这两个等式与时间步结束时的平衡方程(5.1.4)结合，提供了从时刻已知的、和。计算时刻的、和的基础。执行这些计算需要迭代，因为未知的出现在式(5.4.1)的右侧。然而，对于线性体系，修正Newmark法的原始公式可以允许使用式(5.4.1

10、)和(5.1.4)求解时不迭代。在阐述这些修正之前，我们先来论证Newmark法的两种特殊情况，即众所周知的平均加速度法和线性加速度法。5.4.2 特殊情况对于这两种情况，表5. 4.1总结了+1时刻的反应、和与时刻相应量之间的关系。式(5.4.2)描述了加速度在步长内变化是常数（等于平均加速度）或线性的假定。对进行积分，给出时间步长内速度的变化，即式(5.4.3)；将代入，得+1时刻的速度，即式(5.4.4)。对进行积分，给出时间步长内的位移的变化，即式(5.4.5)；将代人，得+l时刻的位移，即式(5.4.6)。将式(5.4.4)和式(5.4.6)与式(5.4.1)比较，可见Newmark

11、方程在=1/2和=1/4时与平均常加速度假定推导的那些方程是相同的；= 1/2和=l/6时的方程则与加速度线性变化的假定相符。 表5.4.1 平均加速度法和线性加速度法5.4.3 非迭代表达式我们现在返回到式(5.4.1)，为避免迭代，对其重新进行列式，并使用增量 (5.4.7) (5.4.8)尽管对于线性体系的分析增量形式是不必要的，但是引入它是因为这种形式可以方便地扩展到非线性体系。式(5.4.1)可以重新写为 (5.4.9)求解第二个等式，可得 (5.4.10)将式(5.4.10)代入式(5.4.9a)中，得 (5.4.11)下面，对具有和的线性体系，由方程(5.1.4)减去(5.1.3

12、)，得到增量运动方程： (5.4.12)将式(5.4.10)和式(5.4.11)代人方程(5.4.12)中，得 (5.4.13)式中 (5.4.14) (5.4.15)由体系的特性,算法参数和,以及时间步长开始的和可知和,则增量位移由下式计算 (5.4.16)一旦求出,则和就可以根据式(5.4.1 1)和式(5.4.10)分别计算出来，从而、和可从式(5.4.7)计算出来。加速度也能从+l时刻的运动方程确定： (5.4.17)这要好于用式(5.4.10)和式(5.4.7)确定的值。式(5.4.17)开始计算需要得到 见式(5.3.10)。在Newmark法中，+l时刻的解由式(5.4.12)确

13、定，这是与使用+1时刻的平衡条件方程(5.1.4)等价的。这种方法称为隐式法。表5.4.2总结了可在计算机上执行的用Newmark法进行时间步进求解的过程。如果时问步长满足下式，则Newmark法是稳定的: (5.4.18)对于 = 1/2，=1/4，这个条件成为 (5.4.19a)这意味着，平均加速度法对任何都是稳定的，无论它有多大；然而,如5.3结束时所讨论的，仅当足够小时才是准确的。对于=12和1/6，式(5，4.1 8)表明，如果下式成立，则线性加速度法是稳定的 (5.4.19b)然而，像中心差分法的情况一样，这个条件任单自山度体系分析中没有什么意义，因为为了获得激励和反应的正确表示，

14、必须使用比短得多的时间步长。 表5.4.2 Newmark法:线性体系5.5 稳定性与计算误差5.5.1 稳定性如果时间步长比稳定性界限短，则称导致有界解答的数值方法为条件稳定方法。不管时间步长的长短，均导致有界结果的方法称为无条件稳定方法。平均加速度法是无条件稳定的，线性加速度法在时是稳定的，中心差分法在时是稳定的。显然，后两种方法是条件稳定方法。在单自由度体系分析中，稳定性准则不是限制性的（即，它们并不规定时间步长的选择）。这是因为，为保证数值结果有适当的精度。必须比稳定性界限小得多（比方说0.1或更小），然而，数值方法的稳定性在多自由度体系的分析中则是重要的，经常需要使用无条件稳定方法（

15、第1 5章）。5.5.2 计算误差误差是运动方程任何数值解答中所固有的，我们不从数学的观点来讨论误差分析。为了对误差的本质有一个感性认识，我们先来分析数值解的两个重要特征，然后介绍一种处理误差的简单而有效的方法。考虑自由振动问题 和 理论解为 (5.5.1)采用四种数值方法：中心差分法、平均加速度法、线性加速度法和Wilson法对这个问题进行求解。其申，最后一种方法将在第15章中介绍。取时所得的数值结果与理论结果进行比较，如图5.5.1所示。这些比较说明有些数值方法可能会预测出位移幅值会随时间衰减，尽管体系是无阻尼的；并且会预测出固有周期被延长或缩短。 图5.5.1 采用四种数值方法()和理论

16、解的真由振动解答图5.5.2示出了四种数值方法中振幅衰减AD和周期延长PE作为函数的图形，AD和PE分别在图b和图c中定义，不过，没有介绍导出这些数据的数学分析。有三种方法显示位移幅值无衰减。但是，Wilson法含有幅值的衰减，意味着该方法在体系中引入了数值阻尼，等效粘滞阻尼比示于图a中。注意，中心差分法的周期误差在接近方法的稳定性界限时迅速增加。中心差分法引入了最大的周期误差，在这个意义上，它被认为是精度最差的方法。对于比它的稳定性界限小的情况，线性加速度法给出最小的周期延长。联合这个性质和幅值无衰减的特点，使得线性加速度法是单自由度体系已介绍的方法当中最合适的方法。然而，对于多自由度体系，

17、因为稳定性需求，我们将得到不同的结论（第15章）。除体系的固有振动周期以外，步长的选择还与动力激励随时间的变化有关。图5.5.2表明，取会给出相当精确的结果。为了保持激励函数的最小失真，时间步长也将是足够短。为了用数字描述地震发生时所记录的高度不规则的地震地面加速度，通常需要非常精细的时间步长，典型的时间步长为 =0. 02秒，对于计算结构反应所选择的时间步长不应该比这更长。一种虽然不复杂，但非常有用的选择时间步长的方法是：先用一个看起来合理的时间步长求解问题，然后以稍小的时间步长重复求解，对比结果，持续这个过程直到连续的两个解足够接近。前述关于稳定性和精度的讨论严格地说适用于线性体系。对于这

18、些问题是怎样影响非线性反应分析的，读者应参考其他资料。5.6 非线性反应分析：中心差分法超过线弹性范围的体泵的动力分析一般不适合用解析方法求解，即使激励随时间的变化可以用一个简单函数来描述。因此，数值方法是非线性体系分析的基本方法。可容易地对中心差分法进行修改，以求解时刻的非线性运动方程(5.1.3)。将速度和加速度的中心差分近似式(5.3.1)代人，得到式(5.3.2)，其中由代替，式(5.3.2)可重新写成求+1时刻反应的表达式 (5.6.1)其中， (5.6.2) (5.6.3)将这些等式与线性体系的那些等式进行比较，可见惟一的差别是对的定义。作此修正后，表5.3.1也适用于非线性体系。

19、抗力看起来是显式的，因为它仅取决于时刻的反应，而不依赖于未知的+1时刻的反应。因此它容易计算，使得中心差分法对非线性体系来说或许是最简单的方法。尽管在这方面吸引人，但是由于有更有效的方法可用，因此在实际应用或者研究应用中这种方法并不流行。5.7 非线性反应分析：Newmark法本节将5.4所述的线性体系的Newmark法扩展到非线性体系。尽管不像中心差分法那样简单，但是由于它有非常高的精度，因此这种方法也许足地震反应分析中最流行的方法。式(5.1.3)和式(5.1.4)之差给m增量平衡方程： (5.7.1)增量抗力为 (5.7.2)其中，图5.7.1中所示的割线刚度不能确定，因为u，是未知的。

20、如果我们做这样的假定，在一个微小的时间间隔 内，割线刚度可用切线刚度代替，如图5. 7.1所示，那么，式(5,7.2)可近似地表示为 图5.7.1 (5.7.3)去掉式(5.7.3)中的下标T，并将该式代入方程(5.7.1)中，得 (5.7.4)这个方程与线性体系的相应方程(5.4.1 2)之间的相似之处启示我们，前面对线性体系介绍的Newmark法的非迭代公式也可以用于非线性反应分析。所有需要做的是将式(5.4.14)中的用每个时问步开始时求出的切线刚度代替。这个变化意味着表5. 4.2中的步骤1.3应该移到步骤2.1的后面。对于非线性体系，步骤2.5和式(5.4.17)会给出最的不同值，后

21、者给出的值更好一些，因为它满足+1时刻的平衡方程。这个方法用不变的时间步长会导致不能接受的错误结果。重大的误差起源于两个原因：(1)用切线剐度代替割线刚度；(2)常数时间步长推迟了力一变形关系中转折点的发现。首先，我们考虑误差的第二个来源，用图5.7. 2a所示的力变形关系来说明。假设时间步开始时的时刻处位移为，速度是正的（即位移是渐增的），对应图中的a点。对时间步应用前面描述的数值方法，求得+l时刻的位移和速度，这是图中的b点。 图5.7.2如果是负的，那么在时间步内的某点，速度为零并将改变符号，位移开始减少。在数值方法中，如果我们不想麻烦找到点，而是在b点开始下一个时间步继续计算，并使用与

22、力一变形图的卸载分支相关的切线刚度，那么这个方法将在下一个时同步结束时定位于c点，位移为，速度为负的。另一方面，如果能够确定与点相关的时间瞬时（当速度实际成为零时），那么下一个时间步的计算将从体系在的状态开始，给出时间步结束时的位移和速度，记为。不定位，具有超越到b，以及不跟随力一变形图上精确路径的影响。这些与精确路线的偏离将发生在每一个速度反向处，从而导致数值结果中的误差。类似的问题也出现在力一变形关系（例如弹塑性体系）的其他尖角处。通过精确定位可以避免这些误差，这可以南折回去以更小的时间步长(比如说)在到时间间隔内的积分来达到。另外，还可以使用一个迭代过程，用比整个时间步长还小的步长，从时

23、间点重新开始积分，对积分步长进行连续调整，以使这样调节的步长在结束时的速度接近于零。现在，我们回到误差的第一种原因上来，它与使用切线刚度代替未知的割线刚度有关，如图5.7.2b中的力一变形关系所示。在时间步开始时刻i的位移为图示的口点，在血点使用切线刚度，从时刻到时刻+1的数值积分导出位移，标识为点b。如果我们能沿若正确曲线，那么结果会是点对应的位移。这个偏差经过一系列时间步的累积，会引入非常大的误差。这些误差可用迭代方法使其最小化。在Newmark法的每一步中要求解的关键方程是式(5.4.13)，对于非线性体系，经改造成为： (5.7.5)式中，由式(5-4.15)给出， （5.7.6）为符号的方便起见，我们将的下标去掉，用T代替，以强调这是切线刚度；还有和的下标也去掉。式(5.7.5)和(5.7.6)成为 (5.7.7) （5.7.8）图5.7.3a给出了式(5.7.7)的示意图。关系是非线性的，因为切线刚度取决于位移，因此斜率不是常数。在非线性体系的静力分析中，=，的非线性特性与的相同。在动力分析中，中的质量和阻尼项的存在减少了非线性特性，因为常数项对于的典型值通常比大许多。 图5.7.3非线性体系一个时间步内的迭代：(a)修正的Newton - Raphson迭代；(b)Newton - Raphson迭代接下来参考图5.7. 3a叙述迭代方法。第一个