数值分析作业、MATLAB作业

基于MATLAB的结构动力特性分析

摘要：本文对结构动力特性分析的两种常用方法：矢量迭代法和子空间迭代法，进行了

简要概述，并在基于MATLAB的平台下，运用矢量迭代法和子空间迭代法以及MATLAB提供的函数eigs对平面61杆桁架实例进行简单的结构动力特性分析。

关键词：矢量迭代法 子空间迭代法 eigs 结构动力特性分析

1 引言

有限元发展至今天，已经成为工程数值分析的有力工具，在理论和实践上均取得了令人瞩目的成就，事实上它已发展成为工程领域中一门不可或缺的技术。本文是基于在当今工程和教育界非常流行的数学软件MATLAB来进行结构的动力学分析。采用MATLAB作为编程平台，利用MATLAB强大的科学计算和符号运算功能，跨越繁琐的公式推导和复杂的编程技巧，提高效率。

2 动力特性分析

2.1矢量迭代法

矢量迭代法是结构动力特性分析中经常使用的方法之一. 矢量迭代法分为逆迭代法和正迭代法两种. 逆迭代法向结构的最低频率和振型收敛，而正迭代法向结构的最高频率和振型收敛. 一般情况下，不使用正迭代法，这一方面是因为在实际工程问题中所感兴趣的是系统的少数低特征对，另一方面是因为以结构的有限单元离散化模型为基础计算结构的高特征对精度很差. 在具体实施矢量迭代法时，可以有各不相同的具体步骤，但它们的基本思想与使用的迭代关系式是相同的.

逆迭代法除了求最小特征值和特征向量外，和Gram-Schmidt(格拉姆-施密特)正交化过程相结合，可用来求取最低的几阶特征对. 对于第j阶特征对的求取，设有初始向量x?j0?及其第s?1,2,?时迭代公式为

Kx?js?1??Mx?js?

(01)

清型(即清理模型)是求取了最小特征对之后求其余特征对时所必须进行的一步，不仅包括特征对求取的迭代初始向量与前j?1阶特征向量正交，且每次迭代过程中都要进行正交化处理，不断把前j?1阶特征向量从迭代向量中清除掉，因为由于实际计算的误差，在迭代过程中不可避免地会产生低阶特征向量，迭代到最后还是得到最低阶的特征向量. 则Gram-Schmidt正交化过程如下 其中

xj?1?xj?1???i?i

i?1?0??0?j(02)

0? ?i??iTMx?j＋1(03)

初始向量选得好与坏，对迭代的效率影响极大. 显然，要求初始向量选得尽可能接近所求的振型. 一种简单的选取方法是，把结构简单地考虑为各个自由度是互不耦合的，即 (04) K?diag?k11k22?knn?

M?diag?m11m22?mnn?

(05) (06)

相应式(04)与式(05)的特征方程为非耦合形式 kii?i??imii?i ?i?1,2?,n,? 则有?i?kiimii,?i?1. 将?i从小到大排序，设?i排在第j位.

一般选取第一阶初始向量为

取第j?1阶初始向量为

Mx1?0???m11m22?mnn?

T(07)

12?i?Tn?Mx?j00?010?0??1??0 (08)

矢量逆迭代法的计算步骤如下：

选取p阶初始n维矢量y?j0?(j?1,2,?,p) 对于j?1,2,?,p进行循环 对于s?1,2,? 计算Kx?js?1??y?js? 计算y?js?1??Mx?js?1? 计算近似特征值 ??js?1??y??? ??????x?yx?js?1?s?1jTsjTs?1j如果特征值的精度满足条件

??js?1?－??js??j?s?1?≤?，则终止内循环，计算特征向量

?j?x?js?1???x?s?1?j?Ty?s?1?j?1/2.

j0?Gram-Schmidt正交化，xj?1?xj?1???i?i，其中?i??iTMx?j＋. 1?0??0?i?1否则计算下次迭代的初始向量yj?s?1??y?js?1???x?s?1?j?Ty?js?1??1/2.

2.2 子空间迭代法

