四边形八节点等参元matlab程序

悬臂钢梁，尺寸如图一所示；v=0.3。

h=1，E=2.1e11.
图一悬臂钢梁
图二单元划分与结点编号
Matlab 输出结果
附录Ⅰ：
有限元ANSYS分析结果
采用PLANE183单元（四边形八节点）单元得出的结构Y向最大位移为-0.216E-04。

约等于matlab平面四边形八节点等参元结点Y向最大位移-2.4024E-5。

附录Ⅱ：
%---------------四边形八节点等参元 matlab计算程序---------------------------- %———————————主程序—————————
%*******************************************************************%* ***********************************
% 2012年
% 11级防灾张子凉
% 学号：2111118088
% 本程序只能处理集中荷载作用下的情况
% 只输出了节点位移、单元中心点的应力
%*******************************************************************%* **************
% 变量说明
% E v h
% 弹性模量泊松比厚度
% NPOIN NELEM NVFIX NNODE NFPOIN
% 总结点数 , 单元数, 约束结点个数, 单元节点数 ,受力结点数
% COORD LNODS
% 结构节点整体坐标数组, 单元定义数组,
% FPOIN FORCE FIXED
% 结点力数组，总体荷载向量, 约束信息数组
% HK DISP
% 总体刚度矩阵,结点位移向量
%******************************
clear all
format short e
FP1=fopen('bjd.txt','rt'); %打开数据文件
%%读入控制数据
E=fscanf(FP1,'%f',1); %弹性模量
v=fscanf(FP1,'%f',1); % 泊松比
h=fscanf(FP1,'%f',1); %厚度
NELEM=fscanf(FP1,'%d',1); %单元数
NPOIN=fscanf(FP1,'%d',1); % 总结点数
NNODE=fscanf(FP1,'%d',1); %单元节点数
NFPOIN=fscanf(FP1,'%d',1); %受力结点数
NVFIX=fscanf(FP1,'%d',1); %约束结点个数
LNODS=fscanf(FP1,'%f',[NNODE,NELEM])'; % 单元定义：单元结点号（逆时针）COORD=fscanf(FP1,'%f',[2,NPOIN])'; % 结点号 x,y坐标(整体坐标下)
FPOIN=fscanf(FP1,'%f',[3,NFPOIN])';
% 节点力：结点号、X方向力(向右正)，Y方向力(向上正)
FIXED=fscanf(FP1,'%d',[3,NVFIX])';
%约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号
%(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0
%******************************************************************* %******************************************************************* %========平面应力问题的求解==============
%
%******************************************************************* %******************************************************************* %—————————————————————
% 刚度矩阵的生成
%计算刚度矩阵，并对约束条件进行处理
Ke=zeros(2*NNODE,2*NNODE); % 单元刚度矩阵并清零
HK=zeros(2*NPOIN,2*NPOIN); % 张成总刚矩阵并清零
%调用子程序生成单元刚度矩阵
for m=1:NELEM %m为单元号
Ke=K(E,v,h,...
COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...
COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...
COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),...
COORD(LNODS(m,7),1),COORD(LNODS(m,7),2)); %调用单元刚度矩阵
a=LNODS(m,:); %临时向量,用来记录当前单元的节点编号
%对总刚度矩阵的处理
for j=1:8
for k=1:8
HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-
1):a(k)*2)+...
Ke(j*2-1:j*2,k*2-1:k*2);
end
end
end
%—————————————————————————————————%对荷载向量进行处理
FORCE=zeros(2*NPOIN,1); % 张成总荷载向量并清零
for i=1:NFPOIN
b1=FPOIN(i,1)*2-1;b2=FPOIN(i,1)*2; %FPION(i,1)为作用点
FORCE(b1)=FPOIN(i,2); %FPION(i,2)为x方向的节点力
FORCE(b2)=FPOIN(i,3); %FPION(i,3)为y方向的节点力
end %—————————————————————————————————%将约束信息加入总刚,总荷载
for i=1:NVFIX
if FIXED(i,2)==1
c1=2*FIXED(i,1)-1;
HK(c1,:)=0; %将一约束序号处的总刚列向量清0
HK(:,c1)=0; %将一约束序号处的总刚行向量清0
HK(c1,c1)=1; %将行列交叉处的元素置为1
FORCE(c1)=0;
end
if FIXED(i,3)==1
c2=2*FIXED(i,1);
HK(c2,:)=0;
HK(:,c2)=0;
HK(c2,c2)=1;
FORCE(c2)=0;
end
end %—————————————————————————————————%=========================================================== %=========================================================== DISP=HK\FORCE %计算节点位移向量
%=========================================================== %=========================================================== %———————————求解单元应力————————————————stress=zeros(3,NELEM);
for m=1:NELEM
u(1:16)=0;
d=LNODS(m,:); %临时向量,用来记录当前单元的节点编号
for i=1:NNODE
u(i*2-1:i*2)=DISP(d(i)*2-1:d(i)*2);
%从总位移向量中取出当前单元的节点位移
end
D=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%弹性矩阵
%形成应变矩阵BM
BM=zeros(3,16);
for i=1:NNODE
J=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),... COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...
COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),...
COORD(LNODS(m,7),1),COORD(LNODS(m,7),2),0,0);
[N_s,N_t]=DHS(0,0);
B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);
B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);
BM(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i]/det(J);
end
stressm=D*BM*u';
stress(:,m)=stressm;
end
stress %输出应力
function Ke=K(E,v,h,x1,y1,x3,y3,x5,y5,x7,y7)
%=========单元刚度矩阵===============
% E 弹性模量
% v 泊松比
% h 厚度
% x1,y1,x3,y3,x5,y5,x7,y7 为4个角结点的坐标
%矩阵尺寸为16 x 16
Ke=zeros(16,16);
D=(E/(1-v*v))*[1 v 0;v 1 0;0 0 (1-v)/2];%弹性矩阵
%高斯积分采用 3 x 3 个积分点书74页
W1=5/9;W2=8/9;W3=5/9; %加权系数
W=[W1 W2 W3];
r=15^(1/2)/5;
x=[-r 0 r];%积分点
for i=1:3
B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j));
J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j));
Ke=Ke+W(i)*W(j)*B'*D*B*det(J)*h;
end
end
function B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,s,t) %调用导函数
[N_s,N_t]=DHS(s,t);
%求Jacobi矩阵
J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t);
%求应变矩阵B
B=zeros(3,16);
for i=1:8
B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);
B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);
B(1:3,2*i-1:2*i)=[B1i 0;0 B2i;B2i B1i];
end
B=B/det(J);
function J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t) %-------Jacobi-----------
%单元坐标
%2,4,6,8点的坐标
x2=(x1+x3)/2;y2=(y1+y3)/2;
x4=(x3+x5)/2;y4=(y3+y5)/2;
x6=(x5+x7)/2;y6=(y5+y7)/2;
x8=(x7+x1)/2;y8=(y7+y1)/2;
x=[x1 x2 x3 x4 x5 x6 x7 x8];
y=[y1 y2 y3 y4 y5 y6 y7 y8];
%%调用形函数对局部坐标的导数
[N_s,N_t]=DHS(s,t);
%求Jacobi矩阵的行列式的值
x_t=0;y_t=0;
for i=1:8
x_s=x_s+N_s(i)*x(i);y_s=y_s+N_s(i)*y(i); x_t=x_t+N_t(i)*x(i);y_t=y_t+N_t(i)*y(i); end
J=[x_s y_s;x_t y_t];
function N=shape(s,t)
%ξ,η
N(1) = (1-s)*(1-t)*(-s-t-1)/4;
N(3) = (1+s)*(1-t)*(s-t-1)/4;
N(5) = (1+s)*(1+t)*(s+t-1)/4;
N(7) = (1-s)*(1+t)*(-s+t-1)/4;
N(2) = (1-t)*(1+s)*(1-s)/2;
N(4) = (1+s)*(1+t)*(1-t)/2;
N(6) = (1+t)*(1+s)*(1-s)/2;
N(8) = (1-s)*(1+t)*(1-t)/2;
function [N_s,N_t]=DHS(s,t)
%形函数求导
%ξ,η
N_s(1)=-1/4*(1-t)*(-s-t-1)-1/4*(1-s)*(1-t);
N_s(3)=1/4*(1-t)*(s-t-1)+1/4*(1+s)*(1-t);
N_s(5)=1/4*(1+t)*(s+t-1)+1/4*(1+s)*(1+t); N_s(7)=-1/4*(1+t)*(-s+t-1)-1/4*(1-s)*(1+t); N_s(2)=1/2*(1-s)*(1-t)-1/2*(1+s)*(1-t);
N_s(4)=1/2*(1+t)*(1-t);
N_s(6)=1/2*(1-s)*(1+t)-1/2*(1+s)*(1+t);
N_s(8)=-1/2*(1+t)*(1-t);
N_t(1)=-1/4*(1-s)*(-s-t-1)-1/4*(1-s)*(1-t);
N_t(3)=-1/4*(1+s)*(s-t-1)-1/4*(1+s)*(1-t);
N_t(5)=1/4*(1+s)*(s+t-1)+1/4*(1+s)*(1+t); N_t(7)=1/4*(1-s)*(-s+t-1)+1/4*(1-s)*(1+t); N_t(2)=-1/2*(1+s)*(1-s);
N_t(4)=1/2*(1+s)*(1-t)-1/2*(1+s)*(1+t);
N_t(6)=1/2*(1+s)*(1-s);
N_t(8)=1/2*(1-s)*(1-t)-1/2*(1-s)*(1+t);
bjd.txt 文件数据
2.1E11 0.3 1 5 28 8 1 3 1 2 3 13 20 19 18 12
3 4 5 14 22 21 20 13
5 6 7 15 24 23 22 14
7 8 9 16 26 25 24 15
9 10 11 17 28 27 26 16 0.0 0.0
0.5 0.0
1.0 0.0
1.5 0.0
2.0 0.0
2.5 0.0
3.0 0.0
3.5 0.0
4.0 0.0
4.5 0.0
5.0 0.0
0.0 0.5
1.0 0.5
2.0 0.5
3.0 0.5
4.0 0.5
5.0 0.5
0.0 1.0
0.5 1.0
1.0 1.0
1.5 1.0
2.0 1.0
2.5 1.0
3.0 1.0
3.5 1.0
4.0 1.0
4.5 1.0
5.0 1.0
17 0 -10000
1 1 1
12 1 1
18 1 1
11 / 11。