79 lines
2.5 KiB
Matlab
79 lines
2.5 KiB
Matlab
clear all;close all;clc
|
||
%---------------------GEOMETRY-----------------------
|
||
a=10; %('长度 in x: ');
|
||
b=10; %('宽度 in y: ');
|
||
%---------------------MESH---------------------------
|
||
ex=15; %input('x向单元个数: ');%如果采用完全积分2*2,对于薄板，受剪切自锁现象影响，网格大小对结果有很大的影响
|
||
ey=15; %input('y向单元个数: ');
|
||
mesh=MeshGenerator(a,b,ex,ey);
|
||
%---------------------ELEMENT------------------------
|
||
ShapeOption='Q4';%Q8 Q9
|
||
quadrature=GaussQuadrature('gauss1');
|
||
% gauss1 为减缩积分，采用一个积分点避免剪切自锁，尤其对于薄板
|
||
%---------------------MATERIAL------------------------
|
||
%SI mm t MPa N
|
||
material=struct();
|
||
material.t=0.5; %注意薄板存在剪切自锁现象，可通过尝试用减缩积分（采用1个积分点）来避免
|
||
material.E=200e9;
|
||
material.v=0.3;
|
||
material.rho=7800;
|
||
material.G=material.E/(2+2*material.v);
|
||
|
||
K=StiffnessMatrix(material,mesh,quadrature,ShapeOption);
|
||
M=MassMatrix(material,mesh,quadrature,ShapeOption);
|
||
%---------------------LOADS----------------------------
|
||
p0=1e3;
|
||
F=ForceVector(mesh,quadrature,ShapeOption,p0);
|
||
% c0=input('\n集中荷载的大小: ');
|
||
% cl=input('\nLato da caricare: ');
|
||
% F=ForceVector(cl,c0,mesh);Lato da caricare
|
||
%---------------------CONSTRAINTS-----------------------
|
||
%nc is nodes of constrained elemment
|
||
nc=[mesh.lato1];%悬臂板
|
||
%nc=[mesh.lato1 mesh.lato2(2:end) mesh.lato3(2:end) mesh.lato4(2:end-1)];%%四边固定
|
||
[K_c,M_c,F_c,nctot]=Constraints(nc,K,M,F);
|
||
%---------------------SOLUTION---------------------------
|
||
w=StaticSolver(K_c,F_c,mesh,nctot);
|
||
sprintf('最大位移为：%f',max(abs(w)))
|
||
%---------------------PLOT-------------------------------
|
||
figure
|
||
surf(mesh.x(1,:),mesh.y(:,1),reshape(w,size(mesh.x,2),size(mesh.x,1))');
|
||
% axis equal
|
||
hold on
|
||
title('Deformed geometry')
|
||
%-------------------MODAL ANALYSIS--------------------
|
||
|
||
numberMode=5;
|
||
|
||
% %原始模态计算方法'sm'
|
||
% [V,D]=eigs(M_c\K_c,numberMode,'sm');%'smallestabs'
|
||
% scal=0.5;
|
||
% f=diag(sqrt(D)/(2*pi));
|
||
% Vf=sortrows([V',f],size(V,1)+1);
|
||
% V=Vf(:,1:size(V,1))';%排序后的特征向量
|
||
|
||
%如果M奇异,反转K与M，求K得逆，并删除‘sm'用'largestabs'替代
|
||
[V,D]=eigs(K_c\M_c,numberMode,'lm'); %'largestabs'
|
||
scal=0.5;
|
||
f= power(diag(sqrt(D)),-1)/(2*pi);
|
||
Vf=sortrows([V',f],size(V,1)+1);
|
||
V=Vf(:,1:size(V,1))';%排序后的特征向量
|
||
|
||
|
||
%加入约束自由度
|
||
index=1:mesh.nn*3;
|
||
Vo=zeros(mesh.nn*3,5);
|
||
index(nctot)=[];%删除被约束的自由度
|
||
for i=1:numel(index)
|
||
Vo(index(i),:)=V(i,:);
|
||
end
|
||
|
||
f=sort(f);%排序后的频率
|
||
sprintf('频率为%f',f)
|
||
|
||
for i=1:numberMode
|
||
figure
|
||
surf(mesh.x(1,:),mesh.y(:,1),reshape(Vo(1:3:end,i),size(mesh.x,2),size(mesh.x,1))');
|
||
title(sprintf('modal %d',i))
|
||
end
|