1、精選優質文檔-傾情為你奉上中心差分法計算單自由度體系動力反映的報告前言基于疊加原理的時域積分法與頻域積分法一樣，都假設結構在在全部反應過程中都是線性的。而時域逐步積分法只是假設結構本構關系在一個微小的時間步距內是線性的，相當于分段直線來逼近實際的曲線。時域逐步積分法是結構動力問題中研究并應用廣泛的課題。中心差分法是一種目前發展的一系列結構動力反應分析的時域逐步積分法的一種，時域逐步積分法還包括分段解析法、平均常加速度法、線性加速度法、Newmarket-和Wilson-法等。中心差分法(central difference method)原理 中心差分法的基本思路將運動方程中的速度向量和加速度

2、向量用位移的某種組合來表示，將微分方程組的求解問題轉化為代數方程組的求解問題，并在時間區間內求得每個微小時間區間的遞推公式，進而求得整個時程的反應。中心差分法是一種顯示的積分法，它基于用有限差分代替位移對時間的求導（即速度和加速度）。如果采用等時間步長，ti=t（t為常數），則速度與加速度的中心差分近似為ui=ui+1+ui-12t (1)ui=ui+1-2ui+ui-1t2 (2)用u表示位移，離散時間點的運動為：ui=uti，ui=uti，ui=uti (i=0,1,2)體系的運動方程為mut+cut+kut=P(t) (3)將速度和加速度的差分近似公式(1)和(2)代入(3)中得出在ti

3、時刻的運動方程，將方程整理得到ui+1由ui和ui-1表示的兩步法的運動方程(4)：mt2+c2tui+1=Pi-k-2mt2ui-mt2-c2tui-1 (4)這樣就可以根據ti及以前的時刻的運動求得ti+1時刻的運動。中心差分法屬于兩步法，用兩步法計算時存在起步問題，必須要給出相鄰的兩個時刻的位移值，才能逐步計算。對于地震作用下結構的反應問題和一般的零初始條件下的動力問題，可以用(4)直接計算，因為總可以假設初始的兩個時間點（一般取i=0,-1）的位移等于零。但是對應于非零初始條件或零時刻外荷載很大時，需要進行一定的分析，建立兩個起步時刻（即i=0,-1）的位移值。假設給定的初始條件為u0

4、=u0u0=u0 (5)根據初始條件來確定u-1。根據中心差分公式u0=u1+u-12tu0=u1-2u0+u-1t2 (6)消去u1得到u-1的公式：u-1=u0-tu0+t22u0 (7)其中零時刻加速度值u0可以由t=0時的運動方程得到即u0=1mP0-cu0-ku0 (8)這樣就可以根據初始條件得到u-1，然后再將初始條件應用于公式(4)中，逐步求出不同時刻的運動。中心差分法分析時的具體計算步驟：(1) 基本數據準備與初始條件計算已知：初始位移u0、u0和初始荷載值P0來計算u0和u-1u0=1mP0-cu0-ku0u-1=u0-tu0+t22u0 (2) 計算等效剛度和中心差分法計算

5、公式中的系數k=mt2+c2ta=k-2mt2b=mt2-c2t因此中心差分法計算公式可以表示為：kui+1=Pi-aui-bui-1(3) 根據ti及以前的時刻的運動求得ti+1時刻的運動Pi=Pi-aui-bui-1ui+1=Pik(4)下一步計算中用i+1代替i，對于線彈性體系重復第3步計算步驟，對于非線性彈性體系，重復第2和第3計算步驟。以上的中心差分法逐步計算公式具有2階精度，即誤差O(t2)；并且是有條件穩定的，穩定條件為：tTn式中，Tn為結構的自振周期，對于多自由度體系則為結構的最小自振周期。算例本算例根據結構動力學48頁算例3.1數據編寫，穩定條件為dt<=0.16s對

6、于一個單層框架結構，假設樓板剛度無限大，且結構質量集中于樓層，其質量M=9240kg、剛度K1460KN/m、阻尼系數C6.41KNs/m，對結構施加動力荷載P=73000sin0.5t假設結構處于線彈性狀態，用中心差分法計算結構的自由振動反應。采用MATLAB語言編程，并以單自由度體系為例進行計算,設初位移u00.05m和初速度v0=0，取不同的步長分別計算，以驗證中心差分法的穩定條件。先計算，由穩定條件，而rad/s，則，所以本次計算取0.2，0.1，0.05分別進行計算。計算結果與分析 1）當dt=0.05s時，可以得到位移u，速度v，加速度ac的時程曲線如下：2）當dt=01s時，可以

7、得到位移u，速度v，加速度ac的時程曲線如下：3）當dt=0.2s時，可以得到如下提示：不滿足穩定條件:dt<=Tn/pi,請重新輸入符合穩定條件的時間步長dt。附錄%m=質量；k=剛度；c=阻尼；u0=初始位移；v0=初始速度；all_time=所用時間；P0=荷載幅值；dt=時間步長；%u=位移；v=速度；ac=加速度；ek=等效剛度；p=荷載；ep=等效荷載；t=時間；clear %A0=input('請按格式和順序輸入初始矩陣，如A0=m,k,c,u0,v0,all_time,P0,dt=9240 6410 0.05 0 30 7300*3 0.05='); A0

