< < 第 17卷第 1期 　　　　　　　　　　　　　　湖 南 工 程 学 院 学 报 　　　　　　　　　　　　Vo1. 17. No. 1 2007年 3月 　　　　　　　　　　　　　Journal of Hunan Institute of Engineering　　　　　　　 　Mar. 2007 MATLAB在结构动力学中的应用 周后志 , 冷辉平 (湖南科技大学 土木工程学院 ,湖南 湘潭 411201) 　　摘 　要 : 通过两个实例介绍了 MATLAB 语言在结构动力学中的应用 ,通过结构的自振频率 、振型以及 动力响应在 MATALB 中的实现 ,说明了 MATLAB在结构动力学计算中的强大功能及其编程的便捷性 , 使科技人员真正地从繁杂的计算中解放出来. 关键词 : MATLAB;自振特性 ;动力响应 中图分类号 : TU311. 3　　文献标识码 : A　　文章编号 : 1671 - 119X (2007) 01 - 0091 - 04 　　 0　引 　言 在科学技术飞速发展的当今社会 ,计算机的应 用能力已成为评价科技人员综合能力的一项重要内 容. 在进行科学研究与工程技术应用的过程中 ,科技 人员往往会遇到大量繁重的数学运算和数值分析 , 一些传统的高级程序语言如 FORTRAN 等虽然能在 一定程度上减轻计算量 ,但它们要求应用人员要具 有较强的编程能力和对算法有深入的研究. 另外 ,在 运用这些高级程序语言进行计算结果的可视化分析 及图形处理方面 ,对非计算机专业的普通用户来说 , 存在着很大的难度. MATLAB 正是在这一应用要求 背景下产生的数学类科技应用软件. MATLAB 是 Matrix和 Laboratoty前三个字母的 缩写 ,是以矩阵计算为基础的程序设计语言 ,其指令 格式与教科书中的数学表达式非常相近 ,用 MAT LAB编写程序犹如在便笺上列写公式和求解 ,因而 被称为“便笺式 ”的编程语言 [ 1 ]. 同时 , MATLAB 具 有功能丰富和完备的数学函数库及工具箱 ,大量繁 杂的数学运算和分析可通过调用 MATLAB 函数直 接求解 ,大大提高了编程效率 ,其程序编译和执行速 度远远超过了传统的 FORTRAN语言 ,因而用 MAT LAB编写程序 ,往往可以达到事半功倍的效果. 在 图形处理方面 ,MATLAB可以给数据以二维 、三维乃 至四维的直观表现 ,并在图形色彩 、视角 、品性等方 面具有较强的渲染和控制能力 ,使科技人员对大量 原始数据的分析变得轻松和得心应手 [ 2 ] ,从根本上 满足了科技人员对工程数学计算的要求 ,将科技人 员及普通用户从繁重的数学运算中解放出来. 下面我们将演示运用 MATLAB 程序语言求解 结构动力学中结构的自振频域 、振型以及结构动力 响应的例子 ,为此我们先回顾一下有关结构的自振 特性以及动力响应的知识. 1　结构的自振特性和特征值 [ 3 ] 结构的自振特性是指结构的振动频率和振型 , 求结构的自振频率和振型也称对结构进行模态分 析 ,是结构动力学计算的主要内容之一. 计算经验指 出 ,结构的阻尼对结构的频率和振型的影响很小 ,所 以求频率振型时可以不考虑阻尼的影响. 此时系统 的自由振动方程如式 (1)所示 ,即 ¨ ( 1) [ K ] { x} + [M ] { x } = 0 当系统做自由振动时 ,各质点做简谐振动 ,各节 点的位移可表示为 } cosωt { x} = { 将 ( 2)代入 (1)式 ,并消去公因子 cosωt得到 ( 2) [ K ] -ω2 [M ] { [ K ] { } =ω2 [M ] { } = 0 或 ( 3) 因此求解式 ( 1)就是寻找式 ( 3)的 ω2 值和非零 } , 这种问题称为广义特征值问题. 记 λ = } 向量 { ω2 ,λ和 { }分别称为广义特征值和特征向量. 式 ( 3)可写成 ( [ K ] -λ[M ] ) { } = 0 ( 4) 收稿日期 : 2006 - 07 - 06 作者简介 :周后志 (1975 - ) ,男 ,硕士研究生 ,研究方向 :工程结构损伤诊断与识别. 

29 　　　湖南工程学院学报 　　　　　　　　　　　　　　　　2007年 这是一个齐次的线性方程组 , 若要有 { }的非 零解 ,系数行列式必须等于零 ,即 | [ K ] | -λ[M ] | = 0 (5) 展开此式可得 K11 -λM 11 K21 -λM 21 K12 -λM 12 … K1n -λM 1n K22 -λM 22 … K2n -λM 2n … … … … Kn1 -λM n1 Kn2 -λM n2 … Knn -λM nn = 0 如果弹性结构的总体刚度矩阵 [ K ]和总体质量矩阵 [M ]的阶数都是 n, 则上述行列式展开后为 λ的 n 次代数方程式 ,由此可求出 n个根 ,即 n个广义特征 值 λi , i = 1, 2, …, n, 从而求出结构的 n个自振频率 ωi = λi , ( i = 1, 2, …n) . 求得广义特征值 λi 后 , 就可利用式 ( 4 )算得对 i } , 它代表 n个质点的振幅构 应的广义特征向量 { 成的振型. 在弹性结构的振动分析中 , 结构的自由度总数 n往往很大 , 因此无法直接从上述代数方程中求解 广义特征值 λi ( i = 1, 2, …, n ) , 而大多采用变换法 和迭代法 ,如熟知的雅可比法 、正迭代和逆迭代法 等 , 采用 这些 方法 计算 的过程 比较 繁琐 , 而应 用 MATLAB程序语言则可使的求解结构的特征值的数 学计算简单的多. 2　用 MATLAB对结构自振频率、振型的分析 用 MATLAB 对结构的自振频率 、振型分析计算 的便捷性用实例来证明 : 例题 1 如图 1所示三层刚架结构 ,各层的楼面 质量分别为 m1 = 180 t, m2 = 270 t, m3 = 360 t;各层 的侧移刚度分别为 k1 = 98 MN /m , k2 = 196 MN /m , k3 = 294 MN /m ,求刚架的固有频率和振型. 运用 matlab语言编程 ,很容易就可得到 : %求结构的自振特性 (自振频率和振型 ) clear; %刚度 、质量输入 ,编号从下至上为 1、2、3… k0 (1) = 98000000; k0 ( 2) = 196000000; k0 ( 3) = 294000000; m0 ( 1 ) = 180000; m0 ( 2 ) = 270000; m0 ( 3 ) = 360000; n = 3; %层数 for i = 1: n; ka ( i, 1) = k0 ( i) ; %组成刚度数组 ; m ( i, i) =m0 ( i) ; %组成质量矩阵 ; end %形成总刚度矩阵 k ( n, n) = ka ( n) ; for i = 1: n - 1; k ( i, i) = ka ( i + 1) + ka ( i) ; end for i = 1: n - 1; k (i, i +1) = - ka (i +1); k (i +1, i) = - ka (i +1); end mn = inv (m ) w2 = eig (mn) ; %求特征值 　w = sqrt (w2) %角 k; 频率 3. 1415926) %频率 　T = 1. / f%周期 f =w / (2 %求振型 for i = 1: n; 　　　L = k - w2 ( i) m; L00 = L ( 2: n, 2: n) ; L01 =L (2: n, 1) ; 　　　X = - inv (L00) L01; xa (: , i) = X; end x1 = ones(1, n) ; x = [ x1′, xa′〗′ 运行程序得 : x = 1. 0000　　　1. 0000　　　1. 0000 0. 6485　　　 - 0. 6066　　　 - 2. 5419 0. 3018　　　 - 0. 6790　　　　2. 4396 此程序可以求得任意多层刚架结构的固有频率 和振型. 3　结构动力响应 图 1　三层刚架 求结构的动力响应 ,在数学上就是要求出运动 方程式 (6)的解答. 式 (6)是一个二阶常系数微分方 

第 1期 　　　　　　　　　　　　周后志等 :MATLAB 在结构动力学中的应用 39 程组 ,可以用数值积分的方法对方程直接求解 ,即按 时间增量 Δt逐步求解运动微分方程 ,直至反应终 了 ,这一方法称作逐步积分法. 逐步积分法既可用于 求解线性结构体系问题 ———在整个动力反应过程中 [ K ] , [M ] , [ C ]矩阵保持不变的问题 ;也可用于求 解非线性结构体系的问题 - [ K ] , [M ] , [ C ]矩阵随 动力反应的过程而变化的问题. 这里只讨论线性结 构体系的问题. 逐步积分法求解运动微分方程的基 本思路是 : (1)把连续的时间过程离散为 t1 , t2 , …, tn有限 个点 ,对于运动微分方程 } + [ K ] { x} = [ F ] · ¨ (6) [M ] { x } + [ C ] { x 只要求它们在上述每个时间离散点上得到满 足 ,也就是说 ,最终求解得到的只是位移 、速度和加 速度在有限个时间离散点上的值 ,而不是连续函数 · {X}、{ x ¨ } 、{ x }. (2)在每个时间间隔 △t内 ,假定位移 、速度和 加速度符合某一简单的关系. 而 △t的选择 ,要求保 证计算的稳定性与精确性. 从这样的基本思路出发 ,这里采用 W ilson - θ 法来求解结构的动力响应. W ilson -θ法的计算步骤 归纳如下 : ( 1)初始计算 1)形成总体刚度矩阵和总体质量矩阵 , 而阻尼 的影响一般不考虑. · 2)形成初始值 { x ( 0 ) } , { x ¨ ( 0 ) } , 并计算 { x (0) } 3)选取时间步长 △t和 θ(一般 θ= 1. 4) , 计算 积分常数 ,得 a0 = a3 = 6 (θ△t) 2 , a1 = θ△t a0 θ , a4 = 2 a2 θ, a6 = 1 - a5 = - 3 θ△t , a2 = 2α1 , 3 θ, a7 = △t 2 , a8 = △t2 6 4)形成有效刚度矩阵 ∧ [ K ] = [ K ] + a0 [M ] +α1 [ C ] ( 2)对于每一时间步长 ,循环计算 1)计算 t + △t时刻的有效载荷向量 ∧ F ] ( t +θ△t) = { F ( t) } +θ( { F ( t +θ△t) } - ( t) } ) · { F ( t) } ) + [M ] ( a0 { X ( t) } + a2 { x ¨ ( t) } + 2{ x ¨ ( t) } + a3 { x ( t) } ) · + [ C ] ( a1 { X ( t) } + 2{ x 2)求解 t +θ△t时刻的位移向量 { X ( t +θ△t) } ∧ [ K ∧ ] { X ( t +θ△t) } = F ( t +θ△t) } 3)计算 t + △t时刻的加速度 、速度 、位移向量 · { x ( t + △t) } = a4 ( { X ( t +θ△t) } - { X ( t) } ) + · a5 { x ¨ ( t) } + a6 { x · { x · ( t + △t) } = { x ( t) } ¨ ( t) } + a7 ( { ¨ ( t + △t) } + { x ( t) } ) · { X ( t + △t) } = { X ( t) } + △t{ x ¨ ( t) } + a8 ( { x ( t ¨ + △t) } + 2{ x ( t) } ) 4　用 MATLAB对结构动力响应的分析 例题 2　有一简单的结构系统 ,不考虑阻尼影 响 ,求此结构在 0～3. 5 s的位移响应. 已知结构初 始状态是静止的 ,受突加不变荷载作用. (见图 2) 图 2　二层刚架 运用 MATLAB 语言编程得 (其动力响应图见 图 3) : 图 3　结构动力响应 %求结构的动力响应 (wilson - θ计算结构动 力响应 ) theta = 1. 4; 　　% theta =θ deltat = 0. 28; % delta = t/ n < 0. 1T ( T:结构系统 

2 3 3 3 3 49 2 　　　湖南工程学院学报 　　　　　　　　　　　　　　　　2007年 自由振动最小周期 ,由上面的程序可得 T = 2. 8 s) %计算初始加速度值 y ( 0) ,静止 : v ( 0) = 0, x (0) = 0:运动方程 my + cv + kx = p ( t) y:加速度 , v: 速度 , x:位移 , p ( t) = F; x0 = [ 0; 0 ]; v0 = [ 0; 0 ]; y0 = [ 0; 10 ]; ylabel(‘幅度 ’) ; 运行程序得 : X1 = Columns 1 through 8 0. 0060 0. 0525 0. 1960 0. 4896 0. 9516 1. 5425 a0 = 6 / ( theta deltat) ^2; a1 = 3 / ( theta deltat) ; 2. 1623 2. 6702 a2 = 2 a1; a3 = theta deltat/2; a4 = a0 / theta; a5 = - a2 / theta; a6 =1 - 3 / theta; a7 = deltat/2; a8 = deltat^2 /6; K = k + a0 for i = 1: 14 F (1: 2, i) = [ 0; 10 ]; Y (1: 2, i) = [ 0; 0 ]; V ( 1: m; %有效刚度矩阵 2, i) = [ 0; 0 ]; X (1: 2, i) = [ 0; 0 ]; Columns 9 through 14 2. 9226　2. 8182 2. 3340 1. 5415 X2 = Columns 1 through 8 0. 3663　1. 3393 2. 6394 3. 9235 4. 8793 5. 3093 5. 1781 4. 6064 Columns 9 through 14 3. 8182　 3. 0605 2. 5233 2. 2862 deltat; t = ( i - 1) f1 = F ( 1: 2, i) + theta v0 + 2 ( a1 a2 = f ( t + theta y0) + c deltat) [ 0; 0 ] + m x0 + 2 v0 + a3 ( a0 x0 + y0) ; % f1 5　结束语 y0) ; deltat) v0 + a6 y0; % v0 + a8 ( y2 + 2 ( x1 - x0) + a5 f1; % x1 = x ( t + theta ( y2 + y0) ; % v2 = v ( t + deltat) x1 = inv ( K) y2 = a4 y2 = y ( t + deltat) v2 = v0 + a7 x2 = x0 + deltat % x2 = x ( t + deltat) Y (1, i) = y2 (1, 1) ; Y (2, i) = y2 (2, 1) ; V (1, i) = v2 (1, 1) ; V (2, i) = v2 (2, 1) ; X (1, i) = x2 (1, 1) ; X (2, i) = x2 (2, 1) ; F (1, i) = 0; F (2, i) = 10; y0 = y2; x0 = x2; v0 = v2; end t = 0: deltat: 3. 64; X1 =X (1,: ) X2 =X (2, : ) plot( t, X1, ’g’, t, X2, ’r’) ; axis( [ 0 5 - 15 15 ] ) ; title (‘结构动力位移响应 ’) ; xlabel(‘时间 ’) ; 从以上两个程序及其运行图我们可以清楚地看 到 ,运用 MATLAB 编程简单 ,大大提高了结构动力 学中振动问题的求解效率 ,而且计算结果可用图像 清楚直观地表达 ,效果良好. 用 MATLAB 编写程序 , 其程序编译和执行速度远远超过了传统高级语言 , 可以达到事半功倍的效果. 使科技人员对大量原始 数据的分析变得轻松自如和得心应手 ,真正将科技 人员从繁重的数学运算中解放了出来. 参 　考 　文 　献 [ 1 ] 　石博强. MATLAB数学计算范例教程 [M ]. 北京 :中国 铁道出版社 , 2004. [ 2 ] 　云舟工作室. MATLAB 教学建模基础教程 [M ]. 北京 : 人民邮电出版社 , 2001. [ 3 ] 　 (美 )克夫 ,彭 　津. 结构动力学 [M ]. 王光远 ,等译. 北 京 :科学出版社 , 1983. [ 4 ] 　侯新录.结构分析中的有限元法与程序设计 ———用 Visual C+ +实现 [M ].北京:中国建材工业出版社 , 2004. Applica tion of MATLAB in Structure Dynam ics ZHOU Hou - zhi, LENG Hui - p ing ( School of Civil Engineering, Hunan Science and Technology University, Xiangtan 411201, China) Abstract: This paper introducesMATLAB app lication in the structure dynam ics through two examp les, W ith the realization of the natural vibration frequency, vibration mode and dynam ic response byMATLAB , MATLAB power ful function and its convenience in the structure dynam ics calculation are illustrated. It can rely help peop le be free from the comp licated calculation. Key words: MATLAB; vibration characteristic; dynam ic response free