中心差分法計算單自由度體系的動力反應
北京工業大學結構工程
組員:胡建華 s201204111
馬恆 s201204112
陳相家s201204083
張力嘉s201204022
0前言時域逐步積分法是數值分析方法，它只假設結構本構關係在乙個微小的時間步距內是線性的。時域逐步積分法是結構動力問題中乙個得到廣泛研究的課題，並在結構動力反應計算中得到廣泛應用。由於引進的假設條件不同，可以有各種不同的方法，比如中心差分法，線性加速度法，平均常加速度法，wilson-θ法等，其中中心差分法精度最高。
在本文中，通過編制中心差分法計算單自由度體系的動力反應的程式來了解其應用及穩定性。
1中心差分法原理
中心差分法的基本思路：是將運動方程中的速度向量和加速度向量用位移的某種組合來表示，將微分方程組的求解問題轉化為代數方程組的求解問題，並在時間區間內求得每個微小時間區間的遞推公式，進而求得整個時程的反應。
中心差分法只在相隔一些離散的時間區間內滿足運動方程，其基於有限差分代替位移對時間的求導（即速度和加速度），如果採用等時間步長，，則速度與加速度的中心差分近似為：
a)b)
而離散時間點的運動為
0,1,2,3,……)
由體系運動方程為c)
將速度和加速度的差分近似公式（a）和式（b）代入式（c）可以得到時刻的運動方程：
d）在（d）式中，假設和是已知的，即在及以前時刻的運動已知，則可以把已知項移到方程的右邊，整理得到：
(e)由式（e）就可以根據及以前時刻的運動，求得時刻的運動，如果需要可以用式（a）和式（b）求得體系的速度和加速度。
假設給定的初始條件為
g）由式（g）確定。在零時刻速度和加速度的中心差分公式為：
h）i）
將式（i）消去得j
而零時刻的加速度值可以用t＝0時的運動方程確定
即k）這樣就可以根據初始條件和初始荷載，就可以根據上式確定的值。
2中心差分法程式編制原理：
（1） 基本資料準備和初始條件計算
（2） 計算等效剛度和中心差分計算公式中的相關係數
（3） 根據及以前時刻的運動，計算時刻的運動
如果需要，可計算
（4）下一步計算用i+1代替i，對於線彈性結構體系，重複第3步，對於非線性結構體系，重複第2步和第3步。
以上為中心差分法逐步計算公式，其具有2階精度，即誤差；並且為有條件穩定，穩定條件為：
上式中，為結構的自振週期，對於多自由度結構體系則為結構的最小自振週期。
3算例對於乙個單層框架結構，假設樓板剛度無限大，且結構質量集中於樓層，其質量m=2000kg、剛度k＝50kn/m、阻尼係數c＝3kns/m，假設結構處於線彈性狀態，用中心差分法計算結構的自由振動反應。
採用matlab語言程式設計，並以單自由度體系為例進行計算,設初位移u0＝0和初速度v0=0，取不同的步長分別計算，以驗證中心差分法的穩定條件。
先計算，由穩定條件，而rad/s,
則所以本次計算取＝0.1, 0.3, 0.4, 0.41, 0.42, 0.45分別進行計算
4 matlab程式清單
function [u,v,ac]=centraldifferent(m,c,k,u0,v0,time,dt)
% 本程式採用中心差分法計算結構的動力響應
% 本程式是既可以計算單自由度體系又可以計算多自由度體系，且均假設結構體系處於線彈性狀態；
輸入引數
% m質量矩陣
% c阻尼矩陣
% k剛度矩陣
% u0初始位移
% v0初始速度
% time---------模擬時間
% dt時間步長
輸出值% u位移
% v速度
% ac加速度
中心差分法主要公式及原理
% mx''+cx'+kx=0
% m*(x(t+dt)-2*x(t)+x(t-dt))/(dt^2)+c*(x(t+dt)-x(t-dt))/(2*dt)+k*x(t)=0
% (m/dt^2+c/2*dt)*(x(t+dt))=-(k-2*m/dt^2)*x(t+dt)-(m/dt^2-c/2*dt)*x(t-dt)
等效剛度ke等效荷載pe和相關係數a,b
% ke=m/dt^2+c/2*dt
% a=k-2*m/dt^2
% b=m/dt^2-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))/(dt^2)
初始條件
% x0''=(-c*x0'-k*x0)
% x(-1)=x0-x0'*dt+x0''*(dt^2)/2
copyright by zhouhuaping(s201004232)-----
clear all
m=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,floor(n)+1設定儲存速度矩陣
ac=zeros(m,floor(n)+1設定儲存加速度矩陣
p=zeros(m,floor(n)+1設定儲存荷載矩陣
u(:,2)=u0給定初位移
v(:,2)=v0給定初速度
ke=m/(dt^2)+((c)/(2*dt等效剛度ke及係數a、b
a=k-2*m/dt^2;
b=m/dt^2-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)*(dt^2))/2; %計算(0-dt)時刻位移
pe= -a*u(:,i-1)-b*u(:,i-2計算等效荷載pe
u(:,i)=ke\pe計算位移
v(:,i)=(u(:,i)-u(:,i-2))/(2*dt計算速度
ac(:,i)=(u(:,i)-2*u(:,i-1)+u(:,i-2))/(dt^2); %計算加速度
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/s^2)'),title('頂層加速度的時程曲線');
end執行檔案
輸入引數：
k=50000; m=2000; c=3000; u0=0; v0=0;time＝20s；dt＝？
5中心差分法計算結果穩定性分析
由以上時程圖可以得到當＝0.1, 0.3, 0.
4時逐步計算結果給出的結構運動趨向收斂的，即計算結果是穩定的；當＝0.41,0.42, 0.
45時逐步計算結果給出的結構運動趨向發散的，即結果是不穩定的，且隨著步長的增加，計算結果發散得越來越快。
工程研究中心申請報告編制提綱
附件3 一 摘要 左右 二 建設背景及必要性 一 本領域在國民經濟建設中的地位與作用。二 國內外技術和產業發展狀況 趨勢與市場分析。三 本領域當前急待解決的關鍵技術問題。四 本領域成果轉化與產業化存在的主要問題及原因。五 建設工程中心的意義與作用。三 申報單位概況和建設條件 一 申報單位及主要發起單...
湖南省工程研究中心組建方案編制提綱
附件1 一 工程中心組建方案摘要 左右 1 專案名稱 2 組建單位法人及自然人情況 3 經營戰略與經營計畫 4 建設內容 規模 方案和地點 5 主要建設條件 6 專案總投資 投資構成及資金籌措方案 7 經濟分析與主要技術經濟指標 二 工程中心建設背景及必要性 1 本領域在國民經濟建設中的地位與作用 ...
風電場工程可行性研究報告編制辦法
第一章總則 第一條為了統一風電場工程可行性研究報告訴編制的原則 內容 深度和技術要求，特制定 風電場工程可行性研究報告編制辦法 以下簡稱本辦法 第二條本辦法適用於規劃建設的風電場工程專案。第二章編制依據和深度 第三條進行可行性研究工作時應對風電場工程的建設條件進行深入調查，取得可靠的基礎資料。收集的...