1、?uUi 1 2Ui Ui 1而離散時間點的運動為? ?Uiu(ti),Uiu(ti), Ui?u(ti)( i? ?由體系運動方程為： 譏山）cu（t）ku（t）(b)(c)將速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到ti時中心差分法計算程序編程姓名：張澤偉 學號：電話:一、中心差分法程序原理說明1.1中心差分法思路中心差分法的基本思路：是將運動方程中的速度向量和加速度向量用位移的某種組合來表示，將微分方程組的求解問題轉化為代數方程組的求解問題，并在時間區間內求得每個微小時間區間的遞推公式，進而求得整個時程的反應。1.2中心差分法原理中心差分法只在相隔 t 一些離散的時間

2、區間內滿足運動方程，其基于有限差分代替位移對時間的求導（即速度和加速度），如果采用等時間步長，tit，則速度與加速度的中心差分近似為:Ui 1 Ui 1(a)2 t刻的運動方程:mU Ac仏亠 kui 0t22 t(d)在（d）式中，假設Ui和Ui 1是已知的，即在ti及ti以前時刻的運動已知，則可以把已知項移到方程的右邊，整理得到:(號無)Ui 1 (k卑出(巴無)Ui 1t 2 ttt 2 t(e)由式（e）就可以根據1及1以前時刻的運動，求得ti 1時刻的運動，如果需 要可以用式（a）和式（b）求得體系的速度和加速度。1.3初始條件轉化假設給定的初始條件為Uo u（0）,? ?U0 U（

3、0），（ g）由式（g）確定u 1。在零時刻速度和加速度的中心差分公式為：U0U1U 1:t(h)2?U12u0U 1U02、t(1)?t2 ?U 1U0t U0U0將式（|）消去U1得：2(j)? ?而零時刻的加速度值U0可以用t = 0時的運動方程mU0 CU0 kU00確定?1 “?U0(CU0kU0)即m(k)?這樣就可以根據初始條件U0,U0和初始荷載Po，就可以根據上式確定U1的1.4中心差分法編程思路 基本數據準備和初始條件計算:?U0(cu0 ku0) m? 2 ?. t -U 1 U0 t U0U02 計算等效剛度和中心差分計算公式中的相關系數:m ct22 t2ma k 2

4、 t2m c孑7t 根據1及1以前時刻的運動，計算ti 1時刻的運動:P aui bui iUi i P kUi 1 Ui 12 t?uUi i 2ui Ui i 下一步計算用i+1代替i，對于線彈性結構體系，重復第3步，對于非線性結 構體系，重復第2步和第3步。1.5中心差分法穩定條件2以上為中心差分法逐步計算公式，其具有 2階精度，即誤差°（ t）；并且為有條件穩定，穩定條件為:二、程序框圖根據中心差分法的原理，可以得出本程序的主要程序思想，以下面框圖的形式展示出來:開始循環體（步數控制）1r勺計算等效剛度和系數 a、b廠1計算位移、速度、加速度并保存1J三、程序清單%m,k,c

5、分別為質量、剛度、阻尼%pO,dt,t分別為外荷載幅值、時間步距、總時間 %uO,vO為初始條件初位移和初速度%u,v,ac分別為位移、速度、加速度反應ek=等效剛度；p=荷載；ep=等效荷載%定義矩陣XO=input('請按格式和順序輸入初始矩陣，如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=zeros(size(t);u(：,2)=u0;v(

6、:,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)-b*u(:,i-1);%得出所需要結果u(:

7、,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,%分別取出其中的參數：%將時間分步，采用等時間步長；%計算t的向量長度，得出步數；%設定存儲u的矩陣；%設定存儲v的矩陣；%設定存儲ac的矩陣；%賦值向量第2項為uO;%賦值向量第2項為vO；%求出初始加速度acO ；%計算初始條件u-1項；%計算等效剛度；%計算方程系數；%給出初始荷載條件；%計算初始等效荷載；%計算位移u仁u(:,3)%從第二項開始進行中心

8、差分法計%給出荷載條件，按照簡諧荷載計算; %計算等效荷載；%計算位移量；%計算速度量；%計算加速度量；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);%繪制位移、速度、加速度時程曲線 %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

9、)速度(m/s)加速度(m/sA2)'),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('速度

10、v的時程曲線');legend('速度v') subplot(3,1,3),plot(t,ac,'r'),grid,xlabel('時間(s)'),ylabel('加速度(m/sA2)'), title('加速度ac的時程曲線');legend('加速度ac')四、輸入數據本程序采用單自由度體系進行計算，主要已知參數信息如下：其質量M=9240kg、剛度K= 1460KN/m、阻尼系數C 6.41kN?s/m，對結構施加動力荷載P 73000sin(0.5 t)N，結構周期T=0.05s,初

11、始位移Uo 0.05m ,初始速度V0 0m/s，假設結構處于線彈性狀態。由中心差分法可知，要使計算結果穩定且不發散，需滿足：時間步長 t Tn 0.159s，本例分別取時間步長為 0.1s、0.15s、0.17s、0.2s分別進行計算，并驗證其穩定條件，取總時間為30s。則:X0=9240 1460000 6410 0.05 0 20 73000 0.05五、計算結果當 dt=0.1s:位移社的時程曲線61015202530吋間抑速度恥的時程曲線當 dt=0.15s 時:時間倒速度甲的時程曲線O時間加速度即的時程曲線加速度L5202530'%時 ®(s)20當 dt=0.17s 時:O氾D£ O 巧£魯位移的吋程曲線|11位移u1ii J2530*15O|»X _iltl O2 0 2-翌SW1520時間倒速度“的時程曲線II遠度¥ 11 1li!幵1535OFX101520時間抑速度機的時程曲線:1加速曲訂1fj 廿廠125101520時間芮)33當 dt=0.2s 時:Illi*a«J電i位移uap11iii1K 100位移口的時程曲線61015202530速度”的時程曲第%f1:-1速度V /:1iii1X 10860520J63015時間 加速度祉的時程曲線%pa加速度丸"7r11