1、精選優質文檔-傾情為你奉上結構動力學大作業班級： 學號：姓名： 目 錄1. Wilson-法原理簡介圖1-1 Wilson-法示意圖Wilson-法是基于對加速度a的插值近似得到的，圖1-1為Wilson-法的原理示意圖。推導由t時刻的狀態求t+t時刻的狀態的遞推公式： (1-1)對積分可得速度與位移的表達式如下： (1-2) (1-3)其中=t，由式(1-2)、(1-3)可以解出： (1-4) (1-5)將式(1-4)、(1-5)帶入運動方程： (1-6) (1-7)注意到此時的式子為和上一個時刻、以及t+t時刻的荷載相關，可以運用迭代的思想來求解，下圖給出線彈性條件下Wilson-法的流程

2、圖：圖1-2 Wilson-法流程圖2. Wilson-程序驗算對線彈性條件下的Wilson-法進行MATLAB編程，源代碼見附錄。選取如下算例進行驗證。對于一個單自由度的無阻尼結構，當其受到一個周期荷載時，其結構響應分為穩態解和瞬態解，由于沒有阻尼的影響，其瞬態解并不會衰減，其理論表達式為： (2-1)式中，為位移響應，為激勵，為剛度，為荷載頻率與固有振動頻率之比，為荷載頻率，為結構固有頻率。現令為1，為1，則為1，取為2/3。程序求得的解與解析解對比如圖2-1所示（由于理論解與程序基本重合，所以將理論解乘以-1，方便比較）： a)位移 b)速度c)加速度圖2-1 Wilson-法結果驗證2

3、.1 t的影響上述算例驗證時選擇的t非常小，因此看不出理論解與Wilson-法的求解區別，以下改變t的取值，探討t對迭代的影響。圖2-2 t對位移曲線的影響可以看出并不是t太大時計算結果很不準確，偏小，反映不出周期特征；當t合適時正好基本和理論解重合，也不是t越小越好越小時越能反映出一些細部特征，但這也不是很準確。2.2 的影響當>1.37時，該算法是無條件穩定的算法，以下探討對算法的影響。圖2-3 對位移曲線的影響由上圖可知隨著值越大，位移的周期變大。3. 非線性問題求解由于實際結構并不一定為線性，其剛度會隨著位移的的變化而改變，下圖為求解非線性問題時的Wilson-法流程。此處要說明

4、的是，剛度矩陣Ky(t)是與位移相關的量，判斷那時候速度的大小是為了確定其是否處于卸載段。具體可能得根據實際情況求解。圖3-1 Wilson-法解非線性問題修改MATLAB程序，并用該程序來計算如下例題：對該問題采用Wilson-法非線性方式計算，采用t=0.1s和t=0.05s兩種方式，計算位移、速度和加速度曲線如下圖所示： a)位移 b)速度 c)加速度圖3-2 非線性分析結果由上圖可知，結構在0.6s時達到位移極值，在t=0.1s和t=0.05s算得的值分別為0.096m和0.108m，速度極值在0.9s取到分別為-0.468和-0.580，加速度極值在t=0.1s時為0.7s時取到，為

5、-2.127，在t=0.05s極值在0.75s時取到，為-2.9。4. 附錄Wilson-法源程序function y_1,y_2,y_3=wilson_theta(p,m,c,k,dt,v0,y0,a_0,theta)%p代表輸入的荷載，c為阻尼矩陣，dt為時間間隔,m為質量矩陣,k為剛度矩陣%v0為初始的速度，y0為初始的位移.a_0為初始加速度%輸出的矩陣y_1代表位移,y_2代表速度，y_3代表加速度if nargin<9 theta=1.4;endL,r=size(p);y_1=NaN(L,r);y_2=NaN(L,r);y_3=NaN(L,r);y_1(:,1)=y0;y_2

6、(:,1)=v0;y_3(:,1)=a_0;%計算積分常數a0=6/(theta*dt)2);a1=3/theta/dt;a2=2*a1;a3=theta*dt/2;a4=a0/theta;a5=-a2/theta;a6=1-3/theta;a7=dt/2;a8=dt2/6;%計算擬剛度矩陣k0=k+a0*m+a1*c;%計算擬荷載for i=1:r-1 R=p(:,i)+theta*(p(:,i+1)-p(:,i)+m*(a0*y_1(:,i)+a2*y_2(:,i)+2*y_3(:,i)+c*(a1*y_1(:,i)+2*y_2(:,i)+a3*y_3(:,i); y_theta=k0R;

7、 y_3(:,i+1)=a4*(y_theta-y_1(:,i)+a5*y_2(:,i)+a6*y_3(:,i); y_2(:,i+1)=y_2(:,i)+a7*(y_3(:,i)+y_3(:,i+1); y_1(:,i+1)=y_1(:,i)+y_2(:,i)*dt+a8*(y_3(:,i+1)+2*y_3(:,i);end上述代碼只適合分析線彈性結構，對于非線性結構，編寫起來比較繁瑣，針對不同的情況可能需要具體處理，所以本文只給出了針對本文例題的代碼。與上述代碼不同的是以下代碼增加了一個判斷的語句。function y_1,y_2,y_3=wilson_theta2(p,m,c,dt,v0

8、,y0,a_0,theta)%p代表輸入的荷載，c為阻尼矩陣，dt為時間間隔,m為質量矩陣%v0為初始的速度，y0為初始的位移.a_0為初始加速度%輸出的矩陣y_1代表位移,y_2代表速度，y_3代表加速度if nargin<9 theta=1.4;endL,r=size(p);y_1=NaN(L,r);y_2=NaN(L,r);y_3=NaN(L,r);y_1(:,1)=y0;y_2(:,1)=v0;y_3(:,1)=a_0;%計算積分常數a0=6/(theta*dt)2);a1=3/theta/dt;a2=2*a1;a3=theta*dt/2;a4=a0/theta;a5=-a2/t

9、heta;a6=1-3/theta;a7=dt/2;a8=dt2/6;for i=1:r-1 if y_2(:,i)>0 k=60*(y_1(:,i)<=0.05)+3/y_1(:,i)*(y_1(:,i)>0.05); else k=60; end %計算擬剛度矩陣 k0=k+a0*m+a1*c; %計算擬荷載 R=p(:,i)+theta*(p(:,i+1)-p(:,i)+m*(a0*y_1(:,i)+a2*y_2(:,i)+2*y_3(:,i)+c*(a1*y_1(:,i)+2*y_2(:,i)+a3*y_3(:,i); y_theta=k0R; y_3(:,i+1)=a4*(y_theta-y_1(: