第三章 离散化结构动力方程解法

第三章 离散化结构动力方程的解法

§3.1 绪 言

对于一个实际结构，由有限元法离散化处理后，应用瞬时最小势能原理可导出动力方程

?M??u???C??u???K??u???F(t)? （3.1） 这里，?u?、?u?、?u?及?F(t)?分别表示加速度、速度、位移及所作用的外力矢量，他们都是与时间有关的。

从数学的角度来看，式（3.1）是一个常系数的二阶线性常微分方程组，对于它的求解原则上并无困难。但是，由于?M?、?C?和?K?的阶数非常高，使得式（3.1）的求解必须花费很大的代价，便促使人们去寻求一些效率高的近似计算方法。目前，用于求解式（3.1）的方法，大致可分为两大类。

一是坐标变换法，它是对结构动力方程式（3.1），在求解之前，进行模态坐标变换，实际上就是一种Ritz变换，即把原物理空间的动力方程变换到模态空间中去求解。现在，普遍使用的方法是模态（振型）迭加法，即用结构的前q阶实际主模态集（主振型阵）构成坐标变换阵进行变换。通过这一变换，实现降阶，求较好的近似解，而且，还用解除耦合的办法，简化方程的计算。还有一种所谓假设模态法，即是用一组假设模态，构成模态坐标变换阵进行变换，获得一组降阶的而不解耦的模态基坐标方程。显然，这种方法的计算精度，取决于所假设的模态。用Ritz矢量法求解的近似模态作为假设模态，可得到满足要求的精度。

二是直接积分法，它是对式（3.1）在求解之前，不进行坐标变换，直接进行数值积分计算。这种方法的特点是对时域进行

1

离散，将式（3.1）分为各离散时刻的方程，然后，将该时刻的加速度和速度用相邻时刻的各位移线性组合而成，于是，式（3.1）就化为一个由位移组成的该离散时刻上的响应值，通常又称为逐步积分法。线性代数方程组的解法与静力时刻的位移来线性组合，就导致了各种不同的方法。主要有中央差分法，Houbolt方法，Wilson-θ法和Newmark方法等。

§3.2 模态（振型）迭加法

设有n个自由度的系统，在外力?F(t)?的作用下，常常被激起较低阶的一部分模态（即振型），而绝大部分高阶模态被激起的分量很小，一般可忽略不计。例如，在地震载荷作用下，通常，只有最低的二阶，三阶模态起主要作用。所以，对于这样的一些问题，采用模态迭加法是有效的。

设有式（3.1）的n阶动力方程，起主要作用的是其前q阶模态，通常取qn。按

Ritz变换，则可将式（3.1）中的?u?用前

qq个模态的线性组合来表示，即

{u}?Y1{?1}?Y2{?2}?...?Yq{?q}??{?j}Yj

j?1 ?[?]{Y} (3.2)

其中，[?]n?q为结构的已知的保留主模态矩阵，而{Y}q×1是维的模 态基坐标矢量，它形成了一个q维的模态空间。它表示在{Y}中，各阶主模态所占有的成分的多少。

假定[?]已用第二章所述的某一方法解出，再将式（3.2）代入（3.1），并左乘以[?]T，可得

2

[M]*{Y}?[C]*{Y}?[K]*{Y}?{F}* (3.3)

式中

[M]*?[?]T[M][?][K]*?[?]T[K][?][C]?[?][C][?]{F}*?[?]T{F}显然，式（3.3）是一个q阶的微分方程组。由于q?n，所以，它比式（3.1）的n阶就小的多了，实现了降阶，因而也就容易求解多了。

若展开上述的[M]*的表达式，根据主模态（主振型）关于[M]的表达式，根据主模态的（主振型）关于[M]的正交性质，可知

*T

mij*?0(i?j)

所以，[M]*是一个对角阵。同理可知[K]*也是一个对角阵。然而，在一般的情况下，[C]*是一个非对角阵，即在模态空间中，系统的的阻尼一般是耦合的。因此，式（3.3）是一个完全解耦的动力学方程。但是，它是一个已降阶的q阶的动力方程，可使用后面即将介绍的直接积分法求解。

当系统的阻尼为比例阻尼时，即[C]可以表示为

