（同济大学）结构动力学教程 第六章 结构动力学中常用的数值方法.pdf

第六章 结构动力学中常用的数值方法第六章 结构动力学中常用的数值方法 6.1 结构动力反应的数值解法（结构动力反应的数值解法（step-by-step法） 一般结构动力学方程： 初始条件： 法） 一般结构动力学方程： 初始条件： 6.6.1 中心差分法中心差分法 位移对时间导数的差分表示： 代入得近似式： 引入记号： 简化上式为 位移对时间导数的差分表示： 代入得近似式： 引入记号： 简化上式为 qxkxcxm=+? ? ; 00 00 xxxx tt ?= )2( )( 1 )( 2 1 2 tttttt ttttt xxx t x xx t x + + + = = ? ? ? tttttt x t c t mx t mkqx t c t m + = + 2 1 )( 1 )( 2 2 1 )( 1 222 021 2 0 2, 2 1 , )( 1 aa t a t a= = = tttt qxkxcxm=+? ? ttttt ttttt xtxtxx xtxtxx 2 1 2 1 2 2 ? ? ? ? + + + ttttt xcamaxmakqq camam = += )()( 102 10 式中， 根据及 式中， 根据及t时刻的向量时刻的向量x的值时刻的值时刻x的方程 求解时刻的向量 的方程 求解时刻的向量xt时刻的向量 的计算：根据初始条件 及 其中，可根据算得。 时刻的向量 的计算：根据初始条件 及 其中，可根据算得。 稳定条件稳定条件：时间步长 （为体系最小自振周期） ：时间步长 （为体系最小自振周期） ttttt xcamaxmakqq camam = += )()( 102 10 tt tt+ ttt qxm= + tt+ xx? ?， t x 00 xx?，)( 2 1 0tt xx t x =? )2( )( 1 0 2 0tt xxx t x + =? ? 0 0 0 1 0 2 1 2 1 x a x a xx t ? ?+= 0 x ? ? 0000 xkxcqxm=? ? n cr t tt= n t ttt qxm= + tt x )2( )( 1 )( 2 1 2 tttttt ttttt xxx t x xx t x + + + = = ? ? ? 计算步骤计算步骤： 1，初始值计算： （ ，初始值计算： （1） 矩阵 （ ） 矩阵 （2） 初始值，确定时间步长 （ ） 初始值，确定时间步长 （3） 计算 （ ） 计算 （4） 计算 （ ） 计算 （5） 计算） 计算 2，对每一时间步 （ ，对每一时间步 （1） 计算： （ ） 计算： （2） 求解位移向量： （ ） 求解位移向量： （3） 求解速度、加速度向量） 求解速度、加速度向量 ,kcm 000 , ,qxx? t 021 2 0 2, 2 1 , )( 1 aa t a t a= = = 0 0 0 1 0 2 1 2 1 x a x a xx t ? ?+= 10 camam+= ttttt xcamaxmakqq =)()( 102 ttt qxm= + )2( )( 0 1 tttttt ttttt xxxax xxax + + += = ? ? ? 0000 xkxcqxm=? ? 6.6.2 newmark法法 设时间间隔内加速度呈线性变化，取设时间间隔内加速度呈线性变化，取常数常数： 若取（平均加速度） 取 ： 若取（平均加速度） 取速度速度的一阶泰勒展开： 将代入 用同样方法处理 的一阶泰勒展开： 将代入 用同样方法处理位移位移 泰勒展开： 类似地设时间间隔内： 泰勒展开： 类似地设时间间隔内： ttt+ ) 10()(+= + tttt xxxx? ? ? ? ? 2/1=)(2/1 ttt xxx? ? ? ?= + tt x + ?txxx ttt += + ? ? )( tttt xxxx? ? ? ? ?+= + txtxxx tttttt += + )1 (? ? ? 2 2 1 txtxxx tttt += + ? ? ttt+)5 . 00()(2+= + tttt xxxx? ? ? ? ? 22 )5 . 0(txtxtxxx ttttttt += + ? ? ? t tt+ tt x + ? ? x ? ? t x ? ? 根据时刻的运动方程： 及 根据时刻的运动方程： 及 三组方程联立三组方程联立，可根据算得 （ ，可根据算得 （1）由用表示： （ ）由用表示： （2）代入用表示 （ ）代入用表示 （3）代入解得： 式中， ）代入解得： 式中， 无条件稳定条件无条件稳定条件： tt+ tttttttt qxkxcxm + =+? ? txtxxx tttttt += + )1 (? ? ? 22 )5 . 0(txtxtxxx ttttttt += + ? ? ? ttt xxx? ? ， tttttt xxx + ? ?， 22 )5 . 0(txtxtxxx ttttttt += + ? ? ? tt x + ? ? tt x + ttttttt xtx t xx t x)1 2 1 ( 1 )( 1 2 ? ? ? = + txtxxx tttttt += + )1 (? ? ? tt x + ? tt x + tttttttt qxkxcxm + =+? ? tt x + tttt qxk + = + + + + += + += + ttt tttttt x t xx t c xx t x t mqq c t m t kk )2( 2 )1( )1 2 1 ( 1 1 1 2 2 2 ? ? ? ? 2 )1 ( 4 1 2 1 +； 计算步骤计算步骤： 1，初始值计算： （ ，初始值计算： （1） 矩阵 （ ） 矩阵 （2） 初始值确定时间步长及参数） 初始值确定时间步长及参数( ) （3） 计算 （ ） 计算 （4） 形成等效刚度矩阵：） 形成等效刚度矩阵： 2，对每一时间步 （ ，对每一时间步 （1） 计算等效荷载： （ ） 计算等效荷载： （2） 求解位移向量： （ ） 求解位移向量： （3） 求解速度、加速度向量：） 求解速度、加速度向量： ,kcm 000 , ,xxx? ? t, 2 )1 ( 4 1 2 1 +； tata t aa a t a t a t a = = = = = = 7654 321 2 0 ),1 (),2( 2 , 1 1 2 1 , 1 , 1 10 camakk+= )( )( 541 326 ttt tttttt xaxaxac xaxaxamqq ? ? ? ? + += + tttt qxk + = ttttttt tttttt xaxaxxax xaxaxx )( 320 76 ? ? ? ? ? ? = += + + 6.6.3 wilson-法法 线性加速度法的推广。设时间间隔内加速度呈 线性变化。增量时刻加速度： 线性加速度法的推广。设时间间隔内加速度呈 线性变化。增量时刻加速度： 加速度变化率加速度变化率： 无条件稳定条件：（为线性加速度法） 速度与位移泰勒展开： 设法求得先求 在、即时刻，有 ： 无条件稳定条件：（为线性加速度法） 速度与位移泰勒展开： 设法求得先求 在、即时刻，有 ttt+ )0(t )()( 1 tttttttt xx t xxxx t x? ? ? ? ? ? ? ? ? += = + 常数 420815. 137. 1=，优化值：1= )( 6 2 1 )( 2 3 2 2 ttttttt tttttt xx t xxxx xx t xxx ? ? ? ? ? ? ? ? += += + + t= tt+ )2( 6 )( 2 22 ttttttt tttttt xx t xtxx xx t xx ? ? ? ? ? ? + += + += + + ttt+ tt x + ? ? t x? ?tt x + ? ? tt+ += + ttxtxxx tttt 2 2 1 ? ? ? ? ) 6 1 2 1 32 += + ttxtxtxxx ttttt ? ? ? ? tt x + ? ? tt x + 由此求得时刻加速度和速度的表达式： 将代入时刻的运动方程， 由此求得时刻加速度和速度的表达式： 将代入时刻的运动方程，荷载也取线性投影荷载也取线性投影： 以位移为未知量建立求解方程，即： 式中， ： 以位移为未知量建立求解方程，即： 式中， ttttttt ttttttt x t xxx t x xx t xx t x 2 2)( 3 2 6 )( 6 22 ? ? ? ? ? = = + + tttt xx + ? ?， tt+ ）（)( tttttt tttttttt qqqq qxkxcxm += =+ + + ? ? tttt rxk + = + + + + += + += + ttt ttttttttt x t xx t c xx t x t mqqqr c t m t kk 2 2 3 2 6 )( 6 )( 3 )( 1 2 2 ? ? ? ? tt x + tt+ 解得后，代入 求得。 然后求时刻的解。 令，将代入： 得： 解得后，代入 求得。 然后求时刻的解。 令，将代入： 得： tt x + ttttttt xx t xx t x2 6 )( 6 22 ? ? ? = + tt x + ? ? t= tt x + ? ? )( 6 2 1 )( 2 )( 3 2 2 ttttttt tttttt ttttt xx t xxxx xx t xxx xx t xx ? ? ? ? ? ? ? ? ? ? ? ? ? += += += + + + )2( 6 )( 2 ) 3 1 ( )( 6 )( )( 6 1 ) 1 1 ( 2 2 ttttttt tttttt tttttttttt xx t xtxx xx t xx xx t xx t xxx ? ? ? ? ? ? ? ? ? ? ? + += + += + =+= + + + tt+ 计算步骤：计算步骤： 1，初始值计算： （ ，初始值计算： （1） 矩阵 （ ） 矩阵 （2） 初始值，确定时间步长及参数 （ ） 初始值，确定时间步长及参数 （3） 计算 （ ） 计算 （4） 形成等效刚度矩阵：） 形成等效刚度矩阵： 2，对每一时间步 （ ，对每一时间步 （1） 计算等效荷载： （ ） 计算等效荷载： （2） 求解位移向量： （ ） 求解位移向量： （3） 求解加速度、速度、位移向量：） 求解加速度、速度、位移向量： ,kcm 000 , ,xxx? ? t4 . 1= 6 , 2 , 3 1, 2 ,2, 3 , )( 6 2 876 2 5 0 4 3121 2 0 t a t aa a a a a t aaa t a t a = = = = = = 10 camakk+= tttttttttttt xaxxacxxaxamqqqr22)( 3120 ? ? ?+= + tttt rxk + = )2( )( )( 8 7 654 ttttttt tttttt tttttt xxaxtxx xxaxx xaxaxxax ? ? ? ? ? ? ? ? ? += += += + + + 6.2 矩阵特征值问题及解法矩阵特征值问题及解法 6.2.1 化广义特征值问题为标准特征值问题化广义特征值问题为标准特征值问题 广义特征值广义特征值问题：问题： 标准特征值标准特征值问题： 由于不一定是 问题： 由于不一定是对称矩阵对称矩阵，使问题复杂化 解决方法：对作 ，使问题复杂化 解决方法：对作cholesky分解分解： （为对角元均不为零的下三角阵）标准特征值问题： 是 ： （为对角元均不为零的下三角阵）标准特征值问题： 是对称矩阵对称矩阵。 212 = kmmk 0)(=xia , 21 = xkma 1 kma = m t llm= l 0)(=xia , 21 xllkla tt = a )()( 21211 2 22 xxllxlkll lxllllk llmk t tttt t = = = 左乘 令： xxa= 的的cholesky分解：需计算： 先由 分解：需计算： 先由m计算计算l。依次取，计算下三角。依次取，计算下三角各列各列元素 的各元素： 再由 元素 的各元素： 再由l计算。依次取，计算下三角计算。依次取，计算下三角各行各行元素 的各元素： 则： 元素 的各元素： 则： m nj, 2 , 1 = = nnnn lll ll l l 21 2221 11 njji lllml lml jj j r jririjij j r jrjjjj , 2, 1 1 1 2/1 1 1 2 += = = = = 1 lni, 2 , 1 = = nnnn vvv vv v l 21 2221 11 1 1, 2 , 1 /1 1 1 = = = = ij lvlv lv ii i r rjirij iiii t ll 1， tt ll）（ 1 = 6.2.2 特征值、特征向量的一些特性（不加证明）特征值、特征向量的一些特性（不加证明） 1 对角阵、三角阵对角阵、三角阵、块对角阵、块三角阵、块对角阵、块三角阵 *1） 对角阵、三角阵的特征值即为对角元素） 对角阵、三角阵的特征值即为对角元素 2）块对角阵、块三角阵的特征值即为对角线上子块的特征值）块对角阵、块三角阵的特征值即为对角线上子块的特征值 2 特殊矩阵特殊矩阵 *1） 实对称矩阵的特征值必为实数，特征向量为实向量） 实对称矩阵的特征值必为实数，特征向量为实向量 2） 反对称矩阵的特征值或为纯虚数、或为零） 反对称矩阵的特征值或为纯虚数、或为零 *3） 正定对称矩阵的特征值全部大于零） 正定对称矩阵的特征值全部大于零 *4） 对称矩阵对应于不同特征值的特征向量彼此正交） 对称矩阵对应于不同特征值的特征向量彼此正交 *3与原矩阵与原矩阵a相关联的矩阵相关联的矩阵 设矩阵设矩阵a的特征值为、对应的特征向量为的特征值为、对应的特征向量为 1）的特征值为、对应的特征向量为）的特征值为、对应的特征向量为 x ia+ x 2）的特征值为、对应的特征向量为）的特征值为、对应的特征向量为 3）的特征值就是）的特征值就是 4）的特征值为、对应的特征向量仍为）的特征值为、对应的特征向量仍为 *5）若非奇异，则的逆矩阵存在为，其特征值 为、 特征向量也为 ）若非奇异，则的逆矩阵存在为，其特征值 为、 特征向量也为 *6）若与相似、即有可逆矩阵存在使， 则的特征值也为、特征向量为 ）若与相似、即有可逆矩阵存在使， 则的特征值也为、特征向量为 *4. 特征值的和与积 设矩阵的特征值为，则有（供校核用） 特征值的和与积 设矩阵的特征值为，则有（供校核用） m a m x a a 1 a /1 x b a p 1 papb = b 1 xp a n , 21 det 21 1 21 a a n n i iin = =+ + = a x t a 5. 特征向量规范化 设矩阵的特征向量为 ，则亦为的特征向量。 常将各阶特征向量规范化标准特征向量，方法： 特征向量规范化 设矩阵的特征向量为 ，则亦为的特征向量。 常将各阶特征向量规范化标准特征向量，方法： 1） 最大值归一法最大值归一法 对第对第r阶特征向量找 第 个元素最大为 阶特征向量找 第 个元素最大为1 2）特殊点归一法特殊点归一法 对第对第r阶特征向量某一点感兴趣。 第 个元素为 阶特征向量某一点感兴趣。 第 个元素为1、但不一定最大、但不一定最大 3）模态质量归一法模态质量归一法 广义第广义第r阶特征向量（模态）。第阶特征向量（模态）。第r阶模态质量为 为使，将 使以下关系成立： 阶模态质量为 为使，将 使以下关系成立： a xxa )( r x = )()( max rr i xx )( )( )( r r i r x x x= i )( r x )( )( )( r r i r x x x= i )( r )()( r t r r mm= 1= r m )()( r r r m = 2 rr k= 6.2.3 标准特征值问题解法简介标准特征值问题解法简介 标准特征值标准特征值问题：问题： 求全部特征值问题矩阵变换方法：求全部特征值问题矩阵变换方法：jacobi法法，广义，广义jacobi法法 求部分特征值问题向量迭代方法求部分特征值问题向量迭代方法：幂迭代法，反幂迭代法，：幂迭代法，反幂迭代法， 子空间迭代法子空间迭代法 jacobi法法 标准特征值问题：，为实对称矩阵。 可通过 标准特征值问题：，为实对称矩阵。 可通过正交矩阵作相似变换化为对角矩阵正交矩阵作相似变换化为对角矩阵： 则的第 ： 则的第n个对角元素即为个对角元素即为a的第的第n 个特征值， 而的第列即为第阶特征向量。 个特征值， 而的第列即为第阶特征向量。 0)(=xia xxa=a s ),( 21n t ddddiagdsas = d si i i 正交矩阵实际上是坐标变换矩阵，正交矩阵实际上是坐标变换矩阵，通过坐标系的旋转确定通过坐标系的旋转确定 的 个正交特征向量的方向。但正交矩阵不容易找。的 个正交特征向量的方向。但正交矩阵不容易找。 jacobi法的基本思想法的基本思想是通过是通过多次多次“小的旋转小的旋转” 逐渐实现最终旋转逐渐实现最终旋转。 “小的旋转小的旋转”的生成法则： 设经过次 的生成法则： 设经过次 “小的旋转小的旋转”后为矩阵。 其 后为矩阵。 其绝对值最大的非对角元素绝对值最大的非对角元素为，则可取如下形式：为，则可取如下形式： s ans sas t i s a1i i a i pq a i s );,(0 ),(1 sin cos jiqpjis qpis ss ss ij ii qppq qqpp = = = = i s = 1 1 cossin 1 1 sincos 1 1 ? ? ? ? ? ? ? ? ? ? i s 列p列q 行q 行p 经变换后得到的新矩阵为 i s ii t ii sasa 1 = + pq a 展开后的元素为 其中， 为 展开后的元素为 其中， 为旋转角旋转角，根据的条件确定。 由上述第 ，根据的条件确定。 由上述第5式可得： 重复运用上述变换，可使一实对称矩阵按任意精度变换成 对角矩阵。 式可得： 重复运用上述变换，可使一实对称矩阵按任意精度变换成 对角矩阵。总的变换矩阵为各总的变换矩阵为各“小的旋转小的旋转” 变换矩阵之积变换矩阵之积 （其他元素） i mn i mn i pq i qq i pp i pq i pp i pq i pp i qq i qq i pq i pp i pp i qj i pj i qj i qj i pj i pj aa aaaa aaaa aaaa qpjaaa qpjaaa = += += += = = + + + + + + 1 221 221 221 1 1 )sin(coscossin)( coscossin2sin sincossin2cos ,cossin ,sincos 0 11 = +i qp i pq aa 4422 1 = i pq i pp i qq a aa arcctg i s n ssss 21 = 1 +i a 例题：设，用例题：设，用jacobi法求解特征值。 任取为第一次旋转准备消去的元素 法求解特征值。 任取为第一次旋转准备消去的元素p=1,q=2。得 按求得 故 。得 按求得 故 = 111 111 111 a 12 a = 100 0cossin 0sincos 1 s = i pq i pp i qq a aa arcctg 22 1 42 11 2 1 = =arcctg = 100 02/12/1 02/12/1 1 s = 120 220 000 )( 1112 sasa t );,(0 ),(1 sin cos jiqpjis qpis ss ss ij ii qppq qqpp = = = = 再依法求得 故 已实现所有非对角元素均为零，即需求的对角矩阵。 相应阵为 一般为一无限迭代过程，计算到 再依法求得 故 已实现所有非对角元素均为零，即需求的对角矩阵。 相应阵为 一般为一无限迭代过程，计算到全部非对角元素平方和 小于某一指定值结束 全部非对角元素平方和 小于某一指定值结束。 = 3/23/10 3/13/20 000 2 s 000 030 000 )( 2223 dsasa t = d s = 3/23/10 6/13/12/1 6/13/12/1 21 sss 6.2.4 广义特征值问题解法简介（广义特征值问题解法简介（子空间迭代法子空间迭代法） jacobi法仅适用于求解法仅适用于求解低阶低阶矩阵特征值问题。难于应用于复 杂系统。同时，工程中有用的往往仅是 矩阵特征值问题。难于应用于复 杂系统。同时，工程中有用的往往仅是低频或某一频段低频或某一频段的固有 频率和振型。此时 的固有 频率和振型。此时子空间迭代法子空间迭代法十分有效。 对 十分有效。 对广义特征值问题广义特征值问题：有：有n个线性无关的特征向量 它们关于正交，张成 个线性无关的特征向量 它们关于正交，张成n维向量空间。 试猜测该系统 维向量空间。 试猜测该系统q个个n维向量维向量， 在其中张成 ， 在其中张成q维子空间维子空间 子空间迭代法就是子空间迭代法就是通过反复迭代，使通过反复迭代，使q个试验向量低阶振型 分量放大，靠拢所张成的子空间 个试验向量低阶振型 分量放大，靠拢所张成的子空间。 即通过迭代 。 即通过迭代收敛于收敛于q个低阶的特征向量个低阶的特征向量。 必须注意的是，在迭代过程中需进行 。 必须注意的是，在迭代过程中需进行关于的正交处理关于的正交处理，即 使不断地旋转，最后指向 ，即 使不断地旋转，最后指向q个特征向量的方向个特征向量的方向 2 mk= )()2()1( , , n aaa m )()2()1( , , q xxx )()2()1( , , q aaa m )()2()1 ( , , q xxx 步骤：求解步骤：求解n阶广义特征值问题： （ 阶广义特征值问题： （1）假定前）假定前q阶振型，阶振型，q应略大于需求的特征向量的阶数应略大于需求的特征向量的阶数p , 推荐 将 推荐 将q阶振型自左至右排列，构成的矩阵 （ 阶振型自左至右排列，构成的矩阵 （2）求解。对作反迭代改善初始）求解。对作反迭代改善初始q阶振型 的品质（增强低阶振型分量）。 首先从开始，求解 阶振型 的品质（增强低阶振型分量）。 首先从开始，求解q个方程组可算得。 （ 个方程组可算得。 （3）计算 （ ）计算 （4）求的特征值问题（降为）求的特征值问题（降为q阶阶,jacobi法 等）。 （ 法 等）。 （5）计算。计算新的迭代向量并进行正交规 范化。判断相邻两次频率差满足精度要求，若满足结束迭代、 否则继续第（ ）计算。计算新的迭代向量并进行正交规 范化。判断相邻两次频率差满足精度要求，若满足结束迭代、 否则继续第（2）步进行新一轮计算直至收敛。）步进行新一轮计算直至收敛。 2 xmxk= )8,2min(+=ppq qn 0 qn x 12 = kk xmxk 0 x 0 1xk= 1 x ktkktk xmxmxkxk)()( * = kk qq k qq amak 2* = kkk axx= )qq （ * 2 amak= 初始向量的选取初始向量的选取直接影响迭代收敛速度，建议如下选取初始 向量： （ 直接影响迭代收敛速度，建议如下选取初始 向量： （1） 第一列全部置） 第一列全部置1。 （ 。 （2） 以后各列按中的大小顺序，分别在该顺序的 位置上置（即第个元素最大），其余元素置 ） 以后各列按中的大小顺序，分别在该顺序的 位置上置（即第个元素最大），其余元素置0 例题例题： （ 表示可取任意数） 初始向量： ： （ 表示可取任意数） 初始向量： 2-4（1/6） 3-3（1/10） 4-1（1/8） iiii km / kk m/1k k = = 4 6 10 9 8 ; 1 1 2 3 2 mk = 4 6 5 3 4 / iiii km = 4/10001 0006/11 0010/101 00001 08/1001 0 x 6.3 结构固有频率的结构固有频率的rayleigh近似解法近似解法 广义特征值广义特征值问题： 设及为第 问题： 设及为第r个特征值及特征向量，有 用乘两边，得 最大势能 个特征值及特征向量，有 用乘两边，得 最大势能=最大动能 设为 最大动能 设为任意向量任意向量，则，则rayleigh商商定义为（为的函数） 若，则 定义为（为的函数） 若，则rayleigh商即为对应于该特征向量的特征值商即为对应于该特征向量的特征值 2 qmqmqk= r rq rrr qmqk= t r q 2 1 r r r t r r t r rrr t rrr t r m k qmq qkq qmqqkq= 2 1 2 1 2 q )( qmq qkq r t t =q q rqq = r t rr qkqv 2 1 max, = r t rrr qmqt 2 1 max, = 设任意向量可表示成特征向量的线性组合 并设特征向量已规格化 将代入 设任意向量可表示成特征向量的线性组合 并设特征向量已规格化 将代入rayleigh商商 考虑与第考虑与第r阶特征向量相似的情况阶特征向量相似的情况 线性组合的系数可写成：线性组合的系数可写成： cqqcqcqcqcq n i iinn =+ += =1 2211 iqmq qkq t t = = cqq = = = = n i i i n i i t t tt tt t t c c cc cc cqmqc cqkqc qmq qkq r 1 2 2 1 )( q r q q ) 1(;, 2 , 1 = irii rinicc nrnrrrrrrrrrr qcqcqcqcqcqcq 11112211 + + += + rinici =;, 2 , 11 将 代入 将 代入rayleigh商 并除