matlab在建筑抗震工程中的应用 [大家看看,上个类似的帖子有误]

江山美人

某三层钢筋混凝土结构，地震波采用200gal El Centro波，采样周期为0.02秒。

用振型分解法求解结构地震反应的MATLAB程序如下:

xs=2*0.287;                 % 地震波数据

dzhbo=load('ELCENTRO.txt', 'r');    

ag=dzhbo*0.01*xs;

dt=0.02;

ndzh=400;

cn=3;                      % 结构参数

m0=[2.762 2.760 2.300]*1e+3;

k0=[2.485 1.9211.522]*1e+5;

l=diag(ones(cn));

m=diag(m0);

[ik]=matrixju(k0,cn);

[x,d]=eig(ik,m);           % 结构动力特性的求解

d=diag(sqrt(d));

for i=1:cn

[d1(i),j]=min(d);

xgd(:,i)=x(:,j);

d(j)=max(d)+1;

end

w=d1;

x=xgd;

a1=2*w(1)*w(2)*(0.05*w(2)-0.07*w(1))/(w(2)^2-w(1)^2);

a2=2*(0.07*w(2)-0.05*w(1))/(w(2)^2-w(1)^2);

for j=1:cn

x(:,j)=x(:,j)/x(cn,j);

znb0(j)=(a1+a2*w(j)^2)/2/w(j);

zhcan(j)=(x(:,j))'*m*l/((x(:,j))'*m*x(:,j));            % 求解振型参与系数

[dlt(j,:),dltacceler(j,:)]=zxzj(znb0(j),w(j),ag);         % 求解 file:///C:/Users/SHISHI~1/AppData/Local/Temp/msohtmlclip1/01/clip_image002.gif和file:///C:/Users/SHISHI~1/AppData/Local/Temp/msohtmlclip1/01/clip_image004.gif

end

for i=1:cn                   % 求解结构各层的地震反应

disp1=0;

accel1=0;

for j=1:cn

disp0=zhcan(j)*dlt(j,:)*x(i,j);

accel0=zhcan(j)*dltacceler(j,:)*x(i,j);

disp1=disp1+disp0;

accel1=accel1+accel0;

end

disp(i,:)=disp1;

accel(i,:)=accel1;

end

t=0:dt:ndzh*dt;

subplot(2,2,1);

plot(t,disp(3,:)*1e+3,'k-');

subplot(2,2,2);

plot(t,accel(3,:),'k-');

% This sub-program is solvingthe dynamic responses of single degree system

function  [bx,acceler]=zxzj(znb,w,dag)

dt=0.02;

n=400;

x(1)=0;

dx(1)=0;

ddx(1)=0;

s=1+znb*dt*w+w^2*dt^2/6;     %中间参数 s

fori=1:n

a(i)=x(i)+dx(i)*dt+ddx(i)*dt^2/3;

b(i)=dx(i)+ddx(i)*dt/2;

ddx(i+1)=-1*(dag(i+1)+2*znb*w*b(i)+w^2*a(i))/s;    % 加速度

dx(i+1)=b(i)+ddx(i+1)*dt/2;                 % 速度

x(i+1)=a(i)+ddx(i+1)*dt^2/6;               % 位移

end

bx=x;

acceler=ddx;

% This sub-program is formatrix aggregation of system

function  [kcju]=matrixju(korc,cn)

kcju=zeros(cn);

for i=1:cn-1

kcju(i,i)=korc(i)+korc(i+1);

kcju(i,i+1)=-korc(i+1);

kcju(i+1,i)=-korc(i+1);

end

kcju(cn,cn)=korc(cn)

运行之后 出现下面提示 ??? Error using ==> loadUnable to read file ELCENTRO.txt: No such file or directory.还有蓝色这个函数matrixju不知道什么意思？

请各位大侠指教。

ChaChing

不是都说了

ELCENTRO.txt: No such file or directory

建议help下load

江山美人

谢谢！晓得了，是我没要把文件放在工作目录里面，随便放在桌面山了

江山美人

本人刚刚开始学习，以后请大侠多多指教

江山美人

我把最后两个子程序存为matrixju.m和zxzj.m两个m文件后，程序如下，但是所得结果与已知结果不同，郁闷啊！！请大侠再指教！！不胜感激！！

1. xs=2*0.287;                 % 地震波数据
   
2. dzhbo=load('ELCENTRO.txt');   
   
3. ag=dzhbo*0.01*xs;
   
4. dt=0.02;
   
5. ndzh=400;
   
6. cn=3;                      % 结构参数
   
7. m0=[2.762 2.760 2.300]*1e+3;
   
8. k0=[2.485 1.921 1.522]*1e+5;
   
9. l=diag(ones(cn));
   
10. m=diag(m0);
    
11. [ik]=matrixju(k0,cn);
    
12. [x,d]=eig(ik,m);           % 结构动力特性的求解
    
13. d=diag(sqrt(d));  
    
14. for i=1:cn
    
15. [d1(i),j]=min(d);
    
16. xgd(:,i)=x(:,j);
    
17. d(j)=max(d)+1;
    
18. end
    
19. w=d1;
    
20. x=xgd;
    
21. a1=2*w(1)*w(2)*(0.05*w(2)-0.07*w(1))/(w(2)^2-w(1)^2);
    
22. a2=2*(0.07*w(2)-0.05*w(1))/(w(2)^2-w(1)^2);
    
23. for j=1:cn
    
24. x(:,j)=x(:,j)/x(cn,j);
    
25. znb0(j)=(a1+a2*w(j)^2)/2/w(j);
    
26. zhcan(j)=(x(:,j))'*m*l/((x(:,j))'*m*x(:,j));           % 求解振型参与系数
    
27. [dlt(j,:),dltacceler(j,:)]=zxzj(znb0(j),w(j),ag);           % 求解 ∆_j和∆ ̈_j
    
28. end
    
29. for i=1:cn                   % 求解结构各层的地震反应
    
30. disp1=0;
    
31. accel1=0;
    
32. for j=1:cn
    
33. disp0=zhcan(j)*dlt(j,:)*x(i,j);
    
34. accel0=zhcan(j)*dltacceler(j,:)*x(i,j);
    
35. disp1=disp1+disp0;
    
36. accel1=accel1+accel0;
    
37. end
    
38. disp(i,:)=disp1;
    
39. accel(i,:)=accel1;
    
40. end
    
41. t=0:dt:ndzh*dt;
    
42. subplot(2,2,1)
    
43. plot(t,disp(3,:)*1e+3,'b-')
    
44. subplot(2,2,2)

复制代码