则[C]*为对角阵。此外，若系统的阻尼是一般的的线性阻尼，并非比例阻尼，但是只要结构的固有频率不相等，而且不十分接近，则可用舍去[C]*阵中的非对角元来实现[C]*的对角阵，也不会引起太大的误差。

在上述两种情况下，可以获得对于模态坐标的完全解耦的动力学方程。即式（3.3）是q个独立的方程，每个方程只包含一个未知量，相互之间不耦合。因而式（3.3）可按单自由度的动

3

[C]*??[M]*??[K]* (3.4)

力学方程写为

mii*yi?cii*yi?kii*yi?Fi*(t)(i?1,2,...q) (3.5)

或

yi+2?i?iyi+?i2yi?fi*(t)(i?1,2,...q) (3.6)

其中

2?i?i?cii*/mii*,fi(t)?Fi*(t)/mii*。式（3.6）可用直接积

分法计算，或用Duhamel积分求得其解为

yi(t)?1?ie?t0fi(?)e??i?i(t??)sin?id??--(3.7)

-??i?it{aisin?it?bicos?it}(i?1,2,...q)式中，?i??i1??i2，而ai，bi由初始条件

{y0}?([M]*)?1[?]T[M]{u0}{y0}?([M]*)?1[?]T[M]{u0}得出的

(3.8)

yi0与yi0决定。

由于有阻尼的存在，由初始条件所激发的振动，随时间的增长而衰减以致消失。因此，常可不计式（3.7）中的第二项，即是由初始条件激发的自由衰减振动。

计算出yi(t)后，便可利用式（3.2），计算出物理坐标的响应

{u(t)}。

整个计算步骤可归纳如下:

第一步：根据结构的离散化模型，建立系统的[M],[K]以及

{F(t)}，并进行结构的固有特性分析，即求解特征值问题

([K]-?2[M]){?}?{0}

4

求出前q阶特征对(?i,{?i})，（,..i.?,2,1标下的动力方程

..q）

第二步：形成模态阵[?]n?q?[{?1} {?2}...{?q}]，并建立模态基坐

yi?2?i?iyi??i2yi?fi(t) (i?1,2,...,q)

其中fi(t)?数据确定各阶主振动中的比例阻尼?i。

11T{?i}{F(t)}，而mii?{?i}T[M]{?i}。根据实验结果或经验mii.第三步：求解主模态基坐标的动力方程，有

yi(t)??i?t0fi(?)e?i?i(t-?)sin{?i(t-?)}d?，其中，

?i??i1-?2i。

第四步：进行坐标变换后，求得动力响应

?u??[?]?Y?

§3.3模态假设法

上节所述的模态迭加法，是用系统的真实主模态组成的模态矩阵，再对系统的物理坐标进行模态坐标变换，从而在主模态空间中得到降阶并解耦的动力学方程，这样来实现简化计算。而这里提出的假设模态法，则是用一组假设模态矩阵，对系统的物理坐标进行模态坐标转换，从而在模态空间中得到一组只降阶的动力学方程。

若令假设模态矩阵为?Φ?n?m，而m<

?X?t??n?1????n?m?q?t??m?1

T (3.10)

把它代入式（3.1），并左乘???，则可得到降阶的动力学方程为

5

?M??q???C??q???K??q???Q?t?? (3.11)

***其中

?M??????M????， ?K??????K????, ?C??????C????, ?Q?t???????F(t)?。

*T*T*TT它们分别对应于假设模态坐标?q?的质量矩阵、阻尼矩阵、刚度矩阵与广义力列阵。因为矩阵???中的各列都是假设模态，它们一般不具有正交性，所以，?M?、?C?和?K?都不是对角阵。于是，方程(3.11)是不能解耦的方程组，但它却是比式(3.1)的阶数要低得多了。显然，对式(3.11)采用直接积分法求解，将比对式(3.1)求解要简便得多。这是假设模态法的优点。

假设模态法的计算精度，很显然地是取决于假设模态阵中模态假设的好坏与质量。因此，应用假设模态法能否成功的关键在于确定出一个适宜的假设模态矩阵。在第五章中，我们介绍了几种构造假设模态的方法。

实际上，Rayleigh-Ritz法可认为是一种假设模态法。它的作用在于降低方程的阶数，简化计算。它的基本思想是，事先假定出若干近似的特征矢量，然后按照这些特征矢量的最佳线性组合，而算得前若干阶特征值的近似值。显然，运用这种方法时，其计算精度与事先假定的特征矢量的近似程度和数量有关。

