结构动力学中的常用数值方法.doc

《结构动力学中的常用数值方法.doc》由会员分享，可在线阅读，更多相关《结构动力学中的常用数值方法.doc（16页珍藏版）》请在三一办公上搜索。

1、第五章 结构动力学中的常用数值方法51结构动力响应的数值算法 当c为比例阻尼、线性问题模态叠加最常用。但当C无法解耦，有非线性存在，有冲击作用（激起高阶模态，此时模态叠加法中的高阶模态不可以忽略）。此时就要借助数值积分方法，在结构动力学问题中，有一类方法称为直接积分方法最为常用。所识直接是为模态叠加法相对照来说，模态叠加法在求解之前，需要对原方程进行解耦处理，而本节的方法不用作解耦的处理，直接求解。（由以力学，工程中的力学问题为主要研究对象的学者发展出来的）中心差分法的解题步骤1. 初始值计算（1） 形成刚度矩阵K，质量矩阵M和阻尼矩阵C。（2） 定初始值，。（3） 选择时间步长，使它满足，并

2、计算 ，（4） 计算（5） 形成等效质量阵（6） 对阵进行三角分解2对每一时间步长（1） 计算时刻t的等效载荷 （2） 求解时刻的位移 （3） 如需要计算时刻t的速度和加速度值，则若系统的质量矩阵和阻尼矩阵为对角阵时，则计算可进一步简化。纽马克法的解题步骤1.初始值计算（1）形成系统刚度矩阵K，质量矩阵M和阻尼矩阵C（2）定初始值，。（3）选择时间步长，参数、。并计算积分常数 ， ， ，（4）形成等效刚度矩阵 （5）矩阵进行三角分解 2. 对第一时间步长（1）计算时刻的等效载荷 （2）求解时刻的位移 （3）计算时刻的加速度和速度 威尔逊-法的解题步骤1. 初始值计算 （1）形成系统刚度矩阵K，

3、质量矩阵M和阻尼矩阵C（2）定初始值，。（3）选择时间步长，并计算积分常数 ， ， ，（4）形成等效刚度 （5）将等效刚度进行三角分解 2.对每一个时间步长 （1）计算时刻的等效载荷 （2）求解时刻的位移 （3）计算在时刻的加速度、速度和位移 52 结构动力响应数值算法性能分析对公式（5.1）描述的线性系统结构动力学问题，已经有证明对整个多自由度的积分，等价于将模态分解后对单自由度的积分的结果进行模态叠加，因此可以通过对单自由度问题的分析，来说明算法的特性，其中阻尼均假设为比例阻尼，这样，模态分解后的单自由度结构动力学方程为： （5-29）以下算法的性能分析，均将算法用于这个方程。521 算法

4、用于结构动力学方程的有限差分表示将数值计算方法应用于（5-29）, 即分别在相邻的不同时刻应用算法可得如下一般形式 （5-30）A为放大矩阵或称逼近算子，为载荷逼近算子。 例如将Newmak方法应用于方程（5-29）有： （5-33） ， 矩阵A的特征多项式为 （5-34）其中A1，A为该矩阵的两个特征向量，分别为矩阵的迹的一半和矩阵的行列式 （5-35） （5-36） （5-37）对Newmak方法有：, （5-38）其中 h为时间步长，。Newmak方法放大矩阵的规模是二维的，因此特征值也只有两个，可以根据它们进行分析。有的算法放大矩阵是三维的，例如Wilson-方法,在无阻尼情况下放大矩

5、阵为： （5-39）放大矩阵A的特征多项式为： （540）其中A1，A2，A3为该矩阵的三个特征向量，分别为矩阵的迹的一半、各阶主子式的和以及矩阵的行列式，对Wilson-方法有 （5-41）此外，在几个不同时刻应用数值算法，然后将方程中的速度和加速度项消去，可得数值算法关于位移的差分方程，例如Newmak方法，有 （5-42）很显然，其特征方程与其放大矩阵A的特征方程是相同的，使用关于位移的线性多步方式和放大矩阵来说明算法性能是一样的，只不过各有方便之处。522 算法的稳定性分析设为放大矩阵A的特征值，则定义为A的谱半径，若特征值互异，则的算法是稳定的，但若有重特征根，则要求。如果算法的稳定

