中心差分法計算程序編程.doc

中心差分法計算程序編程姓名：張澤偉 學號： 電話：一、中心差分法程序原理說明 1.1 中心差分法思路中心差分法的基本思路：是將運動方程中的速度向量和加速度向量用位移的某種組合來表示，將微分方程組的求解問題轉化為代數(shù)方程組的求解問題，并在時間區(qū)間內求得每個微小時間區(qū)間的遞推公式，進而求得整個時程的反應。1.2 中心差分法原理中心差分法只在相隔一些離散的時間區(qū)間內滿足運動方程，其基于有限差分代替位移對時間的求導（即速度和加速度），如果采用等時間步長，則速度與加速度的中心差分近似為： (a) (b)而離散時間點的運動為 ( 0,1,2,3,)由體系運動方程為： (c)將速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到時刻的運動方程： （d）在（d）式中，假設和是已知的，即在及以前時刻的運動已知，則可以把已知項移到方程的右邊，整理得到： (e)由式（e）就可以根據(jù)及以前時刻的運動，求得時刻的運動，如果需要可以用式（a）和式（b）求得體系的速度和加速度。1.3 初始條件轉化假設給定的初始條件為 （g）由式（g）確定。在零時刻速度和加速度的中心差分公式為： （h） （i）將式（i）消去得： （j） 而零時刻的加速度值可以用t0時的運動方程 確定即 （k）這樣就可以根據(jù)初始條件和初始荷載，就可以根據(jù)上式確定的值。1.4 中心差分法編程思路 基本數(shù)據(jù)準備和初始條件計算： 計算等效剛度和中心差分計算公式中的相關系數(shù)： 根據(jù)及以前時刻的運動，計算時刻的運動： 下一步計算用i+1代替i，對于線彈性結構體系，重復第3步，對于非線性結構體系，重復第2步和第3步。1.5 中心差分法穩(wěn)定條件以上為中心差分法逐步計算公式，其具有2階精度，即誤差；并且為有條件穩(wěn)定，穩(wěn)定條件為： 2、 程序框圖根據(jù)中心差分法的原理，可以得出本程序的主要程序思想，以下面框圖的形式展示出來：編輯回調函數(shù)調用輸入數(shù)據(jù)計算初始加速度和位移前一時刻值開始循環(huán)體（步數(shù)控制）計算等效剛度和系數(shù)a、b計算位移、速度、加速度并保存分析穩(wěn)定條件繪圖3、 程序清單%m,k,c分別為質量、剛度、阻尼%p0,dt,t分別為外荷載幅值、時間步距、總時間%u0,v0為初始條件初位移和初速度%u,v,ac分別為位移、速度、加速度反應ek=等效剛度；p=荷載；ep=等效荷載%定義矩陣X0=input(請按格式和順序輸入初始矩陣，如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; %將時間分步，采用等時間步長；mm,nn=size(t); %計算t的向量長度，得出步數(shù)；u=zeros(size(t); %設定存儲u的矩陣；v=zeros(size(t); %設定存儲v的矩陣；ac=zeros(size(t); %設定存儲ac的矩陣；u(:,2)=u0; %賦值向量第2項為u0；v(:,2)=v0; %賦值向量第2項為v0；ac(:,2)=(P0-c*v(:,2)-k*u(:,2)/m; %求出初始加速度ac0；u(:,1)=u(:,2)-dt*v(:,2)+(dt)2)*ac(:,2)/2; %計算初始條件u-1項；ek=m/(dt2)+c/(2*dt); %計算等效剛度；a=k-(2*m)/(dt2);b=m/(dt2)-c/(2*dt); %計算方程系數(shù)；p(:,2)=P0*sin(0); %給出初始荷載條件；ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1); %計算初始等效荷載；u(:,3)=ep(:,2)/ek; %計算位移u1=u(:,3)for i=3:nn %從第二項開始進行中心差分法計算；p(:,i)=P0*sin(.5*pi*(i-2)*dt); %給出荷載條件，按照簡諧荷載計算；ep(:,i)=p(:,i)-a*u(:,i)-b*u(:,i-1); %計算等效荷載；%-得出所需要結果-%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)/(dt2); %計算加速度量；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);%-繪制位移、速度、加速度時程曲線-%plot(t,u,b-o),hold on,plot(t,v,g-p),hold on,plot(t,ac,r:x),grid on,xlabel(時間(s),ylabel(位移(m)速度(m/s)加速度(m/s2),title(頂層u,v,ac的時程曲線);subplot(3,1,1),plot(t,u,b-),grid,xlabel(時間(s),ylabel(位移(m),title(位移u的時程曲線);legend(位移u)subplot(3,1,2),plot(t,v,k),grid,xlabel(時間(s),ylabel(速度(m/s),title(速度v的時程曲線);legend(速度v)subplot(3,1,3),plot(t,ac,r),grid,xlabel(時間(s),ylabel(加速度(m/s2),title(加速度ac的時程曲線);legend(加速度ac)4、 輸入數(shù)據(jù)本程序采用單自由度體系進行計算，主要已知參數(shù)信息如下：其質量M=9240kg、剛度K1460KN/m、阻尼系數(shù)，對結構施加動力荷載，結構周期T=0.05s，初始位移，初始速度，假設結構處于線彈性狀態(tài)。由中心差分法可知，要使計算結果穩(wěn)定且不發(fā)散，需滿足：時間步長，本例分別取時間步長為、分別進行計算，并驗證其穩(wěn)定條件，取總時間為30s。則：X0=9240 1460000 6410 0.05 0 20 73000 0.055、 計算結果當dt=0.1s:當dt=0.15s時：當dt=0.17s時：當dt=0.2s時：6、 結果穩(wěn)定性分析由以上時程圖可以得到當0.1，0.15時逐步計算結果給出的結構運動趨向收斂的，即計算結果是穩(wěn)定