马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
用matlab编程计算悬臂薄板前十五阶固有频率 ，结果为;
frequency =
1.0e+004 *
Columns 1 through 8
0 + 0.0000i 0.0000 0.0000 0.1837 0.5090 0.5852 1.0043 1.1981
Columns 9 through 15
1.6716 1.8640 2.5144 2.6052 3.4420 3.5358 4.3945
是否因为引入边界条件有问题，应如何更改？
程序如下：
function gy
PlaneFrameModel ; % 定义有限元模型
SolveModel ; % 求解有限元模型
return ;
function PlaneFrameModel
global lx ly jdx jdy gNode gElement gMaterial gBC1 ke me gK gM
% 给定几何特征
E=2.1e11; %elastic molulus
poisson =0.3; % poisson ratio
density=7.85e3; %density
t=0.000152; %plate thickness
lx=0.021; %length in x direction
ly=0.004; %length in y direction
jdx=11; %number of nodes in x direction
jdy=11; %number of nodes in y direction
% 计算结点坐标
dx = lx / (jdx-1) ;
dy = ly / (jdy-1) ;
gNode = zeros( jdx*jdy, 2 ) ;
for i=1:jdx
for j=1:jdy
gNode( (i-1)*jdy+j, : ) = [dx*(j-1),dy*(i-1)] ;
end
end
% 定义单元
gElement = zeros( (jdx-1)*(jdy-1), 5 ) ;
for i=1:(jdx-1)
for j=1:(jdy-1)
gElement( (i-1)*(jdx-1)+j, 1:4) = [ (i-1)*jdx+j, ...
(i-1)*jdx+j+1, ...
i*jdx+j+1,...
i*jdx+j ] ;
end
end
gElement( :, 5 ) = 1 ;
% 定义材料
gMaterial = [ E, poisson, t, density] ;
% 确定边界条件
gBC1 = zeros( jdx*3, 3 ) ;
for i=1:jdx
gBC1( i, : ) = [1+(i-1)*jdy, 1, 0 ] ; % x=0的边界上挠度等于零
end
for i=1:jdx
gBC1( jdx+i, : ) = [1+(i-1)*jdy, 2, 0 ] ; % x=0的边界上绕x轴的转角等于零
end
for i=1:jdx
gBC1( jdx*2+i, : ) = [1+(i-1)*jdy, 3, 0 ] ; % x=0的边界上绕y轴的转角等于零
end
% 定义整体刚度矩阵和节点力向量
[node_number,dummy] = size( gNode ) ;
gK = zeros( node_number * 3, node_number * 3 ) ;
gM = zeros( node_number * 3, node_number * 3 ) ;
f = zeros( node_number * 3, 1) ;
% 计算单元刚度和质量矩阵，并集成到整体刚度和质量矩阵中
[element_number,dummy] = size( gElement ) ;
for ie=1:1:element_number
ke = StiffnessMatrix( ie ) ;
me = MassMatrix( ie ) ;
AssembleGlobalMatrix( ie, ke, me ) ;
end
return
% 对gK进行边界条件处理
[bc1_number,dummy] = size( gBC1 ) ;
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
K1im(:,ibc) = gK(:,m) ;
gK(:,m) = zeros( node_number*3, 1 ) ;
gK(m,:) = zeros( 1, node_number*3 ) ;
gK(m,m) = 1.0 ;
end
function SolveModel
global lx ly jdx jdy gNode gElement gMaterial gBC1 ke me gK gM
%solve eigenvalue problem
[v,d] = eig(gK,gM);
tempd=diag(d);
[nd,sortindex]=sort(tempd);
v=v(:,sortindex);
mode_number=1:15;
frequency(mode_number)=sqrt(nd(mode_number))/(2*pi);
frequency
return
function ke = StiffnessMatrix( ie )
% 计算单元刚度矩阵
% 输入参数:
% ie ------- 单元号
% 返回值:
% k ---- 整体坐标系下的刚度矩阵
global lx ly jdx jdy gElement gMaterial
ke = zeros( 12, 12 ) ;
E = gMaterial( gElement(ie, 5), 1 ) ;
poisson = gMaterial( gElement(ie, 5), 2 ) ;
t = gMaterial( gElement(ie, 5), 3 ) ;
density = gMaterial( gElement(ie, 5), 4 ) ;
a=lx/(jdx-1)/2; %element length
b=ly/(jdy-1)/2; %element height
ke=[E*t^3/(360-360*poisson^2)/a/b*(21-6*poisson+30*b^2/a^2+30*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(3*b+12*poisson*b+30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-3*a-12*poisson*a-30*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(-21+6*poisson-30*b^2/a^2+15*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*b-12*poisson*b+15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-3*a+3*poisson*a-30*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(21-6*poisson-15*b^2/a^2-15*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*b+3*poisson*b+15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(3*a-3*poisson*a-15*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(-21+6*poisson+15*b^2/a^2-30*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(3*b-3*poisson*b+30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(3*a+12*poisson*a-15*b^2/a);
E*t^3/(360-360*poisson^2)/a/b*(3*b+12*poisson*b+30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(8*b^2-8*poisson*b^2+40*a^2), -30*E*t^3/(360-360*poisson^2)*poisson, E*t^3/(360-360*poisson^2)/a/b*(-3*b-12*poisson*b+15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-8*b^2+8*poisson*b^2+20*a^2), 0, E*t^3/(360-360*poisson^2)/a/b*(3*b-3*poisson*b-15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(2*b^2-2*poisson*b^2+10*a^2), 0, E*t^3/(360-360*poisson^2)/a/b*(-3*b+3*poisson*b-30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-2*b^2+2*poisson*b^2+20*a^2), 0;
E*t^3/(360-360*poisson^2)/a/b*(-3*a-12*poisson*a-30*b^2/a), -30*E*t^3/(360-360*poisson^2)*poisson, E*t^3/(360-360*poisson^2)/a/b*(8*a^2-8*poisson*a^2+40*b^2), E*t^3/(360-360*poisson^2)/a/b*(3*a-3*poisson*a+30*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(-2*a^2+2*poisson*a^2+20*b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*a+3*poisson*a+15*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(2*a^2-2*poisson*a^2+10*b^2), E*t^3/(360-360*poisson^2)/a/b*(3*a+12*poisson*a-15*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(-8*a^2+8*poisson*a^2+20*b^2);
E*t^3/(360-360*poisson^2)/a/b*(-21+6*poisson-30*b^2/a^2+15*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*b-12*poisson*b+15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(3*a-3*poisson*a+30*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(21-6*poisson+30*b^2/a^2+30*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(3*b+12*poisson*b+30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(3*a+12*poisson*a+30*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(-21+6*poisson+15*b^2/a^2-30*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(3*b-3*poisson*b+30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-3*a-12*poisson*a+15*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(21-6*poisson-15*b^2/a^2-15*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*b+3*poisson*b+15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-3*a+3*poisson*a+15*b^2/a);
E*t^3/(360-360*poisson^2)/a/b*(-3*b-12*poisson*b+15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-8*b^2+8*poisson*b^2+20*a^2), 0, E*t^3/(360-360*poisson^2)/a/b*(3*b+12*poisson*b+30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(8*b^2-8*poisson*b^2+40*a^2), 30*E*t^3/(360-360*poisson^2)*poisson, E*t^3/(360-360*poisson^2)/a/b*(-3*b+3*poisson*b-30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-2*b^2+2*poisson*b^2+20*a^2), 0, E*t^3/(360-360*poisson^2)/a/b*(3*b-3*poisson*b-15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(2*b^2-2*poisson*b^2+10*a^2), 0;
E*t^3/(360-360*poisson^2)/a/b*(-3*a+3*poisson*a-30*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(-2*a^2+2*poisson*a^2+20*b^2), E*t^3/(360-360*poisson^2)/a/b*(3*a+12*poisson*a+30*b^2/a), 30*E*t^3/(360-360*poisson^2)*poisson, E*t^3/(360-360*poisson^2)/a/b*(8*a^2-8*poisson*a^2+40*b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*a-12*poisson*a+15*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(-8*a^2+8*poisson*a^2+20*b^2), E*t^3/(360-360*poisson^2)/a/b*(3*a-3*poisson*a-15*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(2*a^2-2*poisson*a^2+10*b^2);
E*t^3/(360-360*poisson^2)/a/b*(21-6*poisson-15*b^2/a^2-15*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(3*b-3*poisson*b-15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-3*a+3*poisson*a+15*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(-21+6*poisson+15*b^2/a^2-30*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*b+3*poisson*b-30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-3*a-12*poisson*a+15*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(21-6*poisson+30*b^2/a^2+30*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*b-12*poisson*b-30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(3*a+12*poisson*a+30*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(-21+6*poisson-30*b^2/a^2+15*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(3*b+12*poisson*b-15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(3*a-3*poisson*a+30*b^2/a);
E*t^3/(360-360*poisson^2)/a/b*(-3*b+3*poisson*b+15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(2*b^2-2*poisson*b^2+10*a^2), 0, E*t^3/(360-360*poisson^2)/a/b*(3*b-3*poisson*b+30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-2*b^2+2*poisson*b^2+20*a^2), 0, E*t^3/(360-360*poisson^2)/a/b*(-3*b-12*poisson*b-30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(8*b^2-8*poisson*b^2+40*a^2), -30*E*t^3/(360-360*poisson^2)*poisson, E*t^3/(360-360*poisson^2)/a/b*(3*b+12*poisson*b-15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-8*b^2+8*poisson*b^2+20*a^2), 0;
E*t^3/(360-360*poisson^2)/a/b*(3*a-3*poisson*a-15*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(2*a^2-2*poisson*a^2+10*b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*a-12*poisson*a+15*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(-8*a^2+8*poisson*a^2+20*b^2), E*t^3/(360-360*poisson^2)/a/b*(3*a+12*poisson*a+30*b^2/a), -30*E*t^3/(360-360*poisson^2)*poisson, E*t^3/(360-360*poisson^2)/a/b*(8*a^2-8*poisson*a^2+40*b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*a+3*poisson*a-30*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(-2*a^2+2*poisson*a^2+20*b^2);
E*t^3/(360-360*poisson^2)/a/b*(-21+6*poisson+15*b^2/a^2-30*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*b+3*poisson*b-30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(3*a+12*poisson*a-15*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(21-6*poisson-15*b^2/a^2-15*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(3*b-3*poisson*b-15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(3*a-3*poisson*a-15*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(-21+6*poisson-30*b^2/a^2+15*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(3*b+12*poisson*b-15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-3*a+3*poisson*a-30*b^2/a), E*t^3/(360-360*poisson^2)/a/b*(21-6*poisson+30*b^2/a^2+30*a^2/b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*b-12*poisson*b-30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-3*a-12*poisson*a-30*b^2/a);
E*t^3/(360-360*poisson^2)/a/b*(3*b-3*poisson*b+30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-2*b^2+2*poisson*b^2+20*a^2), 0, E*t^3/(360-360*poisson^2)/a/b*(-3*b+3*poisson*b+15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(2*b^2-2*poisson*b^2+10*a^2), 0, E*t^3/(360-360*poisson^2)/a/b*(3*b+12*poisson*b-15*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(-8*b^2+8*poisson*b^2+20*a^2), 0, E*t^3/(360-360*poisson^2)/a/b*(-3*b-12*poisson*b-30*a^2/b), E*t^3/(360-360*poisson^2)/a/b*(8*b^2-8*poisson*b^2+40*a^2), 30*E*t^3/(360-360*poisson^2)*poisson;
E*t^3/(360-360*poisson^2)/a/b*(3*a+12*poisson*a-15*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(-8*a^2+8*poisson*a^2+20*b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*a+3*poisson*a+15*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(2*a^2-2*poisson*a^2+10*b^2), E*t^3/(360-360*poisson^2)/a/b*(3*a-3*poisson*a+30*b^2/a), 0, E*t^3/(360-360*poisson^2)/a/b*(-2*a^2+2*poisson*a^2+20*b^2), E*t^3/(360-360*poisson^2)/a/b*(-3*a-12*poisson*a-30*b^2/a), 30*E*t^3/(360-360*poisson^2)*poisson, E*t^3/(360-360*poisson^2)/a/b*(8*a^2-8*poisson*a^2+40*b^2)];
return
function me = MassMatrix( ie )
% 计算单元质量矩阵
% 输入参数:
% ie ------- 单元号
% 返回值:
% m ---- 整体坐标系下的质量矩阵
global lx ly jdx jdy gElement gMaterial
me = zeros( 12, 12 ) ;
E = gMaterial( gElement(ie, 5), 1 ) ;
poisson = gMaterial( gElement(ie, 5), 2 ) ;
t = gMaterial( gElement(ie, 5), 3 ) ;
density = gMaterial( gElement(ie, 5), 4 ) ;
a=lx/(jdx-1)/2; %element length
b=ly/(jdy-1)/2; %element height
w=a*b*t*density;
syms kx yt kxi yti real;
ni=1/8*(1+kx*kxi)*(1+yt*yti)*(2+kx*kxi+yt*yti-kx^2-yt^2);
nix=-1/8*b*yti*(1+kx*kxi)*(1+yt*yti)*(1-yt^2);
niy=1/8*a*kxi*(1+kx*kxi)*(1+yt*yti)*(1-kx^2);
n(1)=subs(ni,{kxi,yti},{-1,-1});
n(2)=subs(nix,{kxi,yti},{-1,-1});
n(3)=subs(niy,{kxi,yti},{-1,-1});
n(4)=subs(ni,{kxi,yti},{1,-1});
n(5)=subs(nix,{kxi,yti},{1,-1});
n(6)=subs(niy,{kxi,yti},{1,-1});
n(7)=subs(ni,{kxi,yti},{1,1});
n(8)=subs(nix,{kxi,yti},{1,1});
n(9)=subs(niy,{kxi,yti},{1,1});
n(10)=subs(ni,{kxi,yti},{-1,1});
n(11)=subs(nix,{kxi,yti},{-1,1});
n(12)=subs(niy,{kxi,yti},{-1,1});
temp=n'*n;
m1=int(temp,kx,-1,1);
me=int(m1,yt,-1,1);
me=me*w;
me=double(me);
return
function AssembleGlobalMatrix( ie, ke, me )
% 把单元刚度和质量矩阵集成到整体刚度矩阵
% 输入参数:
% ie --- 单元号
% ke --- 单元刚度矩阵
% me --- 单元质量矩阵
% 返回值:
% 无
global gElement gK gM
for i=1:1:4
for j=1:1:4
for p=1:1:3
for q =1:1:3
m = (i-1)*3+p ;
n = (j-1)*3+q ;
M = (gElement(ie,i)-1)*3+p ;
N = (gElement(ie,j)-1)*3+q ;
gK(M,N) = gK(M,N) + ke(m,n) ;
gM(M,N) = gM(M,N) + me(m,n) ;
end
end
end
end
return
|
发表于 2013-1-25 16:21
提升卡
变色卡
显身卡