6、性要求对步长的选取有限制，称算法是有条件稳定的，反之为无条件稳定的。放大矩阵的谱半径小于等于1成立的充分条件是 （5-43）对的放大矩阵 （5-44）上两式是关于算法自由参数的不等式，由它可以判断算法是否无条件稳定，若不是，将给出稳定条件。 例51分析Newmak方法、Wilson-方法的稳定性解： 将（5-38）代入（5-43）有显然，当 （5-45）算法无条件稳定。当且 （5-46）算法稳定，但为条件稳定，其中为临界采样频率。由于（5-43）式仅仅是充分条件，所以可进一步按照稳定性的定义得到5.1.2节叙述的无条件稳定条件。对Wilson-方法，将（5-41）代入（5-44）得 （5-47

7、）容易看出，其中第一，二，五不等式恒成立，对第三，四不等式若希望对任意的均成立，则有：求解上述不等式得 （5-48）实际使用中通常选取=1.45.2.3 算法的相容性和收敛性直接积分算法的相容性、收敛性分析同样要使用其位移型的差分方程，或对应的单步多值形式。在算法（5-30）式中，用精确解代替近似解，即可得到局部截断误差表达式，用符号表示 （5-49）局部截断误差表达式用放大矩阵的特征量以最常用的线性三步法为例可表示为 （5-50）其中分别为对应的的放大矩阵的三个特征向量，然后将 在点进行泰勒展开，然后利用运动平衡方程化简即可。若局部截断误差表达式为步长的s（s0）阶小量，则称算法是s阶相容的

8、。对的放大矩阵，可仿照上述步骤，来验证算法的相容性。在经典的数值算法收敛性分析理论中，一个重要的结论就是相容加稳定等于收敛，其相容的阶数就是算法的精度阶。收敛性的含义也是当时间步长趋于零，算法的数值解趋于精确解。对直接积分算法该定理同样可以证明是成立的。例52分析Newmak法的相容性和精度解：其局部误差仿照(550)式得： (5-51)也可以由(5-42)是直接求得，即将在点泰勒展开，并注意到在时刻的运动方程有： （5-52）显然，当物理阻尼为零时，选择算法是二阶的，即截断误差是步长的二阶小量。物理阻尼的存在，使算法精度降了一阶，但若同时选择，算法精度仍然是二阶的，一般称为Newmak线加速

9、度法。显然Newmak方法中有两个参数待定，每种特定的选取都是一个特定的算法，最常用的几个算法见表5-1表5-1： 常用的Newmak族直接积分算法方法名称稳定条件无阻尼问题精度阶有阻尼问题精度阶类型1/21/4平均加速度方法（梯形法）无条件21隐式1/21/6线性加速度方法22隐式1/20中心差分方法21显式如果在一个时间步内需要求解一个隐式的方程组，则称算法是隐式的，反之不需要求解方程，直接计算即可得到下一时刻的值，则称算法是显示的。从5.1节的Newmak方法的计算步骤可以看出，这类方法是隐式的，但对于中心差分方法，若质量矩阵和阻尼矩阵都是对角矩阵就可以显示地计算。显然显示方法计算量要小

10、得多。读者可自行分析Wilson-方法的精度，不难分析，无论是无阻尼还是有阻尼其精度都是2阶的，它也是隐式方法。5.2.4 算法耗散和弥散特性算法的精度，在小步长的情况下可以通过局部截断误差分析来说明比较，但是，在实际计算过程中，步长的选取可能不是很小，此时如何来度量算法的计算精度，当然可以针对有解析解的问题进行大量的数值计算，将数值解与解析解进行比较来分析算法的计算精度。理论上还可以通过数值耗散(disspation)和弥散（dispersion）来辅助度量与分析，为引出这两个概念的含义，我们仍然以单自由度有阻尼自由振动问题为例，该问题的解析解为： (5-53)当直接积分算法用于这样的问题，