按照Ritz变换的思想，找到了近似的特征矢量?Xi?后，

***?i=0,???,q?，即有

????a1?X1??a2?X2??????aq?Xq???X??A? (3.12)

求解如下的广义特征值问题，即

?K??A?=ρ?M??A?,???2 (3.13)

其中

6

**

?K?=?X??K??X? ?M?=?X??M??X?

*T*T?K?和?M?为原结构离散化之刚度阵和质量阵，它们都是n阶方

阵。求解式(3.13)，得到q个特征矢量，有

1a?A1?=?1?2?A2?=??a1a2a21???aq??21T

2???aq?????aqq??T??????

??qAq=?a1?a2qT再按照Ritz的变换，即式(3.12)，由特征矢量?Aj?，可计算出矢量??1?，??2?，???，??q?，即是

??i???aji?xi? ?i=1,2,???,q? (3.14)

j?1q现在用???n?q?????1???2?构造的假设模态矩阵。

???它就是我们要?????来表示此变换阵，

q§3.4 中心差分法（显示法）

现在开始讨论直接积分法，或称逐步积分法。

前面讨论的模态迭加法，并非总是有效的。当刚度矩阵?K?，或质量矩阵?M?，或阻尼矩阵?C?出现随时间变化时，或当外荷载激起的振型太多，需要计算的特征对太大时，就不宜于采用模态迭加法，在这些情况下，采用逐步积分法是适宜的。中心差分法就是其中的一种。这种方法的特点，是将动力方程在时间域上离散，化成对时间的差分格式，然后根据初始条件，利用直接积分

7

法逐步求解出一系列时刻上的响应值。

?0和??0。假定t?0时，位移、速度和加速度分别为已知的u0，u u再将求解的时间区间划分为n个等分，即?t?。我们要建立的积 分格式就是从已知的0，?t，2?t，…，t的解来计算下一个时间步的解。

在中心差分法中，按中心差分将速度和加速度矢量离散化为

Tn{u}t?1{u}t??t?{u}t??t? ?2?t（3.15） （3.16）

{u}t?1{u}t??t?2{u}t?{u}t??t? 2??t上面二式将t时刻的速度和加速度用相邻时刻的位移表示了。考虑在t时刻的动力方程，有

[M]{u}t?[C]{u}t?[K]{u}t?{F}t

（3.17）

将式（3.15）和（3.16）代入式（3.17）中，得到

121?1????1?[M]?[C]{u}?{F}?[K]?[M]{u}?[M]?[C]{u}t??t??t2?t??tt??t??222?t?t2?t??????t?（3.18）

这样，上式就化为用相邻时刻的位移表示的代数方程组。由它可解出{u}t??t。又由于它是利用t时刻的方程解得{u}t??t的，所以，它称为显示积分。并且，还注意到，在求解{u}t??t时，需要用{u}t，

{u}t??t的值。于是，在计算开始时，即t?0时，要计算{u}?t的值、

就需要{u}??t的值，他是未知的，因此，必须有一个启动的处理，

