计算平板振动模态的程序
本帖最后由 shwning 于 2012-12-8 19:54 编辑本人编的计算平板振动模态的程序，见附件。（部分参考MATLAB有限元结构动力学分析与工程应用），可计算，用到的子函数已给出。计算得到频率域理论值相差很大，振型很好。计算采用集中质量，一致质量矩阵哪位源程序课发给我，谢谢！同时请求找出计算频率过高的原因。
结构理论频率解：
1
2
3
4
5
6
(1,1)
(1,2)
(2,2)
(1,3)
(2,3)
(1,4)
98
245.7
393.1
491.4
638.8
835.4
% 静力学与固有特性分析
clear
clc
Optmass=1;
alpha=0.01;beta=0.02; %比例阻尼：C=alpha*M+beta*K
% input basic data
no_elL=10; %长度方向单元数:10
no_elW=10; %宽度方向单元数:10
no_el=no_elL*no_elW;%单元数
no_nel=4; %每个单元的节点数
no_dof=3; %每个节点自由度数
no_node=(10+1)*(10+1);%系统的节点总数
sys_dof=no_node*no_dof;
% 材料和几何参数
prop(1)=2.1e11; %杨氏模量
prop(2)=0.3; %泊松比
prop(3)=7860; %密度
prop(4)=0.5; %板长度
prop(5)=0.5; %板宽度
prop(6)=0.005; %板厚度
% 节点坐标
xn=zeros(no_node,1); yn=zeros(no_node,1); zn=zeros(no_node,1);
dx=prop(4)/no_elL;dy=prop(5)/no_elW;
k=0;
for i=1:(no_elW+1)
for j=1:(no_elL+1)
yn(j+k*(no_elL+1))=(i-1)*dy; %节点y轴坐标
xn(j+k*(no_elL+1))=(j-1)*dx; %节点x轴坐标
end
k=k+1;
end
gcoord=;
% 节点连接关系
ndconn=zeros(no_el,no_nel+1);
m=1;n=0;
for i=1:no_elW
for j=1:no_elL
ndconn(m,1)=m;
ndconn(m,2)=m+n;
ndconn(m,3)=m+n+1;
ndconn(m,4)=m+n+1+11;
ndconn(m,5)=m+n+11;
m=m+1;
end
n=n+1;
end
% 加载
% ======静载荷
P=;
% 边界条件
for i=1:10
connode(i,:)=;
connode(i+10,:)=;
connode(i+20,:)=;
connode(i+30,:)=;
end
for i=1:10
conval(i,:)=;
conval(i+10,:)=;
conval(i+20,:)=;
conval(i+30,:)=;
end
%---------------------------------------------
% 初始化
k=zeros(no_nel*no_dof,no_nel*no_dof);
m=zeros(no_nel*no_dof,no_nel*no_dof); %单元质量刚度矩阵
K=zeros(sys_dof,sys_dof);
M=zeros(sys_dof,sys_dof); %系统质量刚度矩阵
F=zeros(sys_dof,1); %系统载荷向量
bcdof=zeros(sys_dof,1);
bcval=zeros(sys_dof,1);
index=zeros(no_nel*no_dof,1);
=quadEleStiffmat(prop,2); %单元刚度矩阵
for th=1:no_el
%--------------------------------------------
nd(1)=ndconn(th,2); %第th单元的节点1
nd(2)=ndconn(th,3); %第th单元的节点2
nd(3)=ndconn(th,4); %第th单元的节点3
nd(4)=ndconn(th,5); %第th单元的节点4
%--------------------------------------------------
index=femSysEleindex(nd,no_nel,no_dof); %获得第th单元在系统矩阵中的位置
K=femAssembmat(K,k,index); %将第th单元的刚度矩阵组装到系统矩阵中
M=femAssembmat(M,m,index); %将第th单元的质量矩阵组装到系统矩阵中
end
% 施加边界条件
=femKcheck(K,M,F,bcdof,bcval);
=femSysAppBC(K0,M0,F0,bcdof0,bcval0);
=femMcheck(K1,M1,F1,bcdof0,bcval0);
% 阻尼矩阵
C=alpha*M2+beta*K2;
% 静力学分析
dispplate=K2\F2;
dispp=PostDeal(dispplate,Mi2,Ki0,[],sys_dof);
dispP=reshape(dispp(1:3:end),11,11);
figure(42),surfc(dispP)
% 内力分布
F=K*dispp;
F=reshape(F(1:3:end),11,11);
figure(43),surfc(F)
% 模态分析
=eig(K2,M2);
=sort(diag(D));
omega=sqrt(lambda);
Homega=omega/2/pi
%模态振型
%集中质量矩阵
for i=1:7
Ys=PostDeal(V(:,i),Mi2,Ki0,[],sys_dof);
ys=Ys(1:3:end);
ys=reshape(ys,11,11);
figure(i),surf(ys)
end
%================================
function =femCalConst(ConNode,ConVal,No_dof,Sys_dof)
bcdof=zeros(Sys_dof,1); % initializing the vector bcdof
bcval=zeros(Sys_dof,1); % initializing the vector bcval
=size(ConNode); % calculate the constrained dofs
for ni=1:n1
ki=ConNode(ni,1);
bcdof((ki-1)*No_dof+(1:No_dof))=ConNode(ni,2:(No_dof+1));% the code for constrained dofs
% 施加约束取“1”
% 不施加约束取“0”
bcval((ki-1)*No_dof+(1:No_dof))=ConVal(ni,2:(No_dof+1)); % the value at constrained dofs
end
%==========================================
function =quadEleStiffmat(prop,masschoose)
% //////////////////////////////
% E：杨氏模量
% poisson：泊松比
% density：密度
% a：单元x方向长度/2
% b：单元x方向长度/2
% t：薄板厚度
% masschoose：1-协调质量矩阵；2-集中质量矩阵
% -------------------
% k：刚度矩阵
% m：质量矩阵
% ///////////////////////////////
%刚度矩阵
E=prop(1);
poisson=prop(2);
density=prop(3);
a=prop(4)/2;
b=prop(5)/2;
t=prop(6);
k=cell(4,4); %刚度矩阵
zeta=[-1 1 1 -1];%x方向
eta= [-1 -1 1 1];%y方向
a_b=(a/b)^2;b_a=(b/a)^2;
H=(1/60/a/b)*E*t*t*t/12/(1-poisson*poisson);
for i=1:4
for j=1:4
zeta0=zeta(i)*zeta(j); %x方向
eta0=eta(i)*eta(j); %y方向
%----------------------------------------------------------------
H11=3*H*(15*(b_a*zeta0+eta0*a_b)+(14-4*poisson+5*b_a+5*a_b)*zeta0*eta0);
H12=-3*H*b*((2+3*poisson+5*a_b)*zeta0*eta(i)+15*a_b*eta(i)+5*poisson*zeta0*eta(j));
H13=3*H*a*((2+3*poisson+5*b_a)*eta0*zeta(i)+15*b_a*zeta(i)+5*poisson*eta0*zeta(j));
%------------------------------------------------------------------
H21=-3*H*b*((2+3*poisson+5*a_b)*zeta0*eta(j)+15*a_b*eta(j)+5*poisson*zeta0*eta(i));
H22=H*b*b*(2*(1-poisson)*zeta0*(3+5*eta0)+5*a_b*(3+zeta0)*(3+eta0));
H23=-15*H*poisson*a*b*(zeta(i)+zeta(j))*(eta(i)+eta(j));
%-----------------------------------------
H31=3*H*a*((2+3*poisson+5*b_a)*eta0*zeta(j)+15*b_a*zeta(j)+5*poisson*eta0*zeta(i));
H32=-15*H*poisson*a*b*(zeta(i)+zeta(j))*(eta(i)+eta(j));
H33=H*a*a*(2*(1-poisson)*eta0*(3+5*zeta0)+5*b_a*(3+zeta0)*(3+eta0));
k{i,j}=;
end
end
k=cell2mat(k); %元胞数组转化为矩阵
%========================================
%质量矩阵
switch masschoose
case 1
%协调一致矩阵
disp('采用协调一致矩阵');
case 2
%略去转动自由度上的质量
disp('采用集中质量矩阵')
m=zeros(12,12);
w=a*b*t*density;
m(1,1)=w;
m(4,4)=w;
m(7,7)=w;
m(10,10)=w;
otherwise
disp('????????sorry, you only can input ''1-(协调一致矩阵) or 2-(单元质量矩阵)''');
m=[];
end
%==================================
function index=femSysEleindex(elecond,numnodelement,numdofnode)
% 单元与系统之间的关联节点向量
% elecond：与某单元相连的节点
% numnodelement：单元节点数
% numdofnode：节点自由度数
% -----------------
% index：单元各节点向量在系统矩阵中位置
k=0;
for i=1:numnodelement
start=(elecond(i)-1)*numdofnode;
for j=1:numdofnode
k=k+1;
index(k)=start+j;
end
end
%===========================
function =femAssembmat(KM,km,index)
% 质量矩阵和刚度矩阵的组装
% KM：组装后刚度矩阵和质量矩阵
% km：单元刚度矩阵和质量矩阵
% index：单元与系统关联向量
% %-------------------
% KM：组装后刚度矩阵和质量矩阵
%///////////////////////////////
numdofelement=length(index);
for i=1:numdofelement
for j=1:numdofelement
KM(index(i),index(j))=KM(index(i),index(j))+km(i,j);
end
end
%================================
function =femKcheck(kk,mm,ff,bcdof,bcval)
% 检查总体刚度矩阵KK的主元是否为零
%--------------------------------------------------------------------------
=size(kk);
jk=0;
for ii=1:sdof1 % loop for check of the zero main elements in kk
check=kk(ii,ii);
if check==0
jk=jk+1; % location of the zero main element
kki(jk)=ii; % storring the location of the zero main element
end
end
if jk~=0
disp('main elements of kk:') % display the zero main element
kkii=kki;
disp('equal zero')
end
if jk==0
kki=[];
disp('=======main elements are not equal zero')
end
%--------------------------------------------------------------------------
%(2) eliminate the columns and rows in kk, mm, ff, bcdof, bcval
%--------------------------------------------------------------------------
if jk~=0
for jj=jk:-1:1 % loop for moving the columns and rows associated with
% the zero main element in the system matrix equation
hh=sdof1-kki(jj);
hr=kki(jj);
for i=1:hh % exchanging the rows in the equation
kt=kk(hr+i-1,:);
mt=mm(hr+i-1,:);
ft=ff(hr+i-1,:);
bdt=bcdof(hr+i-1);
bvt=bcval(hr+i-1,:);
kk(hr+i-1,:)=kk(hr+i,:);
mm(hr+i-1,:)=mm(hr+i,:);
ff(hr+i-1,:)=ff(hr+i,:);
bcdof(hr+i-1)=bcdof(hr+i);
bcval(hr+i-1,:)=bcval(hr+i,:);
kk(hr+i,:)=kt;
mm(hr+i,:)=mt;
ff(hr+i,:)=ft;
bcdof(hr+i,:)=bdt;
bcval(hr+i,:)=bvt;
end
for j=1:hh % exchanging the columns in the equation
kt=kk(:,hr+j-1);
mt=mm(:,hr+j-1);
kk(:,hr+j-1)=kk(:,hr+j);
mm(:,hr+j-1)=mm(:,hr+j);
kk(:,hr+j)=kt;
mm(:,hr+j)=mt;
end
end
end
sdof=sdof1-jk;
K=; % eliminating the columns and rows of the kk
M=; % eliminating the columns and rows of the mm
F=; % eliminating the rows of the ff
bc_dof=; % eliminating the rows of the bcdof
bc_val=; % eliminating the rows of the bcval
%--------------------------------------------------------------------------
%===================================
function =femMcheck(kk,mm,ff,bcdof,bcval)
% 检查总体质量矩阵MM的主元是否为零
=size(kk);
jk=0;
for ii=1:sdof1 % loop for check of the zero main elements in mm
check=mm(ii,ii);
if check==0
jk=jk+1; % location of the zero main element
mmi(jk)=ii; % storing the location of the zero main element
end
end
if jk==0
mmi=[];
disp('=======main elements are not equal zero')
end
%--------------------------------------------------------------------------
%(2) eliminate the columns and rows in kk, mm, ff, bcdof, bcval
%--------------------------------------------------------------------------
if jk~=0
for jj=jk:-1:1 % loop for moving the columns and rows associated with
% the zero main element in the system matrix equation
hh=sdof1-mmi(jj);
hr=mmi(jj);
for i=1:hh % exchanging the rows in the equation
kt=kk(hr+i-1,:);
mt=mm(hr+i-1,:);
ft=ff(hr+i-1,:);
bdt=bcdof(hr+i-1);
bvt=bcval(hr+i-1);
kk(hr+i-1,:)=kk(hr+i,:);
mm(hr+i-1,:)=mm(hr+i,:);
ff(hr+i-1,:)=ff(hr+i,:);
bcdof(hr+i-1)=bcdof(hr+i);
bcval(hr+i-1)=bcval(hr+i);
kk(hr+i,:)=kt;
mm(hr+i,:)=mt;
ff(hr+i,:)=ft;
bcdof(hr+i)=bdt;
bcval(hr+i)=bvt;
end
for j=1:hh % exchanging the columns in the equation
kt=kk(:,hr+j-1);
mt=mm(:,hr+j-1);
kk(:,hr+j-1)=kk(:,hr+j);
mm(:,hr+j-1)=mm(:,hr+j);
kk(:,hr+j)=kt;
mm(:,hr+j)=mt;
end
end
end
sdof=sdof1-jk;
K=; % eliminating the columns and rows of the kk
M=; % eliminating the columns and rows of the mm
F=; % eliminating the rows of the ff
bc_dof=; % eliminating the rows of the bcdof
bc_val=; % eliminating the rows of the bcval
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function =femSysAppBC(K,M,F,bcdof,bcval)
ni=length(bcdof);
sdof=size(K);
for ii=1:ni
if bcdof(ii)==1
for j=1:sdof
K(ii,j)=0;
K(j,ii)=0;
M(ii,j)=0;
M(j,ii)=0;
F(j)=F(j)-bcval(ii)*K(j,ii);
end
K(ii,ii)=1;
K(ii)=bcval(ii);
end
end
%--------------------------------------------------------------------------
function Ys=PostDeal(X,Ch3,Ch2,Ch1,SysDof)
% 根据需要将计算结果还原到与原来节点对应位置；
% 应用工程中根据需要做适当修改
% 应该注意原始数据检查的顺序，因此在还原时应该按照相反地顺序进行
% ------------------------------
% Ys-还原后的数据
% --------------------------------
% X-需要还原的数据,以列形式输入
% Ki-刚度矩阵中被舍去的主元位置编号
% Mi-质量矩阵中被舍去的主元位置编号
% BCi-根据边界被舍去的位置编号
% SysDof-系统矩阵的大小
lengX=length(X);
lengCh3=length(Ch3);
lengCh2=length(Ch2);
lengCh1=length(Ch1);
X2=zeros((lengX+lengCh3),1);
Ch3=sort(Ch3);
i=1;j=1;
if lengCh3~=0
for k=1:(lengX+lengCh3)
if k==Ch3(j) & j<lengCh3;
X2(k)=0;
j=j+1;
elseif i<=lengX
X2(k)=X(i);
i=i+1;
end
end
else
X2=X;
end
X1=zeros((lengX+lengCh3+lengCh2),1);
Ch2=sort(Ch2);
i=1;j=1;
if lengCh2~=0
for k=1:(lengX+lengCh3+lengCh2)
if k==Ch2(j) & j<lengCh2;
X1(k)=0;
j=j+1;
elseif i<=(lengX+lengCh3)
X1(k)=X2(i);
i=i+1;
end
end
else
X1=X2;
end
Ys=zeros(SysDof,1);
Ch1=sort(Ch1);
i=1;j=1;
if lengCh1~=0
for k=1:SysDof
if k==Ch1(j) & j<lengCh1;
Ys(k)=0;
j=j+1;
elseif i<=(lengX+lengCh3+lengCh2)
Ys(k)=X1(i);
i=i+1;
end
end
else
Ys=X1;
end
找到了一致质量单元矩阵，计算结果是正确的；只能认为计算频率的差异来源于采用集中质量矩阵带来的误差。 %协调一致矩阵
disp('采用协调一致矩阵');
m=zeros(12,12);
mass=density*t*A/12600;
m=[1727 466*b -461*a 613 199*b 274*a 197 -116*b 116*a 613 -274*b -199*a;
466*b 160*b*b -126*a*b 199*b 80*b*b 84*a*b 116*b -60*b*b 56*a*b 274*b -120*b*b -84*a*b;
-461*a -126*a*b 160*a*a -274*a -84*a*b -120*a*a -116*a 56*a*b -60*a*a -199*a 84*a*b 80*a*a;
613 199*b -274*a 1727 461*b 461*a 613 -274*b 199*a 197 -116*b -116*a;
199*b 80*b*b -84*a*b 461*b 160*b*b 126*a*b 274*b -120*b*b 84*a*b 116*b -60*b*b -56*a*b;
274*a 84*a*b -120*a*a 461*a 126*a*b 160*a*a 199*a -84*a*b 80*a*a 116*a*a -56*a*b -60*a*a;
197 116*b -116*a 613 274*b 199*a 1727 -461*b 461*a 613 -199*b -274*a;
-116*b -60*b*b 56*a*b -274*b -120*b*b -84*a*b -461*b 160*b*b -126*a*b -199*b 80*b*b 84*a*b;
116*a 56*a*b -60*a*a 199*a 84*a*b 80*a*a 461*a -126*a*b 160*a*a 274*a -84*a*b -120*a*a;
613 274*b -199*a 197 116*b 116*a*a 613 -199*b 274*a 1727 -461*b -461*a;
-274*b -120*b*b 84*a*b -116*b -60*b*b -56*a*b -199*b 80*b*b -84*a*b -461*b 160*b*b 126*a*b;
-199*a -84*a*b 80*a*a -116*a -56*a*b -60*a*a -274*a 84*a*b -120*a*a -461*a 126*a*b 160*a*a];
m=mass*m; 头大啊 亲们 乃们说的啥 请问为什么是先检查刚度矩阵的主元是否为0，然后施加边界条件，再次检查质量矩阵的主元是否为0，这是什么意思，我是在徐斌书上看到这个的，他那个好像有点小小的错误吧 {:{39}:} 谢谢啦很详细 好好学习下 好人 好好学习下{:{10}:} 请问如果加的不是静载荷，是集中力载荷该如何编程呢 谢谢楼主分享。 根本看不懂啊!都是直接用ANSYS之类的软件做
页:
[1]