非线性动力学方程的求解方法1、概述在工程实际问题中，我们常常面临这样的选择：我们所遇的问题究竟是静力的还是动力的。

静力问题与动力问题，从力学的角度看就是是否考虑与加速度有关的力，而从数学求解方法看则是一个三维边值问题还是一个四维边值-初值问题。

在这个问题的选择上没有固定的原则，一般取决于我们研究者、分析者对工程问题的判断。

一般认为，实际工程大都是处于动力环境之中，因而属于动力问题。

但是，由于时间、经费等方面的原因的限制，我们不可能把所有的问题都按照动力问题的方法来分析。

对于许多具体的问题，与速度和加速度有关的力足够小，但是又影响结构分析结果的，将采用静力假定来模拟这些力。

线性的动力有限元控制方程如式(1-1)所示。

[]}{}]{[}]{[}{R q K q D qM =++ (1-1) 式中[M ][D ][K ]分别为结构的质量、阻尼和刚度矩阵，{R }为荷载列矢量，}{q、}{q 和}{q 分别是加速度、速度和位移列矢量。

式(1-1)的解法大体上可以分为两类：直接积分法和模态叠加法。

直接积分法在对控制方程进行数值积分之前不对方程做任何形式的变换，直接用数值积分的方法在时域上一步一步地对方程进行积分。

模态叠加法是在求解之前对方程进行某种数学变换，使基底降低，或使矩阵的带宽减小，再进行求解。

这两种方法在形式上不同，但是密切相关。

上述每一类求解方法中又有许多具体的解法，每一种解法又有各自的特点。

因此我们在选择一种方法求解一个问题时，要对该方法的收敛性、稳定性、效率、精度和费用等进行一些分析，讨论它对所求问题的有效性，从而使我们能够针对某一特定的问题，选择合适的方法。

直接积分法基于以下两条：(1)不是在求解时间区间内任意时刻t 都满足式(1-1)，而是在相隔△t 上的一些离散时刻满足式(1-1)。

(2)对位移、速度和加速度在每一时间区间△t 内变化的形式进行假设，事实上若把式(1-1)看成一个常系数微分方程组，便可以用任何一种有限差分格式通过位移来近似表示速度和加速度，因此不同的差分格式就得到不同的方法。

从差分格式上看，分显式与隐式两大类方法，所谓显式差分法就是不必对方程进行求解，而是由前一时刻t 的平衡条件直接可求解△t 增量后t+△t 时刻的各参数解。

而隐式差分法则必须对方程进行求解。

显式差分法有中心差分法等。

隐式差分法有Wilson -θ法和Newmark 法。

不同方法解的精度、稳定性、收敛性、效率不同。

2、中心差分法2.1中心差分法简述中心差分法的差分格式为：}]{}{2}[{1}{2t t t t t t q q q tq ∆+∆-+-∆= (2-1) }]{}{[21}{t t t t t q q tq ∆+∆-+-∆= (2-2) 认为}{t q和}{t q 在任意时刻t 上满足平衡方程(1-1)，即： []}{}]{[}]{[}{t t t t R q K q D qM =++ (2-3) 将式(2-1)和式(2-2)代入到式(2-3)中，可得：}]){[21][1(}]){[1]([}{}]{[21][1(222t t t t t t q D tM t q M t K R q D t M t ∆-∆+∆-∆-∆--=∆+∆ (2-4)2.2中心差分法的特点(1)}{t t q ∆+的解是基于利用在时刻t 的平衡条件，即}{t t q ∆+是在假定t 时刻式(2-3)成立的条件下来计算的，该积分过程称为显式积分法，因此中心差分法为显式差分法。

但是，用式(2-4)计算}{t t q ∆+时，不仅要用到}{t q 而且还要用到}{t t q ∆-。

所以，计算在时刻t+△t 的解，必须存在一个具体的起始过程，即必须利用初始条件。

这样，已知}{0q 和}{0q，可由式(2-4)求得}{0q，由式(2-1)和式(2-3)可求得}{t q ∆-，如式(2-5)所示。

}{2}{}{}{0200q t q t q q t ∆+∆-=∆- (2-5) 我们称式(2-5)为差分格式(2-4)的初始条件。

(2)在求解时不需对刚度矩阵[K ]进行三角分解。

在应用此法时一般采用集束质量矩阵或称对角质量阵，而阻尼矩阵也通常为对角形式，这样，在运用式(2-4)时，就不需要对等号右端的系数矩阵t D t M ∆+∆2/][/]([2)进行三角分解，从而可以节省计算时间。

如果不考虑系统的阻尼，则0][=D ，这时式(2-4)可简化为：}{}]){[1(2t t t R q M t =∆∆- (2-6) 其中： }]{[1}]){[2]([}{}{22t t t t t q M t q M t K R R ∆-∆-∆--= (2-7) 从式(2-6)和式(2-7)可以看出，若质量阵[M ]为对角矩阵，即在系统中只考虑集束质量的情况，则方程(2-6)实际上是解耦的，即方程组中是相互独立的，因此只需进行矩阵相乘就可得到。

把式(2-6)改写为：}{][}{12t t t R M t q -∆-∆=写成分量形式，即：)(2)(i t ii i R m t q tt ∆=∆+ (2-8) 其中，)(i t t q ∆+和)(i t R 分别为}{t t q ∆+和}{t R 的第i 个分量，而ii m 为质量矩阵的第i 个对角线元素。

对于0=ii m ，我们可以用静凝聚的方法进行处理。

若对总体刚度和质量矩阵不必进行三角分角，就意味着不必形成总体的[K ]和[M ]，因此，式(2-7)的计算只需在单元一级进行即可，即：}){2}]({[1})]{([}{}{2t t t e e e e e t t q q M tq K R R -∆--=∆-∑∑ (2-9) 式中[K e ]和[M e ]分别为与待求节点自由度有关的单元刚度矩阵与质量矩阵。

所以，使用式(2-8)和式(2-9)构成的中心差分具有明显的优点。

它不需要计算总体刚度矩阵[K ]和质量矩阵[M ]，求解过程在单元一级进行，计算效率较高。

而且，如果所有单元刚度和质量矩阵都相同时，或有大部分的单元都相同时，可以用一个或少数单元矩阵来求解整个系统。

在工程中，常会遇到重复对称式结构，这时用这种方法来求解将具有很高的效率。

(3)加速度和速度的差分格式都具有(△t 2)阶的误差。

(4)可以证明，中心差分法所要求的步长为： πncr T t t ≤∆≤∆ (2-10)其中，T n 为有限无限无集合体的最小周期，n 为单元系统的阶，cr t ∆为临界时间步长。

周期T n 可利用如下任何一种方法来计算：(1) 广义Jacobi 法；(2) 逆迭代法；(3) 正迭代法；(4) Rayleigh 商迭代法；(5) 多项式迭代法；(6) 斯图姆序列对分法；(7) 行列式搜索法；(8) 子空间迭代法。

事实上，T n 为系统最大特征值的倒数。

要求使用的时间步长△t 小于临界时间步长cr t ∆的差分格式，称为条件稳定的。

若使用一个大于cr t ∆的时间步长，则积分将是不稳定的，此时由数值积分在计算机上导致的舍入误差会增大，从而使结果失去意义。

中心差分法是一种条件稳定的方法，这一缺点导致了它在使用上具有很大的局限性。

首选是时间步长选择困难，当矩阵阶数较高时，而且同时原系统的质量分布不均匀，将使得[M ]矩阵的个别对角上的元素很小，从而使系统的高阶特征趋于无限大，T n 趋于零，这样便使积分过程成为不可能。

所以，仅仅因为一个十分小的质量元素，就会导致计算效率的降低。

当这种情况出现时，通常采用静凝聚的方法使矩阵的阶数降低，但这仅仅能在低频段模拟原系统。

因此，中心差分法是条件稳定的非线性动力学分析方法，有其局限性。

3、Wilson -θ法Wilson -θ法是线性加速度方法的扩展，基于以下两条假定：(1)当时间从t 至变化到时间(t t ∆+θ)时，加速度从t q变化为t t q ∆+θ ； (2)在时间区段t ∆θ内，结构的刚度、阻尼、地面运动加速度都无改变。

则由上述两条假定可得：}){}({}{}{t t t t t q q tq q -∆+=∆++θτθτ (3-1) 式中，θ为参数，可根据积分的稳定性和精度进行选取，一般要求大于1。

对式(3-1)两端的τ进行积分，可得：}){}({2}{}{}{2t t t t t t q q tq q q -∆++=∆++θτθττ (3-2) 再对式(3-2)两端的τ进行积分，可得：}){}({6}{21}{}{}{22t t t t t t t q q tq q q q -∆+++=∆++θτθτττ (3-3) 在式(3-2)和式(3-3)中令t ∆=θτ，则可得在(t t ∆+θ)时刻的速度和位移分别如式(3-4)和式(3-5)所示：}){}({2}{}{t t t t t t q qtq q +∆+=∆+∆+θθθ (3-4) }){2}({6}{}{}{22t t t t t t t q qt t q q q +∆+∆+=∆+∆+θθθθ (3-5) 将加速度}{t t q∆+θ 和速度}{t t q ∆+θ 写成由位移}{t t q ∆+θ表达的形式，如式(3-6)和式(3-7)所示：}{2}{6}){}({6}{22t t t t t t t q q t q q t q -∆--∆=∆+∆+θθθθ (3-6) }{2}{2}{}({3}{t t t t t t t q t q q q t q ∆---∆=∆+∆+θθθθ (3-7) Wilson -θ法认为在(t t ∆+θ)时刻，系统满足动力平衡方程(1-1)，并由Wilson -θ法的第二条假定，我们把式(3-6)和式(3-7)代入到式(1-1)中，可得到式(3-8)。

[]}{}]{[}]{[}{t t t t t t t t R q K q D qM ∆+∆+∆+∆+=++θθθθ (3-8) 其中，}){}({}{}{t t t t t t R R R R -+=∆+∆+θθ (3-9)求(t t ∆+θ)时刻的位移、速度和加速度是我们所要得到的。

求解的过程如下所示：(1) 将式(3-6)和式(3-7)代入到式(3-8)中，并加以整理可得：}]){[3][)(6]([2t t q D tM t K ∆+∆+∆+θθθ =}){2}{6}{)(6]([}){}({}{2t t t t t t t q q t q t M R R R +∆+∆+-+∆+θθθ +}){2}{2}{3]([t t t q t q q t D ∆++∆θθ (3-10) 由式(3-10)可以求得}{t t q ∆+θ。

(2) 将式(3-6)分别代入式(3-1)、式(3-2)和式(3-3)中，并令t ∆=τ，可得：}){31(}{)(6}){}({)(6}{22tt t t t t t q q t q q t q θθθθ-+∆--∆=∆+∆+ (3-11) }){}({2}{}{t t t t t t q q t q q +∆+=∆+∆+ (3-12) }){2}({6}{}{}{2t t t t t t t q q t q t q q +∆+∆+=∆+∆+ (3-13) Wilson -θ法的特点：(1) 该法是一种隐式积分法，因为刚度矩阵[K ]是未知位移矢量的系数矩阵。