?}0和{u??}0是已知的，因而这种算法不是自起步的。由于{u}0,{u所以，

由t?0时的式（3.15）和（3.16），可解得

{u}??t?{u}0??t{u}0??t22{u}0

（3.19）

使用中心差分法的逐步求解过程如下：

8

A． 初始计算

（1）形成刚度矩阵[K]，质量矩阵[M]和阻尼矩阵[C]。

?}0和{u??}0。 （2）给定初始值{u}0,{u（3）选择时间步长?t，?t??tcr，并计算积分常数：

a0?1?t2，a1?12?t，a2?2a0，a3?1a2。

?}0?a3{u??}0。 （4）计算{u}??t?{u}0??t{u~（5）形成有效质量矩阵[M]?a0[M]?a1[C]。 ~（6）三角分解[M]：[M]?[L][D][L]T。

B． 对每个时间步计算

（1）计算t时刻的有效载荷

~{F}t?{F}t?([K]?a2[M]){u}t?(a0[M]?a1[C]){u}t??t。

（2）求解t??t时刻的位移

~[L][D][L]T{u}t??t?{F}t。

（3）如果需要计算t时刻的速度和加速度

{u}t?a0?{u}t??t?2{u}t?{u}t??t?

{u}t?a1?{u}t??t?{u}t??t?

应当指出，这种中央差分算法，左端的系数矩阵只与质量阵

[M]和阻尼阵[C]有关，而与刚度阵[K]无关。如果质量阵和阻尼阵

是对角阵，那么在解方程时，就不需要对系数阵进行三角分解，即不需要解线性代数方程组，从第一步开始逐次直接求得各个时刻{u}t??t的值，这是中央差分格式就是一种显示的格式。此外，由于不求解代数方程组，也就不需要进行组集，它的右端项的形成也只须在单元一级水平上，由每个单元对有效载荷矢量的贡献迭加而成。因此，ADINA程序规定，在用中心差分法时，必须使

9

用对角的质量阵和阻尼阵。从计算稳定性角度来看，中心差分法的缺点，在于它是条件稳定的，即当时间步长?t太大 时，积分是不稳定的。所以，对步长的限制是

?t??tcr?Tn?

这里，?tcr是临界步长值，Tn是有限元系统的最小周期。这样，当Tn很小时，就限制了?t必须很小，所以求解所花的代价就很大。

§3.5 线性加速度法和Wilson-θ法

线性加速度法和Wilson -θ法，都是属于逐步积分法。线性加速度法是假定在[t,t??t]时间间隔内，即在步长?t时间内，加速度

{u(t??)}呈线性变化，其表达式为

{u}t???{u}t??{A}

（3.20）

其中，{A}?({u}t???{u})/?t。但是，这个方法不是无条件稳定的，所以在应用上受到限制。70年代初期，Wilson推广了线性加速度法，他假定在此步长?t更大的时间区间(t,t???t)内，加速度仍保持线性变化，经过证明，当??1.37时，这一方法是无条件稳定的，这就是Wilson -θ方法。

这个方法的加速度表达式为 式中

{A1}?({u}t???t?{u}t)/??t

{u}t???{u}t???{A1}

(3.21)

显然，对比式（3.20）和式（3.21）得知，线性加速度法是Wilson -θ法中，当??1时的一个特例。所以，我们只讨论Wilson -θ法就够了。

在0?????t区间内，对式（3.21）进行积分，得到

10

和

{u}t??{u}t???2?{u}t?{u}t??({u}t???t?{u}t)

2??t（3.22）

1?32?{u}t?{u}t??{u}t???({u}t???t?{u}t)（3.23）

26??t令????t，由上二式，有 和

{u}t???t?{u}t???t{u}t?{u}t???t?{u}t???t2({u}t???t?{u}t) （3.24）

?2?t26({u}t???t?2{u}t) （3.25）

从这二式，可将(t???t)时刻的加速度和速度用位移来表示即 和

{u}t???t66?22({u}t???t?{u}t)?{u}t?2{u}t （3.26） ??t??t{u}t???t?3??t({u}t???t?{u}t)?2{u}t?{u}t ??t2(3.27)

于是，在t??t时刻的动力方程为

[M]{u}t???t?[C]{u}t???t?[K]{u}t???t?{F}t???t

(3.28)

式中，

{F}t???t?{F}t??({F}t??t?{F}t)

将(3.26)和式(3.27)代入式(3.28)，就得到关于{u}t??t的方程为

([K]?63[M]?[C]){u}t???t?{F}t??({F]t??t?[F]t)22??t??t66(3.29) ?[M](22{u}t?{u}t?2{u}t)??t??t3??t?[C]({u}t?2{u}t?{u}t)??t2 11

记[K]?[K]?63[M]?[C] ?2?t2??t{F}t???t?{F}t??({F}t???t?{F}t)?[M](6??t?2{u}t)?[C]({u}t?2{u}t?{u}t)??t266{u}?{u}tt22??t??t

于是，式(3.29)可写为

[K]{u}t???t?{F}t???t

(3.30)

求解方程(3.30)，则得到{u}t???t

将求解得到的{u}t???t，代入(3.26)中，就得到{u}t???t。如在(3.21)

??t，并将式(3.26)代入，有

中，取?

{u}t??t?663({u}?{u})?{u}?(1?){u}t t???ttt?3?t2?2?t??t(3.31)

将(3.21)代入式(3.22)和(3.23)，并取?

，有

(3.32) (3.33)

{u}t??t?{u}t?{u}t??t?t({u}t??t?{u}t) 2

?t2?{u}t??t{u}t?({u}t??t?2{u}t)

6用Wilson??法逐步求解的过程如下： A． 初始计算

（1） 形成刚度矩阵[K]，质量矩阵[M]和阻尼矩阵[C]。 （2） 给出初始值{u}0，{u}0和{u}0。

（3） 选择时间步长?t，取??1.4，并计算积分常数，

a0?6(??t)23??ta0,

a1?， a2?2a1，

a5??a2a3???t2，

a4??，

?，

12

3a6?1?

?,

?ta7?2 ,

?t2a8?6

**??KK（4）形成有效刚度矩阵?：???????K??a0?M??a1?C?

T**???KK?LDL（5）对?作三角分解： ??????????B．对每个时间步计算 （1）计算t??t时刻的有效载荷

?F?t???t??F?t????F?t??t??F?t??M?a0?u?t?a2?u?t?2?u?t?????C?a1?u?t?2?u?t?a3?u?t??

（2）计算t???t时刻的位移

?L??D??L??u?t???t??F?t???t

T（3）计算t??t时刻的位移，速度和加速度

?u?t??t?a4??u?t???t??u?t??a5?u?t?a6?u?t ?u?t??t??u?t?a7??u?t??t??u?t?

?u?t??t??u?t??t?u?t?a8??u?t??t?2?u?t?与中心差分法相比较，Wilson-?法是隐式积分，即每计算一步，必须解一个线性代数方程组。当??1.37时，它是无条件稳定的。此外，这种算法是自起步的，t??t时刻的位移，速度和加速度都可由t时刻的变量表示，不需要特别的起动处理。

§3.6 Newmark方法

Newmark在1959年提出的逐步积分格式，故称为Newmark方法。它的基本假定是

?u?t??t??u?t????1????u?t???u?t??t???t (3.34)

13

??1?2?u?u?u?t????u??u??t??t??t??t??t??t??t ???t??2???? (3.35)

其中?和?是按积分的精度和稳定性要求可以调整的参数。当??,

12??1时，它就是线性加速度法，所以，Newmark6方法也

可以理解为线性加速度法的一个小延伸。Newmark法最初提出作为无条件稳定的一种积分格式是常平均加速度法，即假定从t到t??t时刻，加速度不变，取为常数

??1。 411u?u??。此时，取?????tt??t?22,

常平均加速度法是应用得最广泛的逐步积分方法之一。

研究表明，当??0.5，??0.25?0.5???2时，Newmark方法是无条件稳定的。从式?3.34?和?3.35?可得到?u?t??t，?u?t??t 用?u?t??t及?u?t、?u?t和?u?t表示的表达式，即有

?u?t??t?和

?u?t??t?????u??t-??ut??t1??t2???ut??u???t??t11???1??u????t??t?2??? u t （3.36）

??1-?????????t?u?t(1-??)t? u

2??（3.37） t考虑t??t时刻的动力方程，有

?M??u?t??t??C??u?t??t??K??u?t??t??F?t??t （3.38）

将式（3.36）和（3.37）代入（3.38），就得到关于?u?t??t的方程为

??K???u?t??t??F?t??t （3.39）

其中

1???K?K?M??????C? ????t2??t 14

??Ft??t??F?t??t??M?(1?u?t??t ?1u2??t??t?1???-1??u?t)?2????? ??C?[u???t?-1??u?t???t?????? ??-1??t?u?t]?2??

求解方程（3.39），就可得到?u?t??t，然后，根据式（3.36）和式（3.37）可解出?u?t??t和?u?t??t。 Newmark方法逐步求解的过程如下： A. 初步计算

（1）形成刚度矩阵?K?，质量矩阵?M?和阻尼矩阵?C?。 （2）给定初始值?u?0，?u?0和?u?0

（3）选择时间步长?t，参数?和?，并计算积分常数：

??0.50,??0.25?0.5???2

1?1; ??; ??;122??t??t??t1??t????3?-1; ?4?-1; ?5??-2?;

2??2????0??6??t?1-??; ?7???t; （4）形成有效刚度矩阵??K??: ??K????K???0?M???1?C?

T???KK?LDL （5）对?作三角分解： ?????????? B．对每个时间步计算

（1）计算t??t时刻的有效载荷

?F?t??t??F?t??t??M??0?u???2?u?t?|?3?u?????C??1?u?t??4?u?t??5?u?t??

（2）求解t??t时刻的加速度和速度

15

?u?t??t??0??u?t??t-?u?t?-?2?u?t-?3?u?t

?u?t??t??u?t??6?u?t??7?u?t??t我们注意到Wilson-?法与Newmark法的计算关系式，在形式上是相同的，只是其中的系数取不同的值而已。因此，它们可用同一计算机程序来实现。

§3.7 Houbolt方法

这个差分格式是利用t+Δt、t、t-Δt、t-2Δt四个时刻上位移的三次插值多项式建立起来的。即假定

{u}t??t?1(2{u}t??t-5{u}t?4{u}t-?t-{u}t-2?t) 2?t （3.40）

和

{u}t??t?1(11{u}t??t-18{u}t?9{u}t-?t-2{u}t-2?t) 6?t （3.41）

这里认为{u}t-2?t，{u}t-?t和{u}t是已知的，而{u}t??t是未知的。考虑t+Δt时刻的动力方程，有

[M]{u}t??t?[C]{u}t??t?[K]{u}t??t?{F}t??t

（3.42）

将式（3.40）和（3.41）代入式（3.42）中，就得到求解{u}t+Δt时刻的方程为

11?2?[M]?[C]?[K])??t2?{u}t??t6?t??33?5??4??{F}t??t??2[M]?[C]?{u}t-?2[M]?[C]?{u}t-?t （3.43）

