基于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