子空间迭代法的基本思想是把逆迭代法和Ritz法结合起来，既利用Ritz法来缩减自由度，又在计算过程中利用逆迭代法使振型逐步趋近其精确值，由于它吸收了两种方法的优点，因而计算效果比较好. 经验表明，这是目前求解大型结构自振频率和振型的最有效方法之一. 对于n阶实对称矩阵M与K，迭代式可取为

? (09) KXk?1?MXk?是X关于M的正交归一化矩阵. X?的正交归一式中：Xk＋1是n?m阶矩阵，且m?n；Xkkk化可进行如下：借变换Xk，经Ritz缩聚后，有

Mk?XkTMXk，Kk?XkTKXk

(10)

分别记缩聚系统(Mk，Kk)的特征值矩阵与特征矢量矩阵为?k*与Pk，有

KkPk?MkPk?k*

(11)

式中，?k*是对角阵，Pk是关于Mk的正归一化矩阵.

故在下一步迭代中可取

??XP (12) Xkkk需要指出的是，在子空间迭代法中，初始向量的选取也将直接影响到迭代的收敛速度. 根据经验和实践，按如下方法选取是比较好的，即初始向量的第1列的元素全部置为1，以后各列按miikii中的大小排序，分别在该顺序对应的位置上置1mkk，而该列其余元素皆置0，其中k即是指第k个元素最大.

子空间迭代法的基本流程如下： 选定初始向量

???x(1)X1?1进行迭代运算

(1)(1)?x2?xm?

?，KX?Y Yk?MXk?1kk进行Ritz方法运算

??XP Xk＋1k＋1k＋1求解新的降阶的广义特征值问题

Kk＋1Pk＋1?Mk＋1Pk＋1?k＋1 其中，Mk＋1?XkT＋1MXk＋1，Kk＋1?XkT＋1KXk＋1.

如果特征值的精度满足要求，则终止循环. 计算原特征值问题的特征值和特征向量为

?i??*k＋1,i i ，ψi?XK＋1P2.3 MATLAB专用函数法

MATLAB提供了专门的函数来计算矩阵的特征值，即eig和eigs，它们既可用于求解标准特征值问题，也可用于求解广义特征值问题. 前者一般用来求解维数较小的全部特征值，而后者一般用来获取大型特征值问题的某几个特征对. 其调用格式如下：

E = eig(K, M) [V, D] = eig(K, M)；

[V, D] = eig(K, M, flag) flag ='chol' or 'qz'

[V, D]=eigs(K, M, nummode, SIGMA) SIGMA='SM' or 'LM'

3 结构动力特性实例分析

如图1所示：平面61杆桁架．

图1 平面61杆桁架

平面61杆桁架的材料特性为：弹性模量E?2.1?1011Pa，质量密度??7300kg/m3，所有单元的横截面积为1?10?4m2．分别采用3种方法计算结构的前9阶频率：矢量迭代法、

子空间迭代法和MATLAB提供的函数eigs．用MATLAB编程计算，其中主要用到的MATLAB函数有：

Inviter(k,m,p,epsilon)——用矢量迭代法计算结构的特征值和特征向量. Ssiter(k,m,p,epsilon)——用子空间迭代法计算结构的特征值和特征向量.

3.1 MATLAB程序

M文件如下：

%----------------------------------------------------------------------- % Plane 61-bar truss

%----------------------------------------------------------------------- % To compute the first nine eigenvalues of plane 61-bar truss %

E=2.1e11; A=1e-4;

density=7.3e3; node_number=26; element_number=61; nc=zeros(26,2);

nc(1,1)=0;nc(1,2)=0; nc(2,1)=0;nc(2,2)=1; for ig=1:12 for jg=1:2

nc(ig*2+jg,1)=nc(jg,1)+ig; %node_coordinate nc(ig*2+jg,2)=nc(jg,2); end end

en=zeros(61,2);

en(1,1)=1; en(1,2)=2; %en:element_node en(2,1)=1; en(2,2)=4; en(3,1)=2; en(3,2)=3; en(4,1)=1; en(4,2)=3; en(5,1)=2; en(5,2)=4; for i=1:11 for j=1:5

