基于Matlab的Simulink结构动力仿真方法探究

闫欣宇，杨睿珊

(1.西南科技大学土木工程与建筑学院，四川绵阳621010；2.西南科技大学外国语学院，四川绵阳621010)

近些年来，地震频发，抗震工程中的结构可以通过结构动力学理论进行分析与计算，但其中需要解决庞杂的微分方程。Matlab是一个强大的矩阵计算数学软件，它是基于C语言开发，但是其编程又不需要繁琐的电脑编程语法知识，十分适合工科学生在进行研究时学习与参考，通过集中质量法对模型进行简化，便可以用Matlab Simulink仿真利用空间状态函数进行微分方程求解，通过振型分解的方法求得各个Simulink模块的系数进行仿真[1]。

1 结构动力学方程

Matlab是Matrix Laboratory的简称，意为是矩阵实验室。顾名思义Matlab是以计算矩阵和数组为核心软件，对矩阵迭代编程计算的功能相当强大。Simulink是Matlab中一个用来对动态系统进行建模、仿真和分析的软件包，它具有直观、简便、易于理解的特点。本文采用振型分解法计算其振型，频率，为了简单考虑采用集中质量模型。结构的运动微分方程为：

式中：M0，k0分别为原结构的质量矩阵、刚度矩阵；

在上式中：

其中采用瑞雷阻尼假定：C0=a1M+a2K

将上述的方程化简为:

在此方程左右两边同时乘以XT,经过化简可得

x(t)=Xq(t)

可以这样理解上式体系的位移可以看作是由各个振型分别乘以相应组合系数q1(t)，q2(t)后叠加而成。换句话说，这种方法实际上是位移按振型加以分解，所以求出广义坐标q的列向量，然后与各阶振型相乘即可得出所要求得的x的相应加速度位移坐标。

2 基于Matlab结构动力特性各编程求解

经过上述振型分解法的原理简介，Matlab式可以求得多自由度体系的各阶振型与频率，然后通过振型分解法的原理利用Matlab编写程序对结构进行分析求解，Matlab独有的eig函数可以快捷的求解出结构整体刚度矩阵k，质量矩阵m的特征值特征向量，通过整理和化简即是所需求的的振型和频率。代码为：

根据结构计算出层间刚度与集中质量编程处理

k=[1,1,1,1,1,1,1]*4.2e7; %层间刚度根据D值法可计算出各层刚度

kcj=k

mc=[1,1,1,1,1,1,1]*2.0e5；%集中质量法质量矩阵

m=diag(mc)；

cn=6；

O=zeros(cn)；

I=eye(cn)；

cjgdu

k0=cjgdu(kcj,cn)； %整体刚度矩阵cjgdu为整体刚度矩阵的求法，

[x,d]=eig(kcju,m)；%用Matlab自带的eig函数进行振型和自振频率求解

d1=sqrt(d)；

w=diag(sqrt(d))；

for i=1:7；

for j=1:7；

xgd(:,j)=x(:,j)/abs(x(7,j))；

end

鉴于人文社会科学的研究成果难以转化，且存在社会效益优于经济效益、长期效益优于短期效益的特点，本文研究的人文社会科学科研项目的人员费用主要包括工资费、劳务费、专家咨询费和绩效支出。截至目前的有关内容梳理如下：

end

a1=2*w(1)*w(2)*(0.05*w(2)-0.07*w(1))/(w(2)^2-w(1)^2);%瑞雷阻尼假定求a1,常系数

a2=2*(0.07*w(2)-0.05*w(1))/(w(2)^2-w(1)^2);

c0=a1*m+a2*k0

3 Simulink仿真模型

通过上述代码可以求得各个所需矩阵，然后利用Simulink仿真结构动力学微分方程进行模拟某7层钢筋混凝土框架。层间均高3 m，层间质量为2.0e5 kg，层间剪切刚度根据D值法计算可得4.2e7 N/m，用集中质量法简化参见图1。

图1 模型简化图

建立如图2所示的Simulink仿真模型。

在Matlab工作空间中用下列代码求出各个模块的参数

for i=1:7;

end

A=d’

图2 Simulink仿真模型

B=t

Ca=xgd

Cw=xgd

for i=1:cn;

h(:,i)=xgd(:,i);

end

for i=1:cn;

f1(i)=h(:,i)'*mc';

end

for i=1:cn;

f 2(i)=h(:,i)'*m*h(:,i);

end;

for i=1:cn;

gama(i)=f1(i)/f 2(i);

end

gama=gama’

for i=1:cn;

for j=1:cn;

xgd(:,j)=x(:,j)/abs(x(cn,j));

end

end

for i=1:6;

for j=1:6;

xgd(:,j)=x(:,j)/abs(x(6,j));

end

end

整体刚度矩阵的m函数

for i=1:cn-1

kcju(i,i)=kcj(i)+kcj(i+1);

kcju(i,i+1)=-kcj(i+1);

kcju(i+1,i)=-kcj(i+1);

kcju(cn,cn)=korc(cn)

end

为了方便与sap2000模型进行对比，采取正弦波时程分析，信号源为正弦函数波，从scope中输出顶层加速度的反应参见图3，峰值为1.16。

图3 加速度反应值

4 Sap2000模型分析

根据所示模型进行如下所示Sap2000建立的平面结构 模型，参见图4、图5。

图4 平面模型

图5 加速度反应值

其顶层加速度大小与Matlab建立模型大小基本相同，验证了结论正确性。

可见Matlab的Simulink模块在数值分析编程上的应用是可行的。

5 结论

(1)文中分别用上述两种算法进行求解，然后进行比较,Sap2000计算结果与Matlab编程求解的结果非常相似，考虑到结构体系的一些假设条件从而存在误差,因此,基于Matlab程序Simulink仿真模拟是一个十分可行的方法。

(2)Matlab是对于求解结构的受力特点，动力特性以及求解复杂的结构动力学微分方程有强大的计算能力。

(3)结构编程计算Matlab是一个具有友好界面的编程工具，适合科研计算工作者的学习。

[1] GB 50011-2001建筑抗震设计规范[S]

[2] 徐赵东.Matlab语言在抗震工程中的应用[M].北京:科学出版社,2004

[3] 徐斌.Matlab有限元结构动力学分析与工程应用[M].北京:清华大学出版社,2009