虚拟激励法程序？？？？

suss2009 · 发表于 2010-11-10 13:27

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

那位兄台有虚拟激励法的程序？？？？最好是从有限元里提取信息的？可否分享一下。谢了！

yejet · 发表于 2010-11-12 06:33

clear;
clc;
n=4;
%虚数单位
II=sqrt(-1);
%结构质量矩阵
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);
end
zk(n,n)=K(n,n);
K=zk;
[tzxl,tzz]=eig(zk,M);
d=diag(sqrt(tzz));
%阻尼矩阵
%振型规一化
for i=1:n
[tzz1(i),j]=min(d);
tzxl1(:,i)=tzxl(:,j);
d(j)=max(d)+1;
end
%振型归一化取第一层的振型为1
for j=1:n
tzxl1(:,j)=tzxl1(:,j)/tzxl1(1,j);
end
w0=tzz1;
w=tzz1/(2*pi);
%振型zhx
zhx=tzxl1;
%阻尼比
ks0=0.05;
%各阶模态阻尼比均为0.05
ks=ones(n,1)*ks0;
%求广义质量
Mn=zhx'*M*zhx;
%求广义阻尼矩阵
C=zeros(n);
for i=1:n
C(i,i)=2*ks(i)*w0(i)*Mn(i,i);
end
%通过振型和广义阻尼矩阵求解阻尼矩阵
c=(zhx')\C/zhx;
%过滤白噪声参数
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);
%频率点数循环变量wn
wn=1;
%对频率进行循环，求解各频率点的时间历程
for w=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;
for kkk=1:maxt/dt
%求确定频率下各时间点的功率谱
Yw=Ywt(:,kkk);
y1=conj(Yw);
y2=transpose(Yw);
%确定时间点确定频率下的功率谱
Syyw=y1*y2;
for kk=1:n
Pwt(kk,kkk,wn)=Syyw(kk,kk);
end
end
wn=wn+1;
end
%求解完成后实际上wn为maxw/dw+1
%各层的时变方差，循环变量为：层数，时间点
Fangcha=zeros(n,maxt/dt);
for tn=1:maxt/dt
%求解各层的时变方差
for kk=1:n
xx1=zeros(wn-1,1);
%每一个时刻的方差对各频率点进行积分
for wn0=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);
A=[Z,I;-m\k,-m\c];
AF=-1*[zeros(n,1);ones(n,1)];
%采用科茨公式求解非齐次项
%求解积分点的矩阵指数
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);
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=(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=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));
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*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);
UK1=TUK+TKTK1;
disp(:,i)=UK1(1:n,:);
velp(:,i)=UK1(n+1:2*n,:);
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;
%计算指数矩阵T
for i=1:20
Ta=2*Ta+Ta^2;
end
T=I+Ta;

复制代码

suss2009 · 发表于 2010-11-18 11:17

谢谢你哥们

zhangqiao · 发表于 2011-1-7 11:09

怎么又是时程又是频域的，实现起来还挺复杂的。怀疑他是否真的像说的那样高效。一般的方法只是频域积分就可以了。还希望明白的人说说

linzhidan · 发表于 2011-1-11 13:48

林家浩老师刚来我们所讲虚拟激励法啊
正缺这个啊
谢谢

linzhidan · 发表于 2011-1-11 13:58

精细直接积分法的积分方法选择.doc
打开出现乱码，请楼主重新上传
谢谢

overmatch · 发表于 2011-1-16 18:26

就是，打开时乱码啊

secondye · 发表于 2011-1-17 01:12

呵呵，明天早上林家浩正要来我们这儿开个讲座。。。

jimo.787 · 发表于 2011-5-9 10:17

正在学习中
楼主有没关于结构失效的程序啊

lmf001 · 发表于 2011-6-17 16:44

非常感谢，正需要！

lmf001 · 发表于 2011-6-23 12:46

但是貌似不准确啊!

sysuzzd · 发表于 2012-3-19 07:41

回复 6 # linzhidan 的帖子

看了一下，要用pdf打开的。

linzhidan · 发表于 2012-4-13 16:48

谢谢了，我转化过来了

张如林 · 发表于 2012-5-15 17:13

好东西，谢谢，学习一下

jianan00 · 发表于 2012-9-28 10:19

好东西，谢谢楼主分享！！！

账号		自动登录	找回密码
密码			我要加入

[随机振动] 虚拟激励法程序？？？？

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐:

评分