北京工業(yè)大學建筑工程學院 結構動力學中心差分法求解單自由度體系的自由振動問題前言時域逐步積分法是根據(jù)運動方程，引進某些假設，建立由t時刻狀態(tài)向量、到t時刻的狀態(tài)向量、的遞推關系，從而從t0時刻的初始狀態(tài)向量、出發(fā)，逐步求出各時刻的狀態(tài)向量，由于引進的假設條件不同，可以有各種不同的方法，下面主要介紹一種時域逐步積分方法中心差分法。中心差分法(central difference method)原理1 中心差分法的基本思路：是將運動方程中的速度向量和加速度向量用位移的某種組合來表示，將微分方程組的求解問題轉化為代數(shù)方程組的求解問題，并在時間區(qū)間內(nèi)求得每個微小時間區(qū)間的遞推公式，進而求得整個時程的反應。中心差分法只在相隔一些離散的時間區(qū)間內(nèi)滿足運動方程，其基于有限差分代替位移對時間的求導（即速度和加速度），如果采用等時間步長，則速度與加速度的中心差分近似為： (a) (b)而離散時間點的運動為 ( 0,1,2,3,)由體系運動方程為： (c)將速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到時刻的運動方程： （d）在（d）式中，假設和是已知的，即在及以前時刻的運動已知，則可以把已知項移到方程的右邊，整理得到： (e)由式（e）就可以根據(jù)及以前時刻的運動，求得時刻的運動，如果需要可以用式（a）和式（b）求得體系的速度和加速度。假設給定的初始條件為 （g）由式（g）確定。在零時刻速度和加速度的中心差分公式為： （h） （i）將式（i）消去得： （j） 而零時刻的加速度值可以用t0時的運動方程 確定即 （k）這樣就可以根據(jù)初始條件和初始荷載，就可以根據(jù)上式確定的值。下面給出采用中心差分法分析時的具體計算步驟：（1） 基本數(shù)據(jù)準備和初始條件計算 （2） 計算等效剛度和中心差分計算公式中的相關系數(shù) （3） 根據(jù)及以前時刻的運動，計算時刻的運動 如果需要，可計算（4）下一步計算用i+1代替i，對于線彈性結構體系，重復第3步，對于非線性結構體系，重復第2步和第3步。以上為中心差分法逐步計算公式，其具有2階精度，即誤差；并且為有條件穩(wěn)定，穩(wěn)定條件為： 上式中，為結構的自振周期，對于多自由度結構體系則為結構的最小自振周期。算例對于一個單層框架結構，假設樓板剛度無限大，且結構質(zhì)量集中于樓層，其質(zhì)量M=2000kg、剛度K50KN/m、阻尼系數(shù)C3KNs/m，假設結構處于線彈性狀態(tài)，用中心差分法計算結構的自由振動反應。采用MATLAB語言編程，并以單自由度體系為例進行計算,設初位移u00和初速度v0=0，取不同的步長分別計算，以驗證中心差分法的穩(wěn)定條件。先計算，由穩(wěn)定條件，而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) % 本程序采用中心差分法計算結構的動力響應% 本程序是既可以計算單自由度體系又可以計算多自由度體系，且均假設結構體系處于線彈性狀態(tài)；% -%輸入?yún)?shù)%-% M-質(zhì)量矩陣% 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(t-dt)%- 等效剛度Ke等效荷載Pe和相關系數(shù)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 zhouhuaping(S)-clear allM=input(輸入質(zhì)量矩陣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; %計算步數(shù)u=zeros(m,floor(n)+1); %設定存儲位移矩陣 v=zeros(m,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及系數(shù)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)/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(頂層位移的時程曲線); subplot(2,2,2),plot(t,v(m,:),r-),grid,xlabel(時間(s),ylabel(速度(m/s),title(頂層速度的時程曲線); subplot(2,2,3),plot(t,ac(m,:),b-),grid,xlabel(時間(s),ylabel(加速度(m/s2),title(頂層加速度的時程曲線); %-end運行centraldifferent.M文件輸入?yún)?shù)：K=50000; M=2000; C=3000; u0=0; v0=0;time20s；dt？中心差分法計算結果穩(wěn)定性分析由以上時程圖可以得到當0.1, 0.3, 0.4時逐步計算結果給出的結構運動趨向收斂的，即計算結果是穩(wěn)定的；當0.41,0.42, 0.45時逐步計算結果給出的結構運動趨向發(fā)散的，即結果是不穩(wěn)定的，且隨著步長的增加，計算結果發(fā)散