FEM-Course-Matlab/4.mindlin厚板理论及matlab有限元编程/main_mindlin.m

43 lines
1.4 KiB
Mathematica
Raw Permalink Normal View History

2024-01-28 16:46:36 +00:00
clear all;close all;clc
%---------------------GEOMETRY-----------------------
a=100; %('<EFBFBD><EFBFBD><EFBFBD><EFBFBD> in x: ');
b=100; %('<EFBFBD><EFBFBD><EFBFBD><EFBFBD> in y: ');
%---------------------MESH---------------------------
ex=20; %input('x<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>: ');
ey=20; %input('y<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>: ');
mesh=MeshGenerator(a,b,ex,ey);
%---------------------ELEMENT------------------------
ShapeOption='Q4';%Q8 Q9
quadrature=GaussQuadrature('gauss2');% Q4 element points and weights
%---------------------MATERIAL------------------------
%SI mm t MPa N
material=struct();
material.t=1.6;
material.E=72000;
material.v=0.33;
material.rho=2.7e-9;
material.G=material.E/(2+2*material.v);
K=StiffnessMatrix(material,mesh,quadrature,ShapeOption);
M=MassMatrix(material,mesh,quadrature,ShapeOption);
%---------------------LOADS----------------------------
p0=1e-3;
F=ForceVector(mesh,quadrature,ShapeOption,p0);
% c0=input('\n<EFBFBD><EFBFBD><EFBFBD>к<EFBFBD><EFBFBD>صĴ<EFBFBD>С: ');
% cl=input('\nLato da caricare: ');
% F=ForceVector(cl,c0,mesh);Lato da caricare
%---------------------CONSTRAINTS-----------------------
%nc is nodes of constrained elemment
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);
%---------------------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')