結(jié)構(gòu)動(dòng)力學(xué)電算報(bào)告(附源代碼)

結(jié)構(gòu)動(dòng)力學(xué)學(xué)電算報(bào)告摘要：本次電算以Matlab軟件為平臺(tái)，采用線加速度差值方法對(duì)單自由度體系的彈塑性材料（兩段式恢復(fù)力模型）編寫程序，了解彈塑性材料單自由度體系的運(yùn)動(dòng)過(guò)程與特點(diǎn)及數(shù)值仿真在學(xué)習(xí)中的重要性。關(guān)鍵詞：彈塑性，線加速度，振動(dòng)1.正文軟件簡(jiǎn)介MATLAB是由美國(guó)MATHWORKS公司發(fā)布的主要面對(duì)科學(xué)計(jì)算、可視化以及交互式程序設(shè)計(jì)的高科技計(jì)算環(huán)境。它將數(shù)值分析、矩陣計(jì)算、科學(xué)數(shù)據(jù)可視化以及非線性動(dòng)態(tài)系統(tǒng)的建模和仿真等諸多強(qiáng)大功能集成在一個(gè)易于使用的視窗環(huán)境中，為科學(xué)研究、工程設(shè)計(jì)以及必須進(jìn)行有效數(shù)值計(jì)算的眾多科學(xué)領(lǐng)域提供了一種全面的解決方案，并在很大程度上擺脫了傳統(tǒng)非交互式程序設(shè)計(jì)語(yǔ)言（如C、Fortran）的編輯模式，代表了當(dāng)今國(guó)際科學(xué)計(jì)算軟件的先進(jìn)水平。題型簡(jiǎn)述1.3.計(jì)算公式推導(dǎo)設(shè)加速度方程：ay(t)y(t+t)=y(t)+ tat貝1」：y(t+t)=y(t)+y(t)t+—ay(t)t22at11y(t+t)=y(t)+y(t)t+y(t)t2+ ay(t)t32 6at所以：ay(t)=y(t)at+ ay(t)at22at11ay(t)=y(t)at+y(t)at2+ ay(t)at32 6at故:ay-(t)= 6ay(t)-—y(t)-3y(t)(at)2 at3 atay(t)= ay(t)-3y(t)- y(t)at 2又：may(t)+c(t)ay(t)+k(t)ay(t)=aP(t)my( t)+ c( t) y(+t) k( t) y=( t) Pt且有：ay=AP/k6 at可得：AP(t) =AP(t)+ m[ y(t) +3y(t)] + c(t)[3y(t) + y(t)]At 263k(t)= m+c(t)+k(t)(At)2 Aty(t+At)=y(t)+Ay(t)并且：y(t+At)=y(t)+Ay(t)y(t+At)=y(t)+Ay(t)若已知：y，y，P，P，k，At(0)(0)(0)(1)(0)可算出：y，Ay(0)，Ay(0)，Ay(0)(0)進(jìn)而可算出：y，y，y(At) (At) (At)如此循環(huán)即可求出質(zhì)點(diǎn)各時(shí)程的位移1.4.程序流程圖

口

算出速度反向時(shí)的位移和下次進(jìn)入塑性階段的位移判斷是否進(jìn)入塑形階段？計(jì)算 y'(t+At) ，k(t+At)，AP(t+At)，Ay(t+At)，Ay(t+At)判斷是否繼續(xù)運(yùn)算?繪制位移時(shí)程曲和恢復(fù)力時(shí)程曲線結(jié)束1.5.程序說(shuō)明該程序程序運(yùn)行后，可根據(jù)提示先后輸入：質(zhì)量m，彈性階段剛度k，進(jìn)入塑性階段的位移Y,時(shí)段長(zhǎng)At，阻尼系數(shù)c，外荷載P(t)，初始位移y(0)，初始速度y(0)。其中須注意的是彈性階段剛度k為常數(shù)，阻尼系數(shù)c也為常數(shù)。外荷載P(t)采用一維

