《结构的动力学方程》由会员分享，可在线阅读，更多相关《结构的动力学方程（8页珍藏版）》请在金锄头文库上搜索。
1、结构的动力学方程clear;clc;n=4;II=sqrt(-1); %复数i%主结构质量、阻尼、刚度矩阵M=eye(n)*1.0e+4; %对角阵K=eye(n)*1.6*1.0e+7;%主结构刚度矩阵聚合zk=zeros(n);for j=1:(n-1)zk(j,j)=K(j,j)+K(j+1,j+1);zk(j,j+1)=-K(j+1,j+1);zk(j+1,j)=-K(j+1,j+1);endzk(n,n)=K(n,n);k=zk;m=M;%此处我的结构模型的质量、刚度矩阵应该从其它有限元软件导出来。%求解各阶模态频率tzxl,tzz=eig(k,m);d=diag(sqrt(tzz)
2、;%振型规一化for i=1:n tzz1(i),j=min(d); tzxl1(:,i)=tzxl(:,j); d(j)=max(d)+1;end%振型归一化取第一层振型为1for j=1:n tzxl1(:,j)=tzxl1(:,j)/tzxl1(1,j);endw0=tzz1;w=tzz1/(2*pi);zhx=tzxl1;广义阻尼矩阵 各阶模态阻尼比都取%阻尼比ks0=0.05;ks=ones(n,1)*ks0;第n阶广义质量： %求广义质量Mn=zhx*m*zhx;阻尼矩阵为：%求阻尼矩阵C=zeros(n);for i=1:n C(i,i)=2*ks(i)*w0(i)*Mn(i,i
3、);endc=(zhx)C/zhx;参数eg即%过滤白噪声参数eg=0.6;wg=15.708;S0=0.001574;%按照书上的要求，取频率和时间的最大值和步长%频率间隔dw=0.3;%最大频率范围maxw=45;%最大时间值maxt=40;%时间间隔dt=0.2;%各层各时间点频率点的功率谱密度,循环变量：层数，时间点，频率点Pwt=zeros(n,maxt/dt,maxw/dw);%频率点数循环变量wnwn=1;%对频率进行循环，求解各频率点的时间历程for w=0:dw:maxwx1=1+4*eg2*(w/wg)2;x2=(1-(w/wg)2)2+4*eg2*(w/wg)2;Sgw=
4、x1*S0/x2;s=sqrt(Sgw);%采用精细积分法进行求解时间历程，得到位移和速度时程disp,velp=JINGXI67(M,zk,c,dt,maxt,w,s,n);Ywt=disp; for kkk=1:maxt/dt %求确定频率下各时间点的功率谱 Yw=Ywt(:,kkk); 每一时刻和频率点的位移向量，对其进行求共轭和装置得到协方差矩阵，对角上的元素即是每一时刻的各层的功率谱 y1=conj(Yw); %共轭 y2=transpose(Yw); %转置矩阵 %确定时间点确定频率下的功率谱Yw,取对角线元素 Syyw=y1*y2; for kk=1:n Pwt(kk,kkk,w
5、n)=Syyw(kk,kk); end endwn=wn+1;end%求解完成后实际上wn为maxw/dw+2，实际频率点个数为maxw/dw+1%各层的时变方差，循环变量为：层数，时间点Fangcha=zeros(n,maxt/dt);for tn=1:maxt/dt%求解各层的时变方差 for kk=1:n xx1=zeros(wn-1,1);%每一个时刻的方差对各频率点进行积分，频率点数取maxw/dw+1，即wn-1 for wn0=1:wn-1 xx1(wn0)=Pwt(kk,tn,wn0); end%采用复合梯形求积公式对功率谱进行积分得到方差 Fangcha(kk,tn)=(xx
6、1(1)+xx1(wn-1)+2*sum(xx1(2:wn-1-1)*dw; endend%画图c1=(1:maxt/dt)*dt;d1=Fangcha(1,:)/S0;d2=Fangcha(2,:)/S0;d3=Fangcha(3,:)/S0;d4=Fangcha(4,:)/S0;figure(3)plot(c1,d1,k,c1,d2,r,c1,d3,m,c1,d4,r-)精细积分的程序function disp,velp=JINGXI67(m,k,c,dt,maxt,w,s,n)%虚数单位II=sqrt(-1);% 中的IIW=II*w; I=eye(n);Z=zeros(n);离散化n自
7、由度结构受均匀调制演变随机激励时的运动微分方程可表示为：其中为平稳高斯白噪声随机过程向量，为调制函数。该方程可表示为其中I为单位矩阵，记，程序中的A即上面的矩阵A，AF即GA=Z,I;-mk,-mc;AF=-1*zeros(n,1);ones(n,1);对线性体系来说，方程为线性向量常微分方程其解为对上式进行离散化为时间步长矩阵T 可用数值方法精细算得采用函数EXPHDT计算，H即输入的矩阵A，dt即为时间间隔%采用科茨公式求解非齐次项%求解积分点的矩阵指数EXPDT=EXPHDT(A,dt);EXP0_75DT=EXPHDT(A,0.75*dt);EXP0_50DT=EXPHDT(A,0.5
8、0*dt);EXP0_25DT=EXPHDT(A,0.25*dt);第二项等价于单个函数的积分，采用科茨公式进行积分，具有5阶精度zeroq=zeros(n,maxt/dt);disp=zeroq; velp=zeroq; zerop=zeros(n,1);dispt=zerop;velpt=zerop;UK=dispt;velpt;for i=1:maxt/dt%给出各积分点的X坐标，即时间t=；t0_25=；t0_50=；t0_75=；t1=t=(i-1)*dt;t0_25=t+0.25*dt;t0_50=t+0.50*dt;t0_75=t+0.75*dt;t1=i*dt;%调制函数求解g
9、t=；gt0_25=；gt0_50=；gt0_75=；gt1=。gt=4*(exp(-0.0995*t)-exp(-0.199*t);gt0_25=4*(exp(-0.0995*t0_25)-exp(-0.199*t0_25);gt0_50=4*(exp(-0.0995*t0_50)-exp(-0.199*t0_50);gt0_75=4*(exp(-0.0995*t0_75)-exp(-0.199*t0_75);gt1=4*(exp(-0.0995*t1)-exp(-0.199*t1);即AF*exp(IIW*t)；即AF*exp(IIW*t0_25)FTK=AF*exp(IIW*t)*gt*
10、s;FT0_25K=AF*exp(IIW*t0_25)*gt0_25*s;FT0_50K=AF*exp(IIW*t0_50)*gt0_50*s;FT0_75K=AF*exp(IIW*t0_75)*gt0_75*s;FTK1=AF*exp(IIW*t1)*gt1*s;%即TUK= EXPDT*UK EXPDT即TUK=EXPDT*UK;%采用科茨公式求解非齐次项X1=7*EXPDT*FTK;X2=32*EXP0_75DT*FT0_25K;X3=12*EXP0_50DT*FT0_50K;X4=32*EXP0_25DT*FT0_75K;X5=7*FTK1;TKTK1=1/90*dt*(X1+X2+X3+X4+X5);即TKTK1 UK1=TUK+TKTK1;%保存各时刻的位移和速度向量disp(:,i)=UK1(1:n,:);velp(:,i)=UK1(n+1:2*n,:);%将UK1赋给UK进入下一轮循环UK=UK1;end求解function T=EXPHDT(A,dt)%计算矩阵指数N=20;m=2N;t=dt/m;I=eye(size(A);Ta=A*t+(A*t)2*(I+(A*t)/3+(A*t)2/12)/2;%计算指数矩阵Tfor i=1:20 Ta=2*Ta+Ta2;endT=I+Ta;最后的曲线是这样子，很明显和课本上对不上