11、前小节已经讲述过它可以写成形如(5-42)的关于位移的有限差分形式，就是可以得到一个关于位移的有限差分方程，对于一个收敛的且有一定精度的算法，这个差分方程通常有一对共扼复根， (5-54)其中，该两根称为主根，其它根称为寄生根(spurious roots)。解的一般形式可写为 (5-55)式中的称为算法阻尼比，当有物理阻尼存在时，它还包括了物理阻尼的影响，称为算法频率，对应的称为算法周期。可以看到上式前两项与(5-53)式形式是相同的，这给了我们与精确解进行比较的可能，如果 1） 寄生根的影响较小，即 。2） 解表达式中得常数差别不太大。这样，我们就可以通过比较不同算法的算法阻尼比和相对周期

12、误差，来比较算法的计算精度。此时，对无阻尼问题，可以很明显看到算法阻尼比会使得数值解曲线的幅值与解析解相比要降低而产生振幅衰减，这就是所谓的算法的数值耗散。同时，不同的算法的算法周期与精确的周期会有一定的误差，这个误差一般用相对周期误差来表示，它会使得数值解曲线上产生周期的延长或缩短，即所谓的数值弥散。实际分析时，可以首先通过求解放大矩阵的特征方程得到特征方程的主根和寄生根，若主根可以表示为： (5-56)并注意到式(5-54)有 (5-57) (5-58) (5-59) (5-60)对结构动力学问题，一般总希望算法在低频段有较小的耗散和弥散，而且在低频段，当时，可以很方便地获得它们的近似解析

13、表达式。在高频段算法的耗散特性，用谱半径来说明更适合，一般用来度量算法对高频分量的耗散特性，特别地，当时，称算法具有高频渐进消去特性，即当算法计算一步以后高频极限完全地被耗散掉，而其它高频分量由高到低渐进地被耗散。由前面的叙述可以看到，算法放大矩阵的特征向量决定了算法对应特征方程的根，也就决定了算法的稳定性，同时确定了谱半径、以及算法耗散和弥散特性。这些特性有时也称算法的谱特性。放大矩阵相同的不同算法，称为互相相同的，放大矩阵不同但特征值相同，称算法互相相似，或称算法频谱等价的。对于算法的耗散特性，应该说明的是高频耗散特性对实际的结构动响应求解是有益的，因为实际结构进行有限元离散以后计算出的高

14、频行为并不真正代表系统的物理行为，它是结构系统在空间进行有限元离散的结果，是虚假的行为，而不具备高频耗散特性的算法，是将所有频率上的响应全部进行了积分，尽管步长取得相对较大时，高频的积分不准确，这样的计算结果显然与系统实际的反应不一致。但是同时要注意，它在低频也不同程度地引入了数值耗散，这样这些算法就不适合进行长时间的计算，因为长时间以后应该精确计算的低频响应，由于耗散特性的存在，已经被耗散得面目全非，因此，有耗散特性的直接积分方法只适合计算瞬态的、短时间内的低频动力响应。例5-3 分析Newmak族算法频谱特性解：对Newmak族算法来说，当时，算法才可能有二阶精度，我们仅讨论这一种情况。此

15、时算法放大矩阵的两个特征量为： (5-61)则谱半径不考虑物理阻尼时，对任意的有 (5-62) (5-63)就是无阻尼时，Newmak族算法不存在数值耗散，但有一定的相对周期误差。对有阻尼问题， (5-64)5.2.4 算法的超调特性谱半径这个指标对算法性能的影响还需要进一步说明的是，它只决定算法的长期特性，即可以保证随着算法计算步数的增加计算过程是数值稳定的。但对无条件稳定的算法，由于步长大小选择没有限制，一般在满足指定精度的条件下，尽可能取较大的时间步长，对于非零初始条件问题，在计算开始的几步可能会出现初始数据及其误差（如初始位移，速度的测量误差，初始加速度的计算误差）被放大的现象，这称为