矩陣散點(diǎn)輸入也可以輸入荷載P(t)與時(shí)段長(zhǎng)t之間的函數(shù)關(guān)系。另外，如果塑性階段的位移y很大，其振動(dòng)最大位移y 不可能達(dá)到Y(jié)時(shí)，本程序計(jì)P max算的是線彈性體系。1.6.源程序functionconstruction()clear;m=input('請(qǐng)輸入質(zhì)量m');k=input('請(qǐng)輸入彈性階段剛度k');Y=input('請(qǐng)輸入進(jìn)入塑性階段的位移y');t=input('請(qǐng)輸入時(shí)段長(zhǎng)At')');c=input('請(qǐng)輸入阻尼系數(shù)c');p=input('請(qǐng)輸入外荷載P');T=input('請(qǐng)輸入計(jì)算時(shí)間長(zhǎng)度');q=T/t;n=length(p);xmax3=-Y;xmax1=Y;fori=n+1:q+1p(i)=0;endy(1)=input('請(qǐng)輸入初始位移y(0)');y1(1)=input('請(qǐng)輸入初始速度v(0)');y2(1)=(p(1)-c*y1(1)-k*y(1))/m;pp(1)=p(2)-p(1)+(6*y1(1)/t+3*y2(1))*m+c*(3*y1(1)+t*y2(1)*0.5);kk(1)=k+6*m/(t*t)+3*c/t;yy(1)=pp(1)/kk(1);yy1(1)=3*yy(1)/t-3*y1(1)-t*0.5*y2(1);z=1;i=2;f(1)=k*y(1);fori=2:qy(i)=y(i-1)+yy(i-1);y1(i)=y1(i-1)+yy1(i-1);ifz==1y2(i)=(p(i)-c*y1(i)-k*(y(i)-xmax3/2-xmax1/2))/m;kk(i)=k+6*m/(t*t)+3*c/t;pp(i)=p(i+1)-p(i)+(6*y1(i)/t+3*y2(i))*m+c*(3*y1(i)+t*y2(i)*0.5);yy(i)=pp(i)/kk(i);yy1(i)=3*yy(i)/t-3*y1(i)-t*0.5*y2(i);f(i)=k*y(i);if(y(i)+yy(i))>(xmax3+2*Y)z=2;xmax1=y(i);continueendendifz==2y2(i)=(p(i)-c*y1(i)-k*Y)/m;kk(i)=6*m/(t*t)+3*c/t;pp(i)=p(i+1)-p(i)+(6*y1(i)/t+3*y2(i))*m+c*(3*y1(i)+t*y2(i)*0.5);yy(i)=pp(i)/kk(i);yy1(i)=3*yy(i)/t-3*y1(i)-t*0.5*y2(i);f(i)=k*Y;if(y1(i)+yy1(i))*y1(i)<0z=3;xmax2=y(i);continueendendifz==3y2(i)=(p(i)-c*y1(i)-k*(y(i)-(xmax2-xmax1)))/m;kk(i)=k+6*m/(t*t)+3*c/t;pp(i)=p(i+1)-p(i)+(6*y1(i)/t+3*y2(i))*m+c*(3*y1(i)+t*y2(i)/2);yy(i)=pp(i)/kk(i);yy1(i)=3*yy(i)/t-3*y1(i)-t*0.5*y2(i);f(i)=k*y(i);if(y(i)+yy(i))<xmax2-2*xmax1z=4;continueendendifz==4y2(i)=(p(i)-c*y1(i)+k*Y)/m;kk(i)=6*m/(t*t)+3*c/t;pp(i)=p(i+1)-p(i)+(6*y1(i)/t+3*y2(i))*m+c*(3*y1(i)+t*y2(i)/2);yy(i)=pp(i)/kk(i);yy1(i)=3*yy(i)/t-3*y1(i)-t*0.5*y2(i);f(i)=k*Y;if(y1(i)+yy1(i))*y1(i)<0z=1;xmax3=y(i);continueendendends=size(y);s=s(2);x=0:t:(s-1)*t;subplot(1,2,1);plot(x,y);holdon;subplot(1,2,2);plot(x,f,'r');2.例題.2.1用上述程序計(jì)算克拉夫《結(jié)構(gòu)動(dòng)力學(xué)（第二版）》P103例題E7-2，并與克拉夫手算結(jié)果進(jìn)行比較。2.2程序運(yùn)行輸入：>>shiyan3i青輸入質(zhì)里m0.1請(qǐng)輸入.彈性階段剛度k5i青輸入進(jìn)入塑性階段的位移y1.2請(qǐng)輸入步長(zhǎng)注]0.1請(qǐng)輸入阻尼系數(shù)匚0.2請(qǐng)輸人外荷載P[058753210]i青輸入計(jì)算時(shí)間長(zhǎng)度T5請(qǐng)輸入初始位移0請(qǐng)輸入初始速度巧⑴）02.3計(jì)算結(jié)果進(jìn)行比較（藍(lán)色為電算結(jié)果，紅色為克拉夫手算結(jié)果）位移時(shí)程曲線（y-t）

恢復(fù)力時(shí)程曲線(fs-t)s2.4結(jié)果分析：當(dāng)外荷載撤消后，體系在阻尼作用下，振動(dòng)逐漸削減，最后靜止不動(dòng)。但由于振動(dòng)過(guò)程中，體系會(huì)進(jìn)入塑性階段產(chǎn)生不可回復(fù)的塑性位移，所以體系靜止位置與其平衡位置不一致，這些在位移時(shí)程圖上均可很好地反映出來(lái)。同時(shí)，克拉夫在計(jì)算P(t)存在一定的問題，因