立即咨询有惊喜哦 !

logo
logo
简介
简介
联系方式
正脉科工有限元分析联系电话010-81387990
邮箱xyq@vipstq.com
 
 

结构动力分析频域分析法——模态叠加法

作者:管理员    发布于:2014-12-22 02:10:01    文字:【】【】【

说明:

%a:频域内的模态分析法基本原理可以参考相关文献

%b:程序输入的参数,在前几行也基本上数目,质量矩阵和刚度矩阵全部来自于ANSYS有限元模型

%c:本程序目前只是用beam44,beam18x基本单元构成的有限元模型动力分析计算

%d:本程序目前得到脉动响应结果偏小,原因未找到

%e:希望大家找出其中的原因,一起来讨论

clc;
clear;
flag=input('(SSRS:
输入0CQC:输入
1):');
fs=input('
样本采样频率:
');
dir_in=uigetfile({'*.txt';'*.dat'},'
输入风荷载文件名
');
dof_in=uigetfile({'*.txt';'*.dat'},'
节点约束矩阵(Size:7*..
..*7)');
mass_in=uigetfile({'*.txt';'*.dat'},'
质量矩阵
');
stiff_in=uigetfile({'*.txt';'*.dat'},'
刚度矩阵
');
nodeL_in=uigetfile({'*.txt';'*.dat'},'
加载节点编号
');
node_degree=load(dof_in);
[row,col]=size(node_degree);
if row==7
    node_degree=node_degree';%
转置成node_sum*7的矩阵

    node_sum=col;
else
    node_sum=row;
end
mass_array=load(mass_in);
stifness_array=load(stiff_in);
DOF_sum=sum(sum(node_degree(:,2:7)>0));%
计算总的自由度
mass_array=reshape(mass_array,DOF_sum,DOF_sum);
stifness_array=reshape(stifness_array,DOF_sum,DOF_sum);
%-------------------------------------------------------------------------
%----------------------
导入荷载时程----------------------------------------
node_load=load(nodeL_in);
load_L=length(node_load);
nstep=1000;
%
与节点位移方向对应

load_dof=zeros(6,node_sum,nstep);
fid=fopen(dir_in,'r');
for ii=1:nstep
    data_in=fscanf(fid,'%g',[3 load_L]);
    load_dof(1:3,node_load,ii)=data_in;
end
fclose(fid);
%-------------------------------------------------------------------------
count=0;
load_F=zeros(DOF_sum,nstep);
mean_F=zeros(DOF_sum,1);
std_F=zeros(DOF_sum,1);
for ii=1:node_sum
    for jj=2:7
        if node_degree(ii,jj)>0.5
         count=count+1;
         if jj<5  %
只读XYZ方向的力
         load_F(count,:)=load_dof(jj-1,ii,:);
         mean_F(count)=mean(load_F(count,:));
         std_F(count)=std(load_F(count,:));
         load_F(count,:)=load_F(count,:)-mean_F(count);
         end
        end     
    end
end
%-------------------------------------------------------------------------
%
静风响应和背景响应
y_mean=stifness_array\mean_F;
y_std=stifness_array\std_F;
out=[y_mean y_std];
save
静风背景响应.txt out -ascii
%-----------------------
生成广义荷载谱
-------------------------------------
%
A:自振频率:natural frequency;振型:
mode of vibration
[mode_all,freq_w]=eig(stifness_array,mass_array);%
求结构振型和自振频率

freq_w=diag(freq_w);
LL=freq_w>0;mode_all=mode_all(:,LL);freq_w=freq_w(LL);
[freq_w,Loc]=sort(freq_w,'ascend');%
频率自从低到高排列
mode_all=mode_all(:,Loc);
%N=sum(freq_w<(pi*fs)^2);%
这里有待商榷,计算前N阶模态
N=100;
mode_all=mode_all(:,1:N);
%
B:确定广义力(generalized forces
force_g=mode_all'*load_F;
%
C:生成广义荷载功率谱
%
采用CPSD函数求互功率谱(Cross power spectral density),窗选用矩形窗rectwin
fn_length=floor(nstep/2)+1;
cross_psd=zeros(N,N,fn_length);
if flag
  for ii=1:N
    for jj=1:ii
        [Pxy,f]=cpsd(force_g(ii,:),force_g(jj,:),rectwin(nstep),0,nstep,fs);
        cross_psd(ii,jj,:)=Pxy;
        cross_psd(jj,ii,:)=conj(Pxy);
    end
 end
else
 for ii=1:N
    [Pxy,f]=cpsd(force_g(ii,:),force_g(ii,:),rectwin(nstep),0,nstep,fs);
    cross_psd(ii,ii,:)=Pxy;
 end
end
pause;
%-------------------------------------------------------------------------
clear load_F load_dof force_g
%-----------------------
确定频域函数
---------------------------------------
%
A:广义刚度矩阵和广义质量矩阵

M_g=mode_all'*mass_array*mode_all;
K_g=mode_all'*stifness_array*mode_all;
theta=0.02;%
各阶阻尼比
C=2*theta*sqrt(M_g.*K_g);
cross_yy=zeros(DOF_sum,fn_length);
for ii=1:fn_length
    w=f(ii)*2*pi;
    FrF=K_g+(1i)*w*C-w*w*M_g;
    %FrF=inv(FrF);  %
频响函数
    H_frf=zeros(N,N);
    for jj=1:N
        H_frf(jj,jj)=1/FrF(jj,jj);
    end
    %
位移响应互功率谱
    Sqq=H_frf*cross_psd(:,:,ii)*conj(H_frf);
    for jj=1:DOF_sum
    cross_yy(jj,ii)=mode_all(jj,:)*Sqq*mode_all(jj,:)';
    end
end

yy_std=f(2).*sum(abs(cross_yy'));
yy_std=sqrt(yy_std);

 

 
 
联系我们
 
 
 
图片
图片
多行文字
电话:010-81387990
手机:18301320667
联系人:李老师
网址:www.101315.com
邮箱:bjxu@vip.163.com
地址:北京市房山区良乡绿地启航国际商务办公区14号楼西楼5层
脚注信息
北京正脉科工科技有限公司  电话:010-81387990                            关注微信了解更多信息
版权所有 Copyright(C)2009-2013 京ICP备14038777号-1                                                                    
图片
图片