虚拟激励法程序求助.docx

资源描述

虚拟激励法程序求助.docx

《虚拟激励法程序求助.docx》由会员分享，可在线阅读，更多相关《虚拟激励法程序求助.docx（6页珍藏版）》请在冰点文库上搜索。

虚拟激励法程序求助.docx

虚拟激励法程序求助

hustdd发表于:

2008-12-1411:

34来源:

振动资讯

学习林家浩教授和张亚辉教授的专著《随机振动的虚拟激励法》时，编制的程序和书上的结果对不上，请大家帮忙看看。

clear;clc;n=4;

II=sqrt（-1）;

%主结构质量、阻尼、刚度矩阵

M=eye（n）*1.0e+4;

K=eye（n）*1.6*1.0e+7;

%主结构刚度矩阵聚合zk=zeros（n）;

forj=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））;

%振型规一化fori=1:

[tzz1（i）,j]=min（d）;

tzxl1（:

i）=tzxl（:

j）;d（j）=max（d）+1;

end

%振型归一化取第一层振型为1forj=1:

tzxl1（:

j）=tzxl1（:

j）/tzxl1（1,j）;end

w0=tzz1;w=tzz1/（2*pi）;zhx=tzxl1;

广义阻尼矩阵

各阶模态阻尼比都取

%阻尼比

ks0=0.05;

ks=ones（n,1）*ks0;

第n阶广义质量：

%求广义质量Mn=zhx'*m*zhx;阻尼矩阵为：

%求阻尼矩阵

C=zeros（n）;fori=1:

C（i,i）=2*ks（i）*w0（i）*Mn（i,i）;end

c=（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;

%对频率进行循环，求解各频率点的时间历程forw=0:

dw:

maxw

x1=1+4*eg^2*（w/wg）^2;

x2=（1-（w/wg）^2）^2+4*eg^2*（w/wg）^2;Sgw=x1*S0/x2;

s=sqrt（Sgw）;

%采用精细积分法进行求解时间历程，得到位移和速度时程[disp,velp]=JINGXI67（M,zk,c,dt,maxt,w,s,n）;

Ywt=disp;

forkkk=1:

maxt/dt

%求确定频率下各时间点的功率谱

Yw=Ywt（:

kkk）;

每一时刻和频率点的位移向量，对其进行求共轭和装置得到协方差矩阵，对角上的元素即是每一时刻的各层的功率谱

y1=conj（Yw）;y2=transpose（Yw）;

%确定时间点确定频率下的功率谱Yw,取对角线元素Syyw=y1*y2;

forkk=1:

nPwt（kk,kkk,wn）=Syyw（kk,kk）;end

endwn=wn+1;end

%求解完成后实际上wn为maxw/dw+2，实际频率点个数为maxw/dw+1

%各层的时变方差，循环变量为：

层数，时间点Fangcha=zeros（n,maxt/dt）;

fortn=1:

maxt/dt

%求解各层的时变方差forkk=1:

nxx1=zeros（wn-1,1）;

%每一个时刻的方差对各频率点进行积分，频率点数取maxw/dw+1，即wn-1forwn0=1:

wn-1

xx1（wn0）=Pwt（kk,tn,wn0）;end

%采用复合梯形求积公式对功率谱进行积分得到方差Fangcha（kk,tn）=（xx1

（1）+xx1（wn-1）+2*sum（xx1（2:

wn-1-1）））*dw;end

end

%画图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自由度结构受均匀调制演变随机激励时的运动微分方程可表示为：

其中为平稳高斯白噪声随机过程向量，为调制函数。

该方程可表示为

其中I为单位矩阵，记

，，

程序中的A即上面的矩阵A，AF即GA=[Z,I;-m\k,-m\c];

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.50*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];

fori=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;

%调制函数

求解gt=；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*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*UKEXPDT即

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=2^N;t=dt/m;I=eye（size（A））;

Ta=A*t+（A*t）^2*（I+（A*t）/3+（A*t）^2/12）/2;

%计算指数矩阵Tfori=1:

Ta=2*Ta+Ta^2;end

T=I+Ta;