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

1、中心差分法計(jì)算程序編程姓名：張澤偉 學(xué)號(hào)： 電話(huà)：一、中心差分法程序原理說(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)。1.2 中心差分法原理中心差分法只在相隔一些離散的時(shí)間區(qū)間內(nèi)滿(mǎn)足運(yùn)動(dòng)方程，其基于有限差分代替位移對(duì)時(shí)間的求導(dǎo)（即速度和加速度），如果采用等時(shí)間步長(zhǎng)，則速度與加速度的中心差分近似為： (a) (b)而離散時(shí)間點(diǎn)的運(yùn)動(dòng)為 ( 0,1,2,3,)由體系運(yùn)動(dòng)方程為： (c)將速度和加速度的差分近似公式（a

2、）和式（b）代入式（c）可以得到時(shí)刻的運(yùn)動(dòng)方程： （d）在（d）式中，假設(shè)和是已知的，即在及以前時(shí)刻的運(yùn)動(dòng)已知，則可以把已知項(xiàng)移到方程的右邊，整理得到： (e)由式（e）就可以根據(jù)及以前時(shí)刻的運(yùn)動(dòng)，求得時(shí)刻的運(yùn)動(dòng)，如果需要可以用式（a）和式（b）求得體系的速度和加速度。1.3 初始條件轉(zhuǎn)化假設(shè)給定的初始條件為 （g）由式（g）確定。在零時(shí)刻速度和加速度的中心差分公式為： （h） （i）將式（i）消去得： （j） 而零時(shí)刻的加速度值可以用t0時(shí)的運(yùn)動(dòng)方程 確定即 （k）這樣就可以根據(jù)初始條件和初始荷載，就可以根據(jù)上式確定的值。1.4 中心差分法編程思路 基本數(shù)據(jù)準(zhǔn)備和初始條件計(jì)算： 計(jì)算等效剛度

3、和中心差分計(jì)算公式中的相關(guān)系數(shù)： 根據(jù)及以前時(shí)刻的運(yùn)動(dòng)，計(jì)算時(shí)刻的運(yùn)動(dòng)： 下一步計(jì)算用i+1代替i，對(duì)于線彈性結(jié)構(gòu)體系，重復(fù)第3步，對(duì)于非線性結(jié)構(gòu)體系，重復(fù)第2步和第3步。1.5 中心差分法穩(wěn)定條件以上為中心差分法逐步計(jì)算公式，其具有2階精度，即誤差；并且為有條件穩(wěn)定，穩(wěn)定條件為： 2、 程序框圖根據(jù)中心差分法的原理，可以得出本程序的主要程序思想，以下面框圖的形式展示出來(lái)：編輯回調(diào)函數(shù)調(diào)用輸入數(shù)據(jù)計(jì)算初始加速度和位移前一時(shí)刻值開(kāi)始循環(huán)體（步數(shù)控制）計(jì)算等效剛度和系數(shù)a、b計(jì)算位移、速度、加速度并保存分析穩(wěn)定條件繪圖3、 程序清單%m,k,c分別為質(zhì)量、剛度、阻尼%p0,dt,t分別為外荷載幅值

4、、時(shí)間步距、總時(shí)間%u0,v0為初始條件初位移和初速度%u,v,ac分別為位移、速度、加速度反應(yīng)ek=等效剛度；p=荷載；ep=等效荷載%定義矩陣X0=input('請(qǐng)按格式和順序輸入初始矩陣，如X0=m,k,c,u0,v0,t,P0,dt, m=X0(1,1);k=X0(1,2);c=X0(1,3);u0=X0(1,4); %分別取出其中的參數(shù)：v0=X0(1,5);t=X0(1,6);P0=X0(1,7);dt=X0(1,8)t=0:dt:t; %將時(shí)間分步，采用等時(shí)間步長(zhǎng)；mm,nn=size(t); %計(jì)算t的向量長(zhǎng)度，得出步數(shù)；u=zeros(size(t); %設(shè)定存儲(chǔ)u的

5、矩陣；v=zeros(size(t); %設(shè)定存儲(chǔ)v的矩陣；ac=zeros(size(t); %設(shè)定存儲(chǔ)ac的矩陣；u(:,2)=u0; %賦值向量第2項(xiàng)為u0；v(:,2)=v0; %賦值向量第2項(xiàng)為v0；ac(:,2)=(P0-c*v(:,2)-k*u(:,2)/m; %求出初始加速度ac0；u(:,1)=u(:,2)-dt*v(:,2)+(dt)2)*ac(:,2)/2; %計(jì)算初始條件u-1項(xiàng)；ek=m/(dt2)+c/(2*dt); %計(jì)算等效剛度；a=k-(2*m)/(dt2);b=m/(dt2)-c/(2*dt); %計(jì)算方程系數(shù)；p(:,2)=P0*sin(0); %給出初始

6、荷載條件；ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1); %計(jì)算初始等效荷載；u(:,3)=ep(:,2)/ek; %計(jì)算位移u1=u(:,3)for i=3:nn %從第二項(xiàng)開(kāi)始進(jìn)行中心差分法計(jì)算；p(:,i)=P0*sin(.5*pi*(i-2)*dt); %給出荷載條件，按照簡(jiǎn)諧荷載計(jì)算；ep(:,i)=p(:,i)-a*u(:,i)-b*u(:,i-1); %計(jì)算等效荷載；%-得出所需要結(jié)果-%u(:,i+1)=ep(:,i)/ek; %計(jì)算位移量；v(:,i)=(u(:,i+1)-u(:,i-1)/(2*dt); %計(jì)算速度量；ac(:,i)=(u(:,i+1)-

7、2*u(:,i)+u(:,i-1)/(dt2); %計(jì)算加速度量；endt=t(:,1:end-1);u=u(:,2:end-1);v=v(:,2:end);ac=ac(:,2:end);p=p(:,2:end);ep=ep(:,2:end);%-繪制位移、速度、加速度時(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/s2)'

8、),title('頂層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('

9、;速度v')subplot(3,1,3),plot(t,ac,'r'),grid,xlabel('時(shí)間(s)'),ylabel('加速度(m/s2)'),title('加速度ac的時(shí)程曲線');legend('加速度ac')4、 輸入數(shù)據(jù)本程序采用單自由度體系進(jìn)行計(jì)算，主要已知參數(shù)信息如下：其質(zhì)量M=9240kg、剛度K1460KN/m、阻尼系數(shù)，對(duì)結(jié)構(gòu)施加動(dòng)力荷載，結(jié)構(gòu)周期T=0.05s，初始位移，初始速度，假設(shè)結(jié)構(gòu)處于線彈性狀態(tài)。由中心差分法可知，要使計(jì)算結(jié)果穩(wěn)定且不發(fā)散，需滿(mǎn)足：時(shí)間步長(zhǎng)，本例分別取時(shí)間步長(zhǎng)為、分別進(jìn)行計(jì)算，并驗(yàn)證其穩(wěn)定條件，取總時(shí)間為30s。則：X0=9240 1460000 6410 0.05 0 20 73000 0.055、 計(jì)算結(jié)果當(dāng)dt=0.1s:當(dāng)dt=0.15s時(shí)：當(dāng)dt=0.17s時(shí)：當(dāng)dt=0.2s時(shí)：6、 結(jié)果穩(wěn)定性分析由以上時(shí)程圖可以得到當(dāng)0.1，0.15時(shí)逐步計(jì)算結(jié)果給出的結(jié)構(gòu)運(yùn)動(dòng)趨向收