四边形八节点等参元matlab程序讲解学习文档格式.docx
- 文档编号:19563141
- 上传时间:2023-01-07
- 格式:DOCX
- 页数:12
- 大小:84.60KB
四边形八节点等参元matlab程序讲解学习文档格式.docx
《四边形八节点等参元matlab程序讲解学习文档格式.docx》由会员分享,可在线阅读,更多相关《四边形八节点等参元matlab程序讲解学习文档格式.docx(12页珍藏版)》请在冰豆网上搜索。
%结点力数组,总体荷载向量,约束信息数组
%HKDISP
%总体刚度矩阵,结点位移向量
%******************************
clearall
formatshorte
FP1=fopen('
bjd.txt'
'
rt'
);
%打开数据文件
%%读入控制数据
E=fscanf(FP1,'
%f'
1);
%弹性模量
v=fscanf(FP1,'
%泊松比
h=fscanf(FP1,'
%厚度
NELEM=fscanf(FP1,'
%d'
%单元数
NPOIN=fscanf(FP1,'
%总结点数
NNODE=fscanf(FP1,'
%单元节点数
NFPOIN=fscanf(FP1,'
%受力结点数
NVFIX=fscanf(FP1,'
%约束结点个数
LNODS=fscanf(FP1,'
[NNODE,NELEM])'
;
%单元定义:
单元结点号(逆时针)
COORD=fscanf(FP1,'
[2,NPOIN])'
%结点号x,y坐标(整体坐标下)
FPOIN=fscanf(FP1,'
[3,NFPOIN])'
%节点力:
结点号、X方向力(向右正),Y方向力(向上正)
FIXED=fscanf(FP1,'
[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);
%张成总刚矩阵并清零
%调用子程序生成单元刚度矩阵
form=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,:
%临时向量,用来记录当前单元的节点编号
%对总刚度矩阵的处理
forj=1:
8
fork=1:
HK((a(j)*2-1):
a(j)*2,(a(k)*2-1):
a(k)*2)=HK((a(j)*2-1):
a(k)*2)+...
Ke(j*2-1:
j*2,k*2-1:
k*2);
end
end
%—————————————————————————————————
%对荷载向量进行处理
FORCE=zeros(2*NPOIN,1);
%张成总荷载向量并清零
fori=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方向的节点力
%将约束信息加入总刚,总荷载
fori=1:
NVFIX
ifFIXED(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;
ifFIXED(i,3)==1
c2=2*FIXED(i,1);
HK(c2,:
c2)=0;
HK(c2,c2)=1;
FORCE(c2)=0;
%===========================================================
DISP=HK\FORCE%计算节点位移向量
%———————————求解单元应力————————————————
stress=zeros(3,NELEM);
NELEM
u(1:
16)=0;
d=LNODS(m,:
NNODE
u(i*2-1:
i*2)=DISP(d(i)*2-1:
d(i)*2);
%从总位移向量中取出当前单元的节点位移
D=(E/(1-v*v))*[1v0;
v10;
00(1-v)/2];
%弹性矩阵
%形成应变矩阵BM
BM=zeros(3,16);
J=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),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)=[B1i0;
0B2i;
B2iB1i]/det(J);
stressm=D*BM*u'
stress(:
m)=stressm;
stress%输出应力
functionKe=K(E,v,h,x1,y1,x3,y3,x5,y5,x7,y7)
%=========单元刚度矩阵===============
%E弹性模量
%v泊松比
%h厚度
%x1,y1,x3,y3,x5,y5,x7,y7为4个角结点的坐标
%矩阵尺寸为16x16
Ke=zeros(16,16);
D=(E/(1-v*v))*[1v0;
%高斯积分采用3x3个积分点书74页
W1=5/9;
W2=8/9;
W3=5/9;
%加权系数
W=[W1W2W3];
r=15^(1/2)/5;
x=[-r0r];
%积分点
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;
functionB=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);
B(1:
B2iB1i];
B=B/det(J);
functionJ=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=[x1x2x3x4x5x6x7x8];
y=[y1y2y3y4y5y6y7y8];
%%调用形函数对局部坐标的导数
%求Jacobi矩阵的行列式的值
x_s=0;
y_s=0;
x_t=0;
y_t=0;
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);
J=[x_sy_s;
x_ty_t];
functionN=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.1E110.31528813
1231320191812
3451422212013
5671524232214
7891626252415
910111728272616
0.00.0
0.50.0
1.00.0
1.50.0
2.00.0
2.50.0
3.00.0
3.50.0
4.00.0
4.50.0
5.00.0
0.00.5
1.00.5
2.00.5
3.00.5
4.00.5
5.00.5
0.01.0
0.51.0
1.01.0
1.51.0
2.01.0
2.51.0
3.01.0
3.51.0
4.01.0
4.51.0
5.01.0
170-10000
111
1211
1811
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 四边形 八节 matlab 程序 讲解 学习