结构动力学方程常用数值解法

结构动力学方程常用数值解法对于一个实际结构，由有限元法离散化处理后，动力学方程可写为：...()MxCxKxFt从数学角度看，这是一个常系数的二阶线性常微分方程组，计算数学领域，常微分数值算法常用的有两大类：－、针对一阶微分方程数值积分法发展的欧拉法，中点法，Rugge-kutta（龙格—库塔）方法。二、直接基于二阶动力学方程发展的方法。对结构动力学问题的数值求解，常用的有两大类：一是坐标变换法，它是对结构动力方程式，在求解之前，进行模态坐标变换，实际上就是一种Rize变换，即把原物理空间的动力方程变换到模态空间中去求解。现在，普遍使用的方法是模态（振型）迭加法。二是直接积分法，它是对结构动力方程式在求解之前不进行坐标变换，直接进行数值积分计算。这种方法的特点是对时域进行离散，然后将该时刻的加速度和速度用相邻时刻的各位移线性组合而成。通常又称为逐步积分法。模态迭加方法，比较常用，但如下情况通常使用直接积分方法（即求解之前不进行模态分析）一、非比例阻尼，非线性情况。二、有冲击作用，激起高频模态，力作用持续时间较短，模态迭加计算量太大。一振型迭加法与Duhamel积分数值解按照有限单元法的一般规则,经过边界条件的约束处理,结构在强迫振动时多自由度体系的运动平衡方程可以表示为:MUCUKUR(1)其中,M是体系的质量矩阵,C是体系的阻尼矩阵,而K则是刚度矩阵.R为外荷载向量.U、U和U则分别是体系单元节点的位移、速度和加速度向量.上述动力平衡方程实质上是与加速度有关的惯性力MU和与速度有关的阻尼力CU及与位移有关的弹性力KU在时刻t与荷载的静力平衡。振型叠加法是把多自由度体系的结构的整体振动分解为与振型次数相对应的单自由度体系,求得各个单自由度体系的动力响应后,再进行叠加得出结构整体响应.振型叠加法原理是利用结构无阻尼自由振动的振型矩阵作为变换矩阵,将结构动力方程式(1)式变换成一组非耦合的微分方程.逐个地求解这些方程后,将解叠加即可得到动力方程的解。将体系单元节点的位移向量表示为如下的变换形式：()()UtXt（2）式中的变换矩阵是由动力方程对应的无阻尼自由振动方程解出的前m阶振型矩阵.即12[...]m;()Xt是与时间有关的m阶向量,X的各分量称为广义位移。将式（2）代入动力方程（1）并左乘以T，则可得广义位移为未知数的方程：()()()()MXtCXtKXtRt（3）式中TMM，TCC，TKK，TRR（4）现在进一步考察式（4）.考虑到特征向量的正交性,可得TMI，TK（5）于是对应于振型的广义位移的平衡方程(3)可改写为()()()()TTXtCXtXtRt（6）其中，为特征值212222......im（7）将式(2)稍加运算可得广义位移用有限元位移表示的形式TXMU（8）在(6)式中,当忽略了阻尼的影响,平衡方程为互不耦合的,可以对每个方程逐个地进行时间积分.出于相同的考虑,在对有阻尼的体系进行分析时仍然希望采用相同的计算过程去求解互不耦合的平衡方程式.问题是式(6)中的阻尼阵C通常不能象体系的质量阵和刚度阵那样由单元的刚度阵和质量阵装配而成.但当假定阻尼与固有频率成比例,即假定2TijiiijC（9）式中，i是振型阻尼参数;ij是Kronecker符号(当ij时，ij=1.当ij时,ij=0)。这时式(6)可简化为如下形式的若干个方程式2()2()()()iiiiiixtxtxtrt（10）其中()ixt的初始条件为下式00TitixMU，00TitixMU（11）式（10）表示了一个具有单位质量，刚度为2i的自由度体系当阻尼比为i时的运动平衡控制方程。这个平衡方程的求解可通过计算Duhamel积分求得。01()()sin(sincos)tiiiittiiiiiiiixtrtetdett（12）式中21iii（13）当利用式(9)来考虑阻尼的影响时意味着假设结构的总阻尼是每个振型的阻尼之和，而每个振型上的阻尼是能够量测的，况且在大多数情况下结构的阻尼比更易于量测。因而便于用来近似地反映结构体系的阻尼特性。同时在计算上也避免计算阻尼阵而只需计算刚度阵和质量阵。积分递推公式对以上方程式(10)，考虑某一模态的振动,并略去下标i可写为2()2()()()xtxtxtrt（14）在初始条件00ttxx，00ttxx（15）下的定解为0000()[sincos]ttxtettttx000011sin()sintttttettxretd（16）式中，21，将上式对时间求导，得02200()[sin]ttxtettx0000[cossin]ttettttx01()[sincos]tttrettd（17）将上式中的0t，t分别代以t，2t（其中2为时间步长），并按抛物线法则计算式中的卷积，有下式：2(2)[sin2cos2]()xtext221sin2()sin2()3extert4sin()3ertt（18）222(2)[sin2]()xtext2[cos2sin2]()ext2(cos2sin2)()3ert4(cossin)()(2)33ertrt（19）以上两式以矩阵表示为：11121314152122232425()()(2)()(2)()(2)xtxtaaaaaxtrtaaaaaxtrtrt（20）记1ee，221ee，0c，1cosc，1sins，2112ssc，22121cc（21）则上述A矩阵元素为112202()aeccs，1222/aes，131225aaa，1411254/aesa，150a22112aa，222202()aeccs，232225aaa241101254()aeccsa，25/3a应用上述递推公式，以前一时刻来求后一时刻的结果。计算不重复。当ija求出后，以后在时域中的步进求解只是一些简单的数组相乘。计算速度很快。二Newmark类方法下面讨论直接积分法234111()()()()()()26nnnnnxtxttxttxttxtOt（1）2311()()()()()2nnnnxtxttxttxtOt.......（2）21()()()()nnnxtxttxtOt.......（3）将（3）代入（1），（2）得：1()()()()nnnxtxtxtOtt.......22411()()()()()()36nnnnnttxtxttxtxtxtOt.....311()()()()()22nnnnttxtxtxtxtOt......1、可以直接略去高阶项2、用变权来调节2111[()]2nnnnnxxtxxxt11[(1)]nnnnxxxxt然后假设在1nt时刻近似满足运动方程1111nnnnMxCxKxF通过变换将速度和加速度用位移表示，代入运动方程，只剩n+1时刻位移一个未知数，得11nnKxQ参数不同选取包含着三个经典算法（1）Newmark平均加速度法12，14（2）Newmark线加速度法12，16（3）中心差分法12，0Newmark法的一般步骤：1初始值计算（1）形成系统矩阵K，M和C（2）定初始值0x，0x.和0x..（3）选择时间步长t，参数、。并计算积分常数:021at，1at3112a，41a，5(2)2ta，6(1)at，7at（4）形成等效刚度矩阵K01KKaMaC（5）K矩阵进行三角分解TKLDL2对每一时间步（1）计算tt时刻的等效载荷623145()()tttttttttQQMaxaxaxCaxaxax......（2）求解tt时刻的位移()TttttLDLxQ（3）计算tt时刻的加速度和速度023()tttttttxaxxaxax...........67ttttttxxaxaxWilson-法的一般步骤：1初始值计算（1）形成系统矩阵K，M和C（2）定初始值0x，.0x，..0x。（3）选择时间步长，并计算积分常数1.4，026()at，13at，212aa，32ta，04aa，25aa，631a，72ta，286ta。（4）形成等效刚度K01KKaMaC（5）将等效刚度进行三角分解TKLDL2对每一个时间步长（1）计算tt时刻的等效载荷......0213()(2)(2)ttttttttttttRQQQMaxaxxCaxxax（2）求解tt时刻的位移()TttttLDLxR（3）计算在tt时刻的加速度、速度和位移.....456()tttttttxaxxaxax7()ttttttxxaxx8()tttttttxxtxaxx三结构动力响应数值算法性能分析算法数值计算结果如何评价，针对不同的结构动力响应计算问题应该如何选择更合适的算法等是非常重要的问题。这就需要深入研究算法的数值计算性能，也就是算法的计算精度、稳定性等。对线性结构动力学问题，已经有证明对整个多自由度的积分，等价于将模态分解后对单自由度的积分的结果进行模态叠加，因此可以通过对单自由度问题的分析，来说明算法的特性，其中阻尼均假设为比例阻尼。算法用于结构动力学方程的有限差分表示为：)(22tfxxx以下算法的性能分析，均将算法用于这个方程。分别在相邻的不同时刻应用算法可得如下一般形式：kkkLAyy1其中A为放大矩阵或称逼近算子，kL为载荷逼近算子。TmkkkkTmkkkkxxxyxxxy,,,,1111对自由振动情况有0yAynn显然计算的第n步的值与A直接有关。例如，Newmak方法：dtAAA122221212thhAhh，222[1(12)]1(12)212(1)(1)dhhhAhh矩阵A的特征多项式为02)det(212AAIA)(212122111AAtraceAA,211222112detAAAAAA对Newmak方法有：211[1(21)()]24vAD,221[1(22)()]2AD其中h为时间步长，2,12hDWilson-方法，放大矩阵为：22366)1(13DAhh6)6)3(()6(22232)636)2/3