16、超调(overshoot)。这种现象是放大矩阵A病态，有较大的条件数而产生的。实际应用时，由于当算法是收敛的，不会出现超调。一般为简单起见，只分析当时，在计算的第一步是否会出现超调。例54 分析Newmak平均加速度法的超调特性为分析简便起见，将Newmak平均加速度法用于无阻尼自由振动问题，此时其放大矩阵为： 在时，可得近似等式 (5-65)其中分别表示关于的零次和一次关系式。由此可知算法在位移上无超调，但由于初位移的影响，在速度上有关于线性超调现象。前面提过Wilson-方法有很强的超调现象，对无阻尼问题，从放大矩阵的各元素的表达式中，很容易得到在时，可得近似等式 (5-66) (5-67

17、)其中分别表示关于的二次和一次关系式。由此可知Wilson-方法在位移上关于初位移有二次超调，同时关于初始速度有一次超调。在速度上有关于初位移一次超调。另外，也可以用数值计算的方法对指定的初始条件，计算出近似解，然后与精确解比较，或计算系统能量范数： (5-68)然后将比较。显然，由于直接积分方法适合于短时间的瞬态问题计算，因此超调现象也是必须加以注意的。综上所述，对一个数值积分算法理论上要分析其相容性，稳定性，数值的耗散与弥散特性，对无条件稳定的算法还要分析其超调特性。这样才可能对算法的本质有深入的了解，进而指导数值计算结果的解释与分析。此外，由于直接积分方法对结构运动平衡方程进行数值积分的

18、目的在于估计结构真实的动力响应。为了精确地预计结构的动力响应，要求模态分解以后所有的单自由度平衡方程都必须被精确地积分，但在直接积分法中，对所有的方程积分都相当于采用相同的步长h，所以时间步长的选择必须针对系统的最小周期。如果是系统的最小周期的话，则h选为，其中一般n10。对条件稳定的算法，当然还要同时考虑这个选取是否满足算法稳定性的要求。对有条件稳定的算法，要求，若n取10，多数的条件稳定算法的稳定条件都满足，不难验证表51中的条件稳定算法全都满足。但需要注意的是，对于大型、复杂的实际结构，经过有限元离散以后通常都有上万，甚至几十万个自由度，其最大固有频率通常都很大，也就是系统的最小周期非常

19、小，此时，按来选取步长就非常小，这会大大增加计算量。而实际工程上只关心较低阶的固有频率，同时结构的响应也主要由若干较低阶的响应构成，因此在计算时高频可以不用精确积分，就积分出那些主要的，感兴趣的低频响应就可以了。也就是步长可选择为，比大倍。由于实际情况中可能会非常大，这样条件稳定算法的稳定条件就可能无法满足，而无条件稳定算法对步长的选取就没有稳定性的限制，因此对于实际的结构动力响应计算，多数都使用无条件稳定算法。5.3 矩阵特征值问题及解法531 .问题分类 这原本是广义特征值问题，但可以化为标准特征值问题，前乘，得 令，则有 这是代数里的标准的矩阵特征值问题，对结构动力学问题而言 M对角 对

20、称 M对称非对角 不一定对称 非对称矩阵的结构特征值问题计算量很大，此时作如下处理：1 对M进行Cholesky分解（因为M对称正定） L（下三角阵）（上三角阵）代入方程 前乘得 提出U得 令 ， 则有 可以征明与A有相同的特征值，因为它们是相似的矩阵，但是对称的因为 所以对称。 只不过特征向量有所变化，求出以后 至于Cholesky分解，可直接由M（或K）的元素计算得到 的各非零元素 i=j+1,j+2,.,n上式依次取j=1,2,.,n，即可求得L的下三角部分各列元素。又若记 则其下三角之i=1,2,.,n行元素，依次为 j=1,2,i-15.3.2 特征值、特征向量的一些特性1.对角阵、

21、三角阵、块对角阵、块三角阵对角阵和三角阵的特征值就是这些矩阵对角元素的数值。块对角阵和块三角阵的特征值就是这些矩阵的对角线上各个子块的特征值。2.几种特殊矩阵（1）实对称矩阵的特征值必为实数，其特征向量也可选为实向量。（2）反对称矩阵的特征值或者为纯虚数，或者为零。（3）正定对称矩阵的特征值全部大于零。（4）对称矩阵对应于不同特征值的特征向量彼此正交。3.与原矩阵A相关联的几种矩阵设矩阵A的特征值是，其对应的特征向量是x，则（1）阵的特征值是，特征向量是x；（2）阵的特征值是（），特征向量x；（3）是A的转置矩阵，则的特征值就是；（4）阵的特征值是，特征向量仍为x,m为整数；（5）若A非奇异，