en(i*5+j,1)=en(j,1)+2*i; en(i*5+j,2)=en(j,2)+2*i; end end

en(61,1)=25; en(61,2)=26;

ed(1:node_number,1:2)=1; í:element_displacement constraint=[1,1;1,2;25,2];

for loopi=1:length(constraint);

ed(constraint(loopi,1),constraint(loopi,2))=0; end dof=0;

for loopi=1:node_number for loopj=1:2

if ed(loopi,loopj)~=0 dof=dof+1;

ed(loopi,loopj)=dof; end

end end

ek=E*A*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; %element_stiffness_matrix em=(density*A)/6*[2 0 1 0; ... %element_mass_matrix 0 2 0 1; ... 1 0 2 0; ... 0 1 0 2];

k(1:dof,1:dof)=0; %structural_stiffness_matrix

m=k; %structural_mass_matrix,same size with k

theta(1:61)=0; el(1:61)=0;

e2s(1:4)=0; % index of transform the element displament number to structure

for loopi=1:element_number for zi=1:2

e2s((zi-1)*2+1)=ed(en(loopi,zi),1); e2s((zi-1)*2+2)=ed(en(loopi,zi),2); end

el(loopi)=sqrt((nc(en(loopi,1),1)-nc(en(loopi,2),1))^2 +(nc(en(loopi,1),2)-nc(en(loopi,2),2))^2);

theta(loopi)=asin((nc(en(loopi,1),2)-nc(en(loopi,2),2))/el(loopi)); lmd=[cos(theta(loopi)) sin(theta(loopi)); -sin(theta(loopi)) cos(theta(loopi))];

t=[lmd zeros(2); zeros(2) lmd]; dk=t'*ek*t/el(loopi); dm=t'*em*t*el(loopi); for jx=1:4 for jy=1:4

if(e2s(jx)*e2s(jy)~=0)

k(e2s(jx),e2s(jy))=k(e2s(jx),e2s(jy))+dk(jx,jy); m(e2s(jx),e2s(jy))=m(e2s(jx),e2s(jy))+dm(jx,jy); end end end end

% Inverse iterative method p=9;

epsilon=1e-5;

[v,d]=Inviter(k,m,p,epsilon); %d:eigenvalue v:eigenvector frequency=sqrt(d)/(2*pi); % Suspace iterative method p=9;

epsilon=1e-5;

[v,d]=Ssiter(k,m,p,epsilon); frequency=sqrt(d)/(2*pi); % MATLAB function eigs

[v,d] =eigs(k,m,9,'SM');

frequency=sqrt(diag(d))/(2*pi); [frequency,indexf]=sort(frequency) v=v(:,indexf);

Factor=diag(v'*m*v);

v=v*inv(sqrt(diag(Factor))); % normalize eigenvector

3.2 结论

3种方法运行结果比较如下： frequency =

method 1 method 2 method 3 16.4815 16.4815 16.4815 54.9565 54.9564 54.9564 73.7466 73.7467 73.7467 132.1518 132.1518 132.1518 193.0646 193.0635 193.0635 222.2496 222.2514 222.2514 302.8284 302.8278 302.8278 337.6138 337.6155 337.6155 404.0057 404.0046 404.0042

4 总结

MATLAB具有强大的运算与分析能力。基于MATLAB平台下，大大的提高了结构的动力特性分析的效率。三种方法运行结果与理论分析结论相符，且都有较好的运算精度，完全满足工程需要。

参考文献

[1]柴鹏飞.结构力学的有限元法[M]. 北京：机械工业出版社，2007.203～206.

[2]朱艳英，陈月娥等.理论力学课程教学中的MATLAB应用研究[J].教学研究.2006.5：258～260.

[3]梁兰菊.MATLAB在力学教学中的应用.北京师范大学学报(自然科学版).2008年4月120～122.

[4]陈怀琛.MATLAB及其在理工课程中的应用指南[M].西安电子科技大学出版社.2007.160～163.转 [5] 蒲俊，吉家锋. MATLAB数学手册[M].上海:浦东电子出版社，2002.

本文来源：https://www.bwwdw.com/article/cl3h.html