?t2?t??t???t?1?1???2[M]?[C]?{u}t-2?t3?t??t? 16

由上式解得{u}t??t后，代入式（3.40）和式（3.41）中，便求得了{u}t??t和{u}t??t。这样，逐步求下去，便可求得任意时刻的动力响应值。

应该注意到，这个差分格式不是自起步的，除了利用初始条件{u}0，{u}0和{u}0之外，尚需利用前述的任一种自起步的方法，求得{u}?t、{u}?t和{u}?t后，再利用式（3.40）和式（3.41），即

{u}?t?1(2{u}?t-5{u}0?4{u}-?t-{u}-2?t) 2?t1{u}?t?(11{u}?t-18{u}0?9{u}-?t-2{u}-2?t)

6?t由上二式，可求 {u}-?t、{u}-2?t。这样，利用{u}-?t，{u}0和{u}?t，就可由Δt开始，用式（3.43）就{u}2?t，再代入式（3.40）和式（3.41），即可求得{u}2?t和{u}2?t，如此逐步求下去，即可求得任意时刻的动力响应。

Houbolt方法逐步求解的过程如下：

A. 初始计算

形成刚度矩阵[K]，质量矩阵[M]和阻尼矩阵[C]。 给定初始值{u}0、{u}0和{u}0。 选择时间步长Δt和积分常数：

