马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
本帖最后由 shwning 于 2012-12-8 19:54 编辑
平板声辐射计算源程序（主函数plate_ex3）.rar
(7.51 KB, 下载次数: 49)
本人编的计算平板振动模态的程序，见附件。（部分参考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=[xn yn zn];
% 节点连接关系
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=[500,60,1];
% 边界条件
for i=1:10
connode(i,:)=[i,1,0,0];
connode(i+10,:)=[i*11,1,0,0];
connode(i+20,:)=[i*11+1,1,0,0];
connode(i+30,:)=[i+111,1,0,0];
end
for i=1:10
conval(i,:)=[i,0,0,0];
conval(i+10,:)=[i*11,0,0,0];
conval(i+20,:)=[i*11+1,0,0,0];
conval(i+30,:)=[i+111,0,0,0];
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);
[k,m]=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
% 施加边界条件
[K0,M0,F0,bcdof0,bcval0,sdof0,Ki0]=femKcheck(K,M,F,bcdof,bcval);
[K1,M1,F1]=femSysAppBC(K0,M0,F0,bcdof0,bcval0);
[K2,M2,F2,bcdof2,bcval2,sdof2,Mi2]=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)
% 模态分析
[V,D]=eig(K2,M2);
[lambda,ki]=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 [bcdof,bcval]=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
[n1,n2]=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 [k,m]=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}=[H11 H12 H13;H21 H22 H23;H31 H32 H33];
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 [KM]=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 [K,M,F,bc_dof,bc_val,sdof,kki]=femKcheck(kk,mm,ff,bcdof,bcval)
% 检查总体刚度矩阵KK的主元是否为零
%--------------------------------------------------------------------------
[sdof1,n1]=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=[kk(1:sdof,1:sdof)]; % eliminating the columns and rows of the kk
M=[mm(1:sdof,1:sdof)]; % eliminating the columns and rows of the mm
F=[ff(1:sdof,:)]; % eliminating the rows of the ff
bc_dof=[bcdof(1:sdof)]; % eliminating the rows of the bcdof
bc_val=[bcval(1:sdof,:)]; % eliminating the rows of the bcval
%--------------------------------------------------------------------------
%===================================
function [K,M,F,bc_dof,bc_val,sdof,mmi]=femMcheck(kk,mm,ff,bcdof,bcval)
% 检查总体质量矩阵MM的主元是否为零
[sdof1,n2]=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=[kk(1:sdof,1:sdof)]; % eliminating the columns and rows of the kk
M=[mm(1:sdof,1:sdof)]; % eliminating the columns and rows of the mm
F=[ff(1:sdof,:)]; % eliminating the rows of the ff
bc_dof=[bcdof(1:sdof)]; % eliminating the rows of the bcdof
bc_val=[bcval(1:sdof)]; % eliminating the rows of the bcval
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [K,M,F]=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
|
发表于 2012-12-8 16:26
提升卡
变色卡
显身卡
楼主
