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

5章 动力反应的数值计算如果激励[作用力)(t p 或地面加速度)(t ug ]是随时间任意变化的，或者体系是非线性的，那么对单自由度体系的运动方程进行解析求解通常是不可能的。

这类问题可以通过数值时间步进法对微分方程进行积分来处理。

在应用力学广阔的学科领域中，有关各种类型微分方程数值求解方法的文献（包括几部著作中的主要章节）浩如烟海，这些文献包括这些方法的数学进展以及它们的精度、收敛性、稳定性和计算机实现等问题。

然而，本章仅对在单自由度体系动力反应分析中特别有用的很少几种方法进行简要介绍，这些介绍仅提供这些方法的基本概念和计算算法。

尽管这些对许多实际问题和应用研究已经足够了，但是读者应该明白，有关这个主题存在大量的知识。

5.1 时间步进法对于一个非弹性体系，欲采用数值求解的运动方程为)(),(t p u u f u c um s =++ 或者 )(t u m g - (5.1.1) 初始条件)0(0u u = )0(0u u= 假定体系具有线性粘滞阻尼，不过，也可以考虑其他形式的阻尼（包括非线性阻尼），后面会明显看到这一点。

然而由于缺乏阻尼信息．因此很少这样做，特别是在大振幅运动时。

作用力)(t p 由一系列离散值给出: )(i i t p p = ,0=i到N 。

时间间隔i i i t t t -=∆+1 （5.1.2）图5.1.1 时间步进法的记号通常取为常数，尽管这不是必需的。

在离散时刻i t （表示为i 时刻）确定反应，单自由度体系的位移、速度和加速度分别为i u 、i u 和i u 。

假定这些值是已知的，它们在i 时刻满足方程i i s i i p f u c um =++)( （5.1.3） 式中，i s f )(是i 时刻的抗力，对于线弹性体系，i i s ku f =)(，但是如果体系是非弹性的，那么它会依赖于i 时刻以前的位移时程和速度。

将要介绍的数值方法将使我们能够确定i +1时刻满足方程(5.1.1)的反应1+i u 、1+i u 和1+i u ，即在i +1时刻1111)(++++=++i i s i i p f u c um （5.1.4） 对于i =0，1，2，3，…，连续使用时间步进法，即可给出i =0，l ，2，3，…所有瞬时所需的反应。

已知的初始条件)0(0u u =)0(0u u =和提供了起动该方法的必要信息。

从i 时刻到i +1时刻的步进一般不是精确的方法，许多在数值上可以实现的近似方法是可能的。

对于数值方法，有三个重要的要求：(1)收敛性一随着时间步长的减少，数值解应逼近精确解；(2)稳定性一在存在数值舍入误差的情况下，数值解应是稳定的；(3)精度一数值方法应提供与精确解足够接近的结果。

这些重要的问题在本书中均作简要的讨论，全面的论述可在着重微分方程数值解法的书中找到。

本章介绍三种类型的时间步进法：(1)基于激励函数插值的方法；(2)基于速度和加速度有限差分表达的方法；(3)基于假设加速度变化的方法。

前两类中各只介绍一种方法，第三类中介绍两种方法。

5.2 基于激励插值的方法对于线性体系，通过在每个时间间隔里对激励进行插值，并利用第4章的方法进行精确求解，能推导出一种非常有效的数值方法。

如果时间间隔较短，则线性插值是令人满意的。

图5. 2.1所示的时间间隔1+≤≤i i t t t ，激励函数为ττiii t p p p ∆∆+=)( （5.2.1a ） 其中i i i p p p -=∆+1 （5.2.1b ）时间变量τ从0到i t ∆变化。

为数学上简单起见，我们首先考虑无阻尼体系，后面再将该方法扩展到有阻尼体系。

待求解的方程为τiii t p p ku um ∆∆+=+ （5.2.2） 在时间间隔i t ∆≤≤τ0内，反应)(τu 为三部分之和：(1) τ=0时刻的初位移iu 和初速度i u引起的自由振动；(2)零初始条件下对阶跃力i p 的反应；(3)零初始条件下对斜坡力τ)(i i t p ∆∆的反应。

对这三种情况分别采用来自§2.1、§4.3和§4.4中已有的解答，得)sin ()cos 1(sin cos )(i n n i i n i n ni n i t t k p k p uu u ∆-∆∆+-++=ωτωττωτωωτωτ （5.2.3a ） )cos 1(1sin cos sin )(τωωτωτωωτωωτn in i n i n ni n i nt k p k p uu u-∆∆+++-= （5.2.3b ） 计算i t ∆=τ时的这些等式，得i +1时刻的位移1+i u 和速度1+i u ：[][])sin(1)cos(1)sin()cos(1i n i n i n i i n ii n ni i n i i t t t k p t k p t ut u u ∆-∆∆∆+∆-+∆+∆=+ωωωωωωω (5.2.4a)[])cos(11)sin()cos()sin(1i n in ii n i i n ni i n i ni t t k p t k p t ut u u ∆-∆∆+∆+∆+∆-=+ωωωωωωω(5.2.4b)将式（5.2.1b ）代入后，可将这些等式重写为如下的递推公式：11+++++=i i i i i Dp Cp uB Au u (5.2.5a) 11++'+'+'+'=i i i i i p D pC u B u A u(5.2.5b) 对于欠临界阻尼体系（即1<ζ），重复上面的推导，表明式(5.2.5)也适用于有阻尼体系，系数A ，B ，…,D '的表达式由表5.2.1给出，这些系数取决于体系的参数n ω、k 和δ以及时间间隔i t t ∆≡∆。

表5.2.1 递推公式中的系数（1<ζ）因为递推公式是从运动方程的精确解推导出的，因此对时问步长t ∆大小的唯一限制条件是，允许它对于激励函数有一个接近的逼近，并以较密的时间间隔提供反应结果，以使反应峰值不会被漏掉。

这种数值方法对于激励由紧密的时间间隔定义的情况（例如对于地震地面加速度的情况）特别有用，从而使得线性捕值即可得到较完美的结果。

如果时间步长t ∆是常数，则系数A ，B ，…，D '仅需计算一次。

这种数值方法所要求的运动方程精确解答仪对线性体系是可行的。

如上所述，这种方法用于单白由度体系较便利，但是对于多自由度体系则是不切实际的，除非它们的反应由振型反应的叠加（第12章和第13章）来获得。

5.3 中心差分法这种方法是基于对位移时间导数（即速度和加速度）的有限差分近似进行的。

步长t t i ∆=∆，则时刻的速度和加速度的中心差分表达式为t u u ui i i ∆-=-+211 tu u u u i i i i ∆+-=-+2211 (5.3.1)将速度和加速度的这些近似表达式代人方程(5.1.3)中，对线弹性体系，得 i i i i i i i p ku t u u c t u u u m=+∆-+∆+--+-+2)(211211 (5.3.2)在这个方程中，i u 和1-i u 假定是已知的（来自于前面时间步内方法的执行）。

将这些已知量移到右侧，导得112222()2()2()i i i i m c m c m u p u k u t t t t t +-⎡⎤⎡⎤⎡⎤+=----⎢⎥⎢⎥⎢⎥∆∆∆∆∆⎣⎦⎣⎦⎣⎦ (5.3.3) 或写成1ˆˆi i ku p += (5.3.4) 其中()2ˆ2mcktt =+∆∆ （5.3.5） ()()1222ˆ2i i i im c m p p u k u t t t -⎡⎤⎡⎤=----⎢⎥⎢⎥∆∆∆⎢⎥⎢⎥⎣⎦⎣⎦（5.3.6） 则未知的1+i u 由下式给出 1ˆˆi i pu k+=（5.3.7） i +l 时刻的解答是1+i u 根据i 时刻的平衡条件即方程(5.1.3)确定的，而不是以时刻i +1的平衡条件式(5.1.4)确定的，这种方法称为显式方法。

观察（5.3.6），为了计算1+i u ，需要已知的位移1-i i u u 和因而，为了确定i u ，需要0u 和1-u 。

特定的初始位移0u 是已知的，为了确定1-u ，我们将式（5.3.1）专门用于i =0的情况，得1102u u u t --=∆ ()101022u u u u t --+=∆ （5.3.8） 从第一个式子解出1u ,然后代入第二个式子,给出 ()()210002t u u t u u -∆=-∆+（5.3.9）其中初位移0u 和初速度0u 是已知的,由0时刻(00=t )的运动方程 0000p ku u c um =++ 可以得到0时刻的加速度为:mku uc p u0000--= （5.3.10）表5.3.1总结了可在计算机上执行的上述方法。

表5.3.1 中心差分法如果时间步长选取得不够短，那么由于数字舍入误差的存在，中心差分法将会“放大”，而给出无意义的结果。

为了稳定性，特别要求π1<∆n T t （5.3.11） 对于单自由体系，上式永远不会是一个约束，因为为了获得准确的结果，选择的时间步长将是非常小的。

为了充分地定义反应，通常选择1.0≤∆n T t ；在大多地震反应分析中甚至选择更短时间步长，为了准确地定义地面加速度)(t u g ，通常选取 t ∆=0.01到0.02秒。

5.4 Newmark 法5.4.1 基本方法1959年，N ．M. Newmark 发展了一类时间步进法，它们基于下面的公式：()[]()111++∆+∆-+=i i i i u t u t u uγγ （5.4.1a ） ()()()[]()[]12215.0++∆+∆-+∆+=i i i i i u t u t ut u u ββ （5.4.1b ） 参数β和γ定义了时间步内加速度的变化，并决定方法的稳定性与精度特征。

对于γ为1/2和1/6≤β≤l/4的典型选择，从包括精度的所有观点来看都是令人满意的。

这两个等式与时间步结束时的平衡方程(5.1.4)结合，提供了从i 时刻已知的i u 、i u和i u 。

计算1+i 时刻的1+i u 、1+i u 和1+i u 的基础。

执行这些计算需要迭代，因为未知的1+i u出现在式(5.4.1)的右侧。

然而，对于线性体系，修正Newmark 法的原始公式可以允许使用式(5.4.1)和(5.1.4)求解时不迭代。

在阐述这些修正之前，我们先来论证Newmark 法的两种特殊情况，即众所周知的平均加速度法和线性加速度法。

5.4.2 特殊情况对于这两种情况，表5. 4.1总结了i +1时刻的反应1+i u 、1+i u和1+i u 与i 时刻相应量之间的关系。