a0=a321153;a=;a=;a=;a=-; a=-2a;12354022Δt6ΔtΔtΔt2a6?a0a; a7?3 29（4）使用一种自起步的方法计算出?u??t，

????。

u.和

u..?t?t~~???（5）计算出有效刚度矩阵?K?:?K???K??a0?M??a1?C? ?????~~T???（6）对?K?作三角分解：?K? ?LDL??????????? 17

B：对每一个时间步计算 （1）计算t??t时刻的有效载荷

?~??F???F?t??t??M?a2?u?t?a4?u?t-?t?a6?u?t-2?t??C?a3?u?t?a5?u?t-?t?a7?u?t-2?t??t??t????（2）求解计算t??t时刻的位移

?L??D??L??u?t??tT?~???F? ??t??t（3）如有需要，计算t??t时刻的加速度和速度

?..??a0?u?t??t-a2?u?t-a4?u?t-?t-a6?u?t-2?t ?u???t??t?.??a1?u?t??t-a3?u?t-a5?u?t-?t-a7?u?t-2?t ?u???t??t此方法对线性动力问题是无条件稳定的，所以，在选取?t时，仅需要考虑精度要求。 此外，如果?M???0?,?C???0?，这种求解方法可提供与时间有关载荷的静力解法。

18

本文来源：https://www.bwwdw.com/article/vzwx.html