中心差分法計(jì)算程序編程

1、中心差分法計(jì)算程序編程姓名：張澤偉 學(xué)號(hào)：電話:一、中心差分法程序原理說(shuō)明1.1中心差分法思路中心差分法的基本思路：是將運(yùn)動(dòng)方程中的速度向量和加速度向量用位移的 某種組合來(lái)表示，將微分方程組的求解問(wèn)題轉(zhuǎn)化為代數(shù)方程組的求解問(wèn)題，并在時(shí)間區(qū)間內(nèi)求得每個(gè)微小時(shí)間區(qū)間的遞推公式，進(jìn)而求得整個(gè)時(shí)程的反應(yīng)。12中心差分法原理中心差分法只在相隔 氏一些離散的時(shí)間區(qū)間內(nèi)滿足運(yùn)動(dòng)方程，其基于有限 差分代替位移對(duì)時(shí)間的求導(dǎo)（即速度和加速度），如果采用等時(shí)間步長(zhǎng)，：t ：t則速度與加速度的中心差分近似為:u 1 -山(a)2:tUi i -25 Uij(b)而離散時(shí)間點(diǎn)的運(yùn)動(dòng)為* w(i =0,1,23)Ui =

2、u(ti),Ui =u(ti),Ui = u(ti)由體系運(yùn)動(dòng)方程為：mu（ti） cu（t） ku（t） =0（c）將速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到ti時(shí)刻的運(yùn)動(dòng)方程:5 1 -25 Ui4At2Ui 1 - Uy2=t(d)在（d）式中，假設(shè)Ui和Ui4是已知的，即在ti及1以前時(shí)刻的運(yùn)動(dòng)已知，則 可以把已知項(xiàng)移到方程的右邊，整理得到：(號(hào) 丄)Uii(k-卑)Ui-(號(hào)-丄 m(e)=t2 2=t：t2-t2 2=t由式（e）就可以根據(jù)1及1以前時(shí)刻的運(yùn)動(dòng)，求得ti1時(shí)刻的運(yùn)動(dòng)，如果需要可以用式（a）和式（b）求得體系的速度和加速度1.3初始條件轉(zhuǎn)化假設(shè)給

3、定的初始條件為Uo 二 u（0）,Uo = u(0),(g)由式（g）確定u。在零時(shí)刻速度和加速度的中心差分公式為:UoUi -U 12-t(h)將式（i）消去Ul得:UoUi - 2uo u itLt Uo.：t2 1Uo2(j)mUo cuo kUo = o(k)而零時(shí)刻的加速度值Uo可以用t = o時(shí)的運(yùn)動(dòng)方程 確定1uo ( _cuo _ kUo)m這樣就可以根據(jù)初始條件Uo,Uo和初始荷載Po，就可以根據(jù)上式確定的 值。1.4中心差分法編程思路 基本數(shù)據(jù)準(zhǔn)備和初始條件計(jì)算:1uo (-cuo-kuo)mt2 ”U 4 二 Uo - Uouo2 計(jì)算等效剛度和中心差分計(jì)算公式中的相關(guān)系

4、數(shù):kX2ma = k 2c2：t 根據(jù)1及1以前時(shí)刻的運(yùn)動(dòng)，計(jì)算ti 1時(shí)刻的運(yùn)動(dòng):P 二-as -bUiui 1 ui J2：tUi i -2ui At2 下一步計(jì)算用i+1代替i，對(duì)于線彈性結(jié)構(gòu)體系，重復(fù)第3步，對(duì)于非線性結(jié) 構(gòu)體系，重復(fù)第2步和第3步。1.5中心差分法穩(wěn)定條件以上為中心差分法逐步計(jì)算公式，其具有 2階精度，即誤差 “ （鳳2）；并且為M Tn有條件穩(wěn)定，穩(wěn)定條件為：二二、程序框圖根據(jù)中心差分法的原理，可以得出本程序的主要程序思想，以下面框圖的形式展示出來(lái)：編輯回調(diào)函數(shù)I）調(diào)用輸入數(shù)據(jù)計(jì)算初始加速度和位移前一時(shí)刻值開始循環(huán)體（步數(shù)控制）1ri計(jì)算等效剛度和系數(shù) a、b廠1

5、計(jì)算位移、速度、加速度并保存1J三、程序清單%m,k,c分別為質(zhì)量、剛度、阻尼%pO,dt,t分別為外荷載幅值、時(shí)間步距、總時(shí)間 %uO,vO為初始條件初位移和初速度%u,v,ac分別為位移、速度、加速度反應(yīng)ek=等效剛度；p=荷載；ep=等效荷載%定義矩陣XO=input(請(qǐng)按格式和順序輸入初始矩陣，如m=X0(1,1);k=X0(1,2);c=X0(1,3);u0=X0(1,4);v0=X0(1,5);t=X0(1,6);P0=X0(1,7);dt=X0(1,8) t=O:dt:t;mm, nn =size(t);u=zeros(size(t);v=zeros(size(t);ac=zer

6、os(size(t);u(：,2)=u0;v(:,2)=v0;ac(:,2)=(P0-c*v(:,2)-k*u(:,2)/m;u(:,1)=u(:,2)-dt*v(:,2)+(dt)A2)*ac(:,2)/ 2;ek=m/(dtA2)+c/(2*dt);a=k-(2*m)/(dtA2);b=m/(dtA2)-c/(2*dt);p(:,2)=P0*si n(0);ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1); u(:,3)=ep(:,2)/ek;for i=3:nn算；p(:,i)=P0*si n(.5*pi*(i-2)*dt);ep(:,i)=p(:,i)-a*u(:,i)

7、-b*u(:,i-1);%得出所需要結(jié)果u(:,i+1)=ep(:,i)/ek;v(:,i)=(u(:,i+1)-u(:,i-1)/(2*dt);ac(:,i)=(u(:,i+1)-2*u(:,i)+u(:,i-1)/(dtA2);endXO=m,k,c,uO,vO,t,PO,dt,%分別取出其中的參數(shù)：%將時(shí)間分步，采用等時(shí)間步長(zhǎng)；%計(jì)算t的向量長(zhǎng)度，得出步數(shù)；%設(shè)定存儲(chǔ)u的矩陣；%設(shè)定存儲(chǔ)v的矩陣；%設(shè)定存儲(chǔ)ac的矩陣；%賦值向量第2項(xiàng)為uO;%賦值向量第2項(xiàng)為vO；%求出初始加速度acO ；%計(jì)算初始條件u-1項(xiàng)；%計(jì)算等效剛度；%計(jì)算方程系數(shù)；%給出初始荷載條件；%計(jì)算初始等效荷載；%

8、計(jì)算位移u仁u(:,3)%從第二項(xiàng)開始進(jìn)行中心差分法計(jì)%給出荷載條件，按照簡(jiǎn)諧荷載計(jì)算; %計(jì)算等效荷載；%計(jì)算位移量；%計(jì)算速度量；%計(jì)算加速度量；t=t(:,1:e nd-1);u=u(:,2:e nd-1);v=v(:,2:e nd);ac=ac(:,2:e nd);p=p(:,2:e nd);ep=ep(:,2:e nd);%繪制位移、速度、加速度時(shí)程曲線 %plot(t,u,b-o),hold on,plot(t,v,g-p),hold on,plot(t,ac,r:x),grid on,xlabel(時(shí)間(s),ylabel( 位移(m)速度(m/s)加速度(m/sA2),titl

9、e(頂層u,v,ac的時(shí)程曲線);subplot(3,1,1),plot(t,u,b-),grid,xlabel(時(shí)間(s),ylabel(位移(m),title(位移u的時(shí)程曲線);legend(位移u)subplot(3,1,2),plot(t,v,k),grid,xlabel(時(shí)間(s),ylabel(速度(m/s),title(速度v的時(shí)程曲線);legend(速度v)subplot(3,1,3),plot(t,ac,r),grid,xlabel(時(shí)間(s),ylabel(加速度(m/sA2),title(加速度ac的時(shí)程曲線);legend(加速度ac)四、輸入數(shù)據(jù)本程序采用單自由度

10、體系進(jìn)行計(jì)算，主要已知參數(shù)信息如下：其質(zhì)量M=9240kg、剛度K= 1460KN/m、阻尼系數(shù)C = 6.41kN s/m，對(duì)結(jié)構(gòu) 施加動(dòng)力荷載P =73000sin(0.5二t)N，結(jié)構(gòu)周期T=0.05s,初始位移u。= 0.05m , 初始速度 =0m/s，假設(shè)結(jié)構(gòu)處于線彈性狀態(tài)。由中心差分法可知，要使計(jì)算結(jié)果穩(wěn)定且不發(fā)散，需滿足：時(shí)間步長(zhǎng)Tn = 0.159s，本例分別取時(shí)間步長(zhǎng)為 0.1s、0.15s、n0.17s、0.2s分別進(jìn)行計(jì)算，并驗(yàn)證其穩(wěn)定條件，取總時(shí)間為30s。則：X0=9240 1460000 6410 0.05 0 20 73000 0.05五、計(jì)算結(jié)果當(dāng) dt=0.

11、1s:時(shí)間速度的時(shí)程曲線O 4時(shí)間力口速度恥的時(shí)程曲線1030 苦2015冋時(shí)O 15旦 STSS當(dāng) dt=0.15s 時(shí):時(shí)間速度的時(shí)程曲線01 0趣121煩時(shí)間 加速度兀的時(shí)程曲線500510152011I1”-扣速度孔hillfwifVy叭p即山: * F丄126H604320當(dāng) dt=0.17s 時(shí):位移u的時(shí)程曲線T位移U112530101520時(shí)間速度“的時(shí)程曲線0L1速度訂j1j12皿525300101520時(shí)間加速度韭的時(shí)程曲線i加速度禮11; Tix 105361020253016時(shí)間當(dāng) dt=0.2s 時(shí):0-50I1i位移ui: L5 *涉位移u的時(shí)程曲線5101520時(shí)間*) 速度的時(shí)程曲線253011!D1速度粧1ii常k一- 卜1x 10Se05101520時(shí)間 加速度機(jī)的時(shí)程曲線2530加速度ac I 1105x 1評(píng)101520時(shí)間2530六、結(jié)果穩(wěn)定性分析由以上時(shí)程圖可以得到當(dāng) 氏=0.1, 0