22、则A的逆矩阵存在，记为，是的特征值为，特征向量为x；（6）若矩阵B与A相似，即有可逆阵P存在，使 则B的特征值也是，特征向量是4.特征值的积与积若矩阵A的特征值为、.，则有 它们可作为校核、估计甚至计算特征值的手段，读者可以自行验证。533 特征值问题的基本计算方法 由于高于四次的一元代权方程无法求精确解，所以必然要采用迭代方法求计算求近似解。目前的迭代方法，只不过迭代的方法和技巧不同。已有很多成熟方法。 方法选择时主要取决于 1) 所要求的特征对数目 2）K、M的阶，3）K、M的带宽以及是否带状。 根据用到的基本关系大致可分四类 1.向量迭代法（幂迭代法Power iteration met

23、hod） 基本关系 根据迭代格式不同又分： 正向迭代 （幂迭代）逆迭代 （反幂迭代）主要用于求特征向量。 2.变换法 基本关系： 质量归化为的模态矩阵 雅可比迭代（Jacobi）（小模型实对称阵标准特征值问题的全部特征对方法） 广义雅可比迭代（求全部特征对，小模型的广义特征值问题，或者有大量非对角线零元素和少量对角线零元素问题） 豪斯霍尔德（Householder）3.多项式迭代法（不单独使用）基本关系： 显式多项式迭代 隐式多项式迭代4.斯图姆（Sturm）序列法 利用的Sturm序列性质来求解。 实际工程结构的动力学问题绝大多数使用有限元法方法离散求解。随着计算机存储和速度能力的成倍提高，

24、所建模型越来越精确，几十万、上百万个自由度的计算成为可能。为此发展了一些针对大型特征值问题的方法，它们综合上述典型方法和技巧，常用的有：多项式迭代 行列式探索法（带宽较窄的低阶特征值问题。 子空间迭代法（Subspace iteration）(与行列式法相似，但计算量小) 兰索斯（Lanczos）法向量迭代 HQRI(HouseholderQR 迭代方法)（针对大型满阵的大多数或全部的特征值和向量）逆迭代方法任意给定初始的第一阶振型向量和第一阶特征值，一般常用的选取为 ，利用关系式按下列方式实施迭代1）选，2）令K1，2，.i) K正定利用K的乔莱斯基分解，进行线性代数代权方程组的求解，解出i

25、i)对实施对M的正定化，得新的下一次迭代向量 然后返回i)使进行下轮迭代，直至相邻的两次迭代值的差 norm雅可比（Jacbi）法雅可比（Jacbi）法是求解实对称矩阵全部特征值和特征向量的简单有效方法。该法自1864年问世以来，至今仍被广泛应用。 雅可比（Jacbi）法概念清楚，方法简单，易于编程，对低阶实对称矩阵求解效率较高，精度可以控制，所以在求高阶矩阵时将它和其它的求解方法联合使用。下面给出其计算步骤。（1） 记，给定精度指标，按式（5-47）和（5-48）算出（m=1,2.）（2） 巡视中的上三角非对角元素，若对所有，则转入（6），若对i=p,j=q，有，则记下p,q；（3） 按（4

26、-45）式计算及，；（4） 按（4-44）计算的各元素；（5） 按（4-46）式计算的各元素，将k加上1，返加（2）；（6） 若，则令，k+1，返回（2）；若，则迭代结束。 最后的对角线位置获得A的近似特征值，从中获得对应的特征向量。524 广义特征值问题解法简介子空间迭代法 当要求一个复杂系统的固有频率和固有振型时，常归结为解一个阶数很高的广义特征值问题，这是非常困难的。但工程中有用的是低频或某一频段内的固有频率和固有振型。此时，子空间迭代法是极其有效的，不但可以节省计算时间，而且所需的计算机内存也大为减少。下面先简单地介绍一下子空间迭代法的基本概念。 子空间迭代实质上就是对一组试验向量反复地使用里兹法和反迭代。此时可按下列步骤进行：（1） 形成系统的刚度矩阵K，质量矩阵M；（2） 形成初始迭代向量矩阵；（3） 求解；（4） 计算 （5） 求的特征值问题；（6） 计算