閆欣宇，楊睿珊

(1.西南科技大學(xué)土木工程與建筑學(xué)院，四川綿陽621010；2.西南科技大學(xué)外國語學(xué)院，四川綿陽621010)

近些年來，地震頻發(fā)，抗震工程中的結(jié)構(gòu)可以通過結(jié)構(gòu)動力學(xué)理論進(jìn)行分析與計算，但其中需要解決龐雜的微分方程。Matlab是一個強(qiáng)大的矩陣計算數(shù)學(xué)軟件，它是基于C語言開發(fā)，但是其編程又不需要繁瑣的電腦編程語法知識，十分適合工科學(xué)生在進(jìn)行研究時學(xué)習(xí)與參考，通過集中質(zhì)量法對模型進(jìn)行簡化，便可以用Matlab Simulink仿真利用空間狀態(tài)函數(shù)進(jìn)行微分方程求解，通過振型分解的方法求得各個Simulink模塊的系數(shù)進(jìn)行仿真[1]。

1 結(jié)構(gòu)動力學(xué)方程

Matlab是Matrix Laboratory的簡稱，意為是矩陣實(shí)驗室。顧名思義Matlab是以計算矩陣和數(shù)組為核心軟件，對矩陣迭代編程計算的功能相當(dāng)強(qiáng)大。Simulink是Matlab中一個用來對動態(tài)系統(tǒng)進(jìn)行建模、仿真和分析的軟件包，它具有直觀、簡便、易于理解的特點(diǎn)。本文采用振型分解法計算其振型，頻率，為了簡單考慮采用集中質(zhì)量模型。結(jié)構(gòu)的運(yùn)動微分方程為：

式中：M0，k0分別為原結(jié)構(gòu)的質(zhì)量矩陣、剛度矩陣；

在上式中：

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

將上述的方程化簡為:

在此方程左右兩邊同時乘以XT,經(jīng)過化簡可得

x(t)=Xq(t)

可以這樣理解上式體系的位移可以看作是由各個振型分別乘以相應(yīng)組合系數(shù)q1(t)，q2(t)后疊加而成。換句話說，這種方法實(shí)際上是位移按振型加以分解，所以求出廣義坐標(biāo)q的列向量，然后與各階振型相乘即可得出所要求得的x的相應(yīng)加速度位移坐標(biāo)。

2 基于Matlab結(jié)構(gòu)動力特性各編程求解

經(jīng)過上述振型分解法的原理簡介，Matlab式可以求得多自由度體系的各階振型與頻率，然后通過振型分解法的原理利用Matlab編寫程序?qū)Y(jié)構(gòu)進(jìn)行分析求解，Matlab獨(dú)有的eig函數(shù)可以快捷的求解出結(jié)構(gòu)整體剛度矩陣k，質(zhì)量矩陣m的特征值特征向量，通過整理和化簡即是所需求的的振型和頻率。代碼為：

根據(jù)結(jié)構(gòu)計算出層間剛度與集中質(zhì)量編程處理

k=[1,1,1,1,1,1,1]*4.2e7; %層間剛度根據(jù)D值法可計算出各層剛度

kcj=k

mc=[1,1,1,1,1,1,1]*2.0e5；%集中質(zhì)量法質(zhì)量矩陣

m=diag(mc)；

cn=6；

O=zeros(cn)；

I=eye(cn)；

cjgdu

k0=cjgdu(kcj,cn)； %整體剛度矩陣cjgdu為整體剛度矩陣的求法，

[x,d]=eig(kcju,m)；%用Matlab自帶的eig函數(shù)進(jìn)行振型和自振頻率求解

d1=sqrt(d)；

w=diag(sqrt(d))；

for i=1:7；

for j=1:7；

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

end

鑒于人文社會科學(xué)的研究成果難以轉(zhuǎn)化，且存在社會效益優(yōu)于經(jīng)濟(jì)效益、長期效益優(yōu)于短期效益的特點(diǎn)，本文研究的人文社會科學(xué)科研項目的人員費(fèi)用主要包括工資費(fèi)、勞務(wù)費(fèi)、專家咨詢費(fèi)和績效支出。截至目前的有關(guān)內(nèi)容梳理如下：

end

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

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

c0=a1*m+a2*k0

3 Simulink仿真模型

通過上述代碼可以求得各個所需矩陣，然后利用Simulink仿真結(jié)構(gòu)動力學(xué)微分方程進(jìn)行模擬某7層鋼筋混凝土框架。層間均高3 m，層間質(zhì)量為2.0e5 kg，層間剪切剛度根據(jù)D值法計算可得4.2e7 N/m，用集中質(zhì)量法簡化參見圖1。

圖1 模型簡化圖

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

在Matlab工作空間中用下列代碼求出各個模塊的參數(shù)

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函數(shù)

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模型進(jìn)行對比，采取正弦波時程分析，信號源為正弦函數(shù)波，從scope中輸出頂層加速度的反應(yīng)參見圖3，峰值為1.16。

圖3 加速度反應(yīng)值

4 Sap2000模型分析

根據(jù)所示模型進(jìn)行如下所示Sap2000建立的平面結(jié)構(gòu) 模型，參見圖4、圖5。

圖4 平面模型

圖5 加速度反應(yīng)值

其頂層加速度大小與Matlab建立模型大小基本相同，驗證了結(jié)論正確性。

可見Matlab的Simulink模塊在數(shù)值分析編程上的應(yīng)用是可行的。

5 結(jié)論

(1)文中分別用上述兩種算法進(jìn)行求解，然后進(jìn)行比較,Sap2000計算結(jié)果與Matlab編程求解的結(jié)果非常相似，考慮到結(jié)構(gòu)體系的一些假設(shè)條件從而存在誤差,因此,基于Matlab程序Simulink仿真模擬是一個十分可行的方法。

(2)Matlab是對于求解結(jié)構(gòu)的受力特點(diǎn)，動力特性以及求解復(fù)雜的結(jié)構(gòu)動力學(xué)微分方程有強(qiáng)大的計算能力。

(3)結(jié)構(gòu)編程計算Matlab是一個具有友好界面的編程工具，適合科研計算工作者的學(xué)習(xí)。

[1] GB 50011-2001建筑抗震設(shè)計規(guī)范[S]

[2] 徐趙東.Matlab語言在抗震工程中的應(yīng)用[M].北京:科學(xué)出版社,2004

[3] 徐斌.Matlab有限元結(jié)構(gòu)動力學(xué)分析與工程應(yīng)用[M].北京:清華大學(xué)出版社,2009