结构动力学数值算法.

结构动力学第五章结构动力学数值算法5.1结构动力学中常用数值方法...0.0()(0)(0)MxCxKxFtxxxv基本思想：首先给定待求时间长度T，[0,]tT，在其中取一系列离散点（i=0,1,…,n），我们不去求x(t)，只求()ixt即可，即给出待求响应在各离散时刻的近似值。两个离散时刻间隔称为步长，步长可以相同（等步长），也可以不同（变步长）。计算数学领域，常微分数值算法常用的有两大类1）针对一阶微分方程数值积分法发展的欧拉法，中点法，Rugge-kutta（龙格—库塔）方法。2）直接基于二阶动力学方程发展的方法对结构动力学问题的数值求解，常用的有两大类1）模态迭加2）直接积分模态迭加方法，比较常用，但如下情况通常使用直接积分方法（即求解之前不进行模态分析）i)非比例阻尼，非线性情况。ii)有冲击作用，激起高频模态，力作用持续时间较短，模态迭加计算量太大。Newmark类方法1）234111()()()()()()26nnnnnxtxttxttxttxtOt2）2311()()()()()2nnnnxtxttxttxtOt.......3）21()()()()nnnxtxttxtOt.......由3）得1()()()()nnnxtxtxtOtt.......代入1），2）得22411()()()()()()36nnnnnttxtxttxtxtxtOt.....311()()()()()22nnnnttxtxtxtxtOt......1）可以直接略去高阶项2）用变权来调节2111[()]2nnnnnxxtxxxt11[(1)]nnnnxxxxt然后假设在1nt时刻近似满足运动方程1111nnnnMxCxKxF通过变换将速度和加速度用位移表示，代入运动方程，只剩n+1时刻位移一个未知数，得11nnKxQ参数不同选取包含着三个经典算法1）12，14Newmark平均加速度法，梯形公式2）12，16Newmark线加速度法3）12，0中心差分法纽马克法的解题步骤初始值计算（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5.1.2威尔逊-法的解题步骤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时刻的等效载荷...02()(2)tttttttttRQQQMaxaxx...13(2)tttCaxxax（2）求解tt时刻的位移()TttttLDLxR（3）计算在tt时刻的加速度、速度和位移.....456()tttttttxaxxaxax7()ttttttxxaxx8()tttttttxxtxaxx5．2结构动力响应数值算法性能分析算法数值计算结果如何评价，针对不同的结构动力响应计算问题应该如何选择更合适的算法等是非常重要的问题。这就需要深入研究算法的数值计算性能，也就是算法的计算精度、稳定性等。对线性结构动力学问题，已经有证明对整个多自由度的积分，等价于将模态分解后对单自由度的积分的结果进行模态叠加，因此可以通过对单自由度问题的分析，来说明算法的特性，其中阻尼均假设为比例阻尼。5.2.1算法用于结构动力学方程的有限差分表示)(22tfxxx以下算法的性能分析，均将算法用于这个方程。分别在相邻的不同时刻应用算法可得如下一般形式kkkLAyy1A为放大矩阵或称逼近算子，kL为载荷逼近算子。TmkkkkTmkkkkxxxyxxxy,,,,1111有的也ky定义为Tkkhvx,或T2k,,kkahhvx对自由振动情况有0yAynn显然计算的第n步的值与A直接有关。例如，Newmak方法：dtAAA12221thAh2212hh，2221(12)2(1)dhAh[1(12)]12(1)hhh矩阵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/336()26(2232222322232hh)6(22D放大矩阵A的特征多项式为：02)det(32213AAAIA其中A1，A2，A3为该矩阵的三个特征向量，分别为矩阵的迹的一半、各阶主子式的和以及矩阵的行列式，分别表达如下：232222122223222222222323221863332(6)4361812(6)6633(6)AAA此外，在几个不同时刻应用数值算法，然后将方程中的速度和加速度项消去，可得数值算法关于位移的差分方程，例如Newmak方法，有2211(12)2[1(21)()]24nnxx211[1(22)()]02nx5．2．2算法的稳定性分析稳定性定义：设mii2,1,为放大矩阵A的特征值，则imax定义为A的谱半径，若特征值互异，则1的算法是稳定的，但若有重特征根，则要求1。如果算法的稳定性要求对步长的选取有限制，称算法是有条件稳定的，反之为无条件稳定的。判断方法：放大矩阵的谱半径小于等于1成立的充分条件是0102102122121AAAAA对33的放大矩阵0)2(10210323032302213132321321321321AAAAAAAAAAAAAAAA上两式是关于算法自由参数,的不等式，由它可以判断算法是否无条件稳定，若不是，将给出稳定条件。算法数值稳定性的物理解释：物理上，对一个无阻尼或者有阻尼自由振动系统，系统的能量随着时间不应该增加，有阻尼情况还应该减小。因此,一个数值方法的计算结果也不应该放大初始能量，如果经过若干步的数值计算以后，计算结果远比初始条件大，那就是数值算法本身计算是不稳定的。例5－1分析Newmak方法、Wilson-方法的稳定性解：将Newmak方法放大矩阵特征量代入稳定性分析表达式0)21(2201)21()2(2显然，当2,21，算法无条件稳定。对Wilson-方法有222232242606(21)012(166)0(416)24120(213)0容易看出，其中第一，二，五不等式恒成立，对第三，四不等式若希望对任意的均成立，则有：01640661232求解上述不等式得37.1231实际使用中通常选取=1.45.2.3算法的精度（相容性和收敛性）直接积分算法的相容性、收敛性分析同样要使用其位移型的差分方程，用精确解代替近似解，得到局部截断误差表达式，用符号)(kte表示)()()(1kkkktheLtAyty以最常用的线性三步法为例，局部截断误差用放大矩阵的特征量可表示为2321/)]2()()(2)([)(hhtxAhtxAtxAhtxtekkkkk然后将),(htxk)2(htxk在kt点进行泰勒展开，利用运动平衡方程化简即可。算法精度定义：若局部截断误差表达式为步长的s阶小量，则称算法是s阶相容的。相容加稳定等于收敛，其相容的阶数就是算法的精度阶。收敛性的含义是当时间步长趋于零，算法的数值解趋于精确解。例5－2分析Newmak法的相容性和精度解：其局部误差表达式得：221)]()(2)([)(hhtAtxAhtxtekkkk即将)(),(htxhtxkk在kt时刻点泰勒展开，并注意到在该时刻的运动方程有：hxtek)3(2)]21()612()21[()()(])412(1211212)612[(32)4(2hohx显然，当物理阻尼为零时，选择21算法是二阶的。物理阻尼的存在，使算法精度降了一阶。Newmak方法中有两个参数待定，每种特定的选取都是一个特定的算法，最常用的几个算法见表5-1表5-1：常用的Newmak族直接积分算法方法名称稳定条件无阻尼精度阶有阻尼精度阶类型1/21/4平均加速度方法（梯形法）无条件21隐式1/21/6线性加速度方法46.332c22隐式1/20中心差分方法2c21显式隐式、显式算法如果在一个时间步内需要求解一个隐式方程组，则称算法是隐式的，反之不需要求解方程，直接计算即可得到下一时刻的值，则称算法是显示的。Newmak方法是隐式的，但对于中心差分方法，若质量矩阵和阻尼矩阵都是对角矩阵就可以显示地计算。显然显示方法计算量要小得多。Wilson-方法的精度，不难分析，无论是无阻尼还是有阻尼其精度都是2阶，隐式方法。5.2.4算法耗散和弥散特性实际计算过程中，步长的选取可能不是很小，此时如何来度量算法的计算精度，当然可以针对有解析解的问题进行大量的数值计算，将数值解与解析解进行比较来分析算法的计算精度。理论上还可以通过数值耗散(disspation)和弥散（dispersion）来辅助度量与分析。以单自由度问题为例，该问题的解析解为：)sincos()(21tctcetxddt01xc，ddvc002，21d对于一个收敛的且有一定精度的算法，特征差分方程通常有一对共扼复根，hide)(2,1其中21d，该两根称为主根，其它根称为寄生根(spuriousroots)。解的一般形式可写为nimiindndtnctctcexn321)sincos(式中的称为算法阻尼比，称为算法频率，对应的/2T称为算法周