有限元三角形单元程序设计

sdof=nnode*ndof; %total system dofs
edof=nnel*ndof; %degrees of freedom per element %-------------------------------------------------------
%input data for nodal coordinate values
%--------------------------------------------
lengthx=4;
%length of x-axis side of problem
lengthy=2;
%length of y-axis side of problem
A(0, 0)
emodule=1.0;
kk=sparse(sdof,sdof);
%system matrix
disp=sparse(sdof,1);
%system displacement vector
eldisp=sparse(edof,1);
%element displacement vector
stress=zeros(nel,3);
D(4, 0) F 1
B(0, 2)
C(4, 1)
21:16
三角形单元程序
clear all
first_time=cputime;
format long
%--------------------------------------------
%input data for control parameters
%matrix containing stress components
strain=zeros(nel,3);
%matrix containing strain components
index=sparse(edof,1);
%index vector
kinmtx=sparse(3,edof);
%kinematic matrix
matmtx=sparse(3,3);
%constitutive matrix
%-------------------------------------------------
%compute material matrices
%-------------------------------------------------
%-------------------------------------------------------
x0=[];
for i=1:lx+1
for j=1:ly+1
x0=[x0; (i-1)*lengthx/lx -0.5*lengthy*(1+(lx+1-i)/lx)*(1-(j-1)/ly)];
bcdof=[bcdof 1+2*(i-1) 2+2*(i-1)]; bcval=[bcval 0 0]; end
A(0, 0) B(0, 2)
D(4, 0) C(4, 1)
21:16
三角形单元程序
%-------------------------------------------------
%initialization of matrices and vectors
%-------------------------------------------------
ff=sparse(sdof,1);
%system force vector
k=sparse(edof,edof);
%initialization of element matrix
%eln=0.0;
%Poisson's ratio
fload=-1;
% the total load
B(0, 2)
21:16
D(4, 0) F 1 C(4, 1)
三角形单元程序
lx=16;
% number of element in x-axis
ly=8;
% number of element in y-axis
end
end
21:16
三角形单元程序
%---------------------------------------------------------------------%input data for nodal connectivity for each element %nodes(i,j) where i->element no. and j->connected nodes %---------------------------------------------------------------------nodes=[]; for i=1:lx
nel=2*lx*ly;
% number of element
nnel=3;
%number of nodes per element
ndof=2;
%number of dofs per node
nnode=(lx+1)*(ly+1); %total number of nodes in system
for j=1:ly nodes=[nodes; (ly+1)*(i-1)+j (ly+1)*i+j (ly+1)*(i-1)+j+1;]; nodes=[nodes; (ly+1)*i+j (ly+1)*i+j+1 (ly+1)*(i-1)+j+1;];
end end
21:16
三角形单元程序
%---------------------------------%input data for boundary conditions %---------------------------------bcdof=[]; bcval=[]; for i=1:ly+1
三角形单元程序设计
问题描述
考虑一个平面应力问题如图所示，假设厚度h=1，材料为各项 同性，杨氏模量为E=1，泊松比为ν=0，相关力和位移边界条 件如图中所示，问题左端为固定约束。试用两个三角形单元 分析此问题，三角形单元的网格划分如图所示。试求问题各 节点位移u、v和应力σx，σy和σxy。
A(0, 0)