8、=9240 6410 0.05 0 20 73000 0.05;m=A0(1,1);k=A0(1,2);c=A0(1,3);u0=A0(1,4);v0=A0(1,5);all_time=A0(1,6);P0=A0(1,7);dt=A0(1,8); if dt>2*sqrt(m/k) %判斷時間步長是否滿足穩定條件 disp('不滿足穩定條件:dt<=Tn/pi,請重新輸入符合穩定條件的時間步長dt') return elseif 0<dt<=2*sqrt(m/k) disp('滿足穩定條件為：dt<=Tn/pi') endt=0:d

9、t:all_time; %將時間分步，采用等時間步長；mm,nn=size(t); %計算t的向量長度，得出步數；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

10、*dt); %計算等效剛度；a=k-(2*m)/(dt2);b=m/(dt2)-c/(2*dt); %計算方程系數；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(

11、:,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

12、,'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,x

13、label('時間(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')中心差分法求解單自由度體系的自由振動問題算例對于一個單層框架結構，假設樓板剛度無限大，且結構質量集中于樓層，其質量M

14、=2000kg、剛度K50KN/m、阻尼系數C3KNs/m，假設結構處于線彈性狀態，用中心差分法計算結構的自由振動反應。采用MATLAB語言編程，并以單自由度體系為例進行計算,設初位移u00和初速度v0=0，取不同的步長分別計算，以驗證中心差分法的穩定條件。先計算，由穩定條件，而rad/s,則所以本次計算取0.1, 0.3, 0.4, 0.41, 0.42, 0.45分別進行計算MATLAB程序清單function u,v,ac=centraldifferent(M,C,K,u0,v0,time,dt) % 本程序采用中心差分法計算結構的動力響應% 本程序是既可以計算單自由度體系又可以計算多自

15、由度體系，且均假設結構體系處于線彈性狀態；% -%輸入參數%-% M-質量矩陣% C-阻尼矩陣% K-剛度矩陣% u0-初始位移% v0-初始速度% time-模擬時間% dt-時間步長% -%輸出值%-% u-位移% v-速度% ac-加速度% -%中心差分法主要公式及原理%-% MX''+CX'+KX=0 % M*(X(t+dt)-2*X(t)+X(t-dt)/(dt2)+C*(X(t+dt)-X(t-dt)/(2*dt)+K*X(t)=0% (M/dt2+C/2*dt)*(X(t+dt)=-(K-2*M/dt2)*X(t+dt)-(M/dt2-C/2*dt)*X(

16、t-dt)%- 等效剛度Ke等效荷載Pe和相關系數a,b-% Ke=M/dt2+C/2*dt% a=K-2*M/dt2% b=M/dt2-C/2*dt% Pe=-a*X(t)-b*X(t-dt)% X(t+dt)=Pe/Ke% X(t)'=(X(t+dt)-X(t-dt)/(2*dt)% X(t)''=(X(t+dt)-2*X(t)+X(t-dt)/(dt2)% -初始條件-% X0''=(-C*X0'-K*X0)% X(-1)=X0-X0'*dt+X0''*(dt2)/2% -Copyright by zhouhuapi

17、ng(S)-clear allM=input('輸入質量矩陣M ：');C=input('輸入阻尼矩陣C：');K=input('輸入剛度矩陣K：');u0=input('輸入初始位移 u0： ');v0=input('輸入初始速度 v0： ');time=input('輸入模擬時間 time：');dt=input('輸入時間步長dt ：');m,m=size(K);n=time/dt; %計算步數u=zeros(m,floor(n)+1); %設定存儲位移矩陣 v=zeros(m

18、,floor(n)+1); %設定存儲速度矩陣 ac=zeros(m,floor(n)+1); %設定存儲加速度矩陣 P=zeros(m,floor(n)+1); %設定存儲荷載矩陣 u(:,2)=u0; %給定初位移v(:,2)=v0; %給定初速度Ke=M/(dt2)+(C)/(2*dt); %等效剛度Ke及系數a、ba=K-2*M/dt2;b=M/dt2-C/(2*dt);for i=3:1:floor(n)+1; t=(i-2)*dt; ac(:,2)=M(-K*u(:,2)-C*v(:,2); %計算初加速度 u(:,1)=u(:,2)-v(:,2)*dt+(ac(:,2)*(dt2

19、)/2; %計算(0-dt)時刻位移 Pe= -a*u(:,i-1)-b*u(:,i-2); %計算等效荷載Pe u(:,i)=KePe; %計算位移 v(:,i)=(u(:,i)-u(:,i-2)/(2*dt); %計算速度 ac(:,i)=(u(:,i)-2*u(:,i-1)+u(:,i-2)/(dt2); %計算加速度end%-%繪制位移、速度、加速度時程曲線%- t=0:dt:time; subplot(2,2,1),plot(t,u(m,:),'k-'),grid,xlabel('時間(s)'),ylabel('位移(m)'),title('頂層位移的時程曲