平面四节点等参单元matlab实现教材.docx
《平面四节点等参单元matlab实现教材.docx》由会员分享，可在线阅读，更多相关《平面四节点等参单元matlab实现教材.docx（13页珍藏版）》请在冰豆网上搜索。
平面四节点等参单元matlab实现教材
计算力学报告
平面四节点等参单元
学生姓名：
朱彬
学号：
20100311
一、问题描述及分析
在无限大平面内有一个小圆孔。
孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。
根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。
由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。
二、有限元划分描述
在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。
具体操作如下：
打开ansys，在单元类型中选择solid->Quad4node182单元；建立类半径为0.5外半径为10的圆环；使用mashtool中的智能划分和将单元退化成三角形单元；使用工具栏中List中的Nodes和Elements选项将节点和单元数据导出并导入Excle中，总共得到了207个单元和229个节点。
如下图：
图1
三、有限元程序及求解
程序求解使用了matlab语言。
具体如下：
程序：
clc
clear
E=2e11;%弹性模量
NU=0.3;%泊松比
t=0.1;%厚度
X=xlsread（'D:
\data','nodes'）;%读取节点坐标
elem=xlsread（'D:
\data','elements'）;%读取单元编号
w=[1,2,3,4,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28];%有位移约束的节点
n=size（X,1）;%节点数
m=size（elem,1）;%单元数
K=zeros（2*n）;%初始总体刚度矩阵
fori=1:
m
symsKsEtxyI1I2abB;%定义可能存在的变量
e=[1,1;-1,1;-1,-1;1,-1];
forj=1:
4%形函数
N（j）=0.25*（1+e（j,1）*Ks）*（1+e（j,2）*Et）;
end
x=0;y=0;
forj=1:
4%标准母单元映射到真实单元
x=x+N（j）*X（elem（i,j）,1）;
y=y+N（j）*X（elem（i,j）,2）;
end
J1=jacobian（[x;y],[Ks;Et]）;%雅克比矩阵及其转置
J=J1';
forj=1:
4
I1=diff（N（j）,Ks）;%形函数分别对Ks和Et的偏导数
I2=diff（N（j）,Et）;
C=（J^-1）*[I1;I2];
a=C
（1）;%形函数对x，y的偏导数
b=C
（2）;
B（1,2*j-1）=a;%组成B阵
B（1,2*j）=0;
B（2,2*j-1）=0;
B（2,2*j）=b;
B（3,2*j-1）=b;
B（3,2*j）=a;
end
D=（E/（1-NU*NU））*[1,NU,0;NU,1,0;0,0,（1-NU）/2];%D阵
k=zeros（8,8）;
Kss=[-0.906179,-0.538469,0,0.538469,0.906179];%5*5高斯积分点
Ett=[-0.906179,-0.538469,0,0.538469,0.906179];
H=[0.236926,0.478628,0.568888,0.478628,0.236926];%高斯积分权系数
forj=1:
5%高斯积分求单元刚度阵
forl=1:
5
Ks=Kss（j）;Et=Ett（l）;
B=subs（B）;
J=subs（J）;
k=k+H（j）*H（l）*B'*D*B*det（J）;
end
end
G=zeros（8,2*n）;%初始总刚变换矩阵
G（1,2*elem（i,1）-1）=1;%总刚变换矩阵
G（2,2*elem（i,1））=1;
G（3,2*elem（i,2）-1）=1;
G（4,2*elem（i,2））=1;
G（5,2*elem（i,3）-1）=1;
G（6,2*elem（i,3））=1;
G（7,2*elem（i,4）-1）=1;
G（8,2*elem（i,4））=1;
K=K+G'*k*G;%总体刚度矩阵合成
end
KK=K;
b=size（w,1）;
fori=1:
b
K（2*w（i）-1,2*w（i）-1）=1e20;
K（2*w（i）,2*w（i））=1e20;
end%置大数法
f=zeros（2*n,1）;%初始载荷矩阵
f（10）=-10e3;%加载荷10kN
U=K\f;%节点位移
fori=1:
m%将每个单元各个节点位移集合
u（:
i）=[U（2*elem（i,1）-1）;U（2*elem（i,1））;U（2*elem（i,2）-1）;U（2*elem（i,2））;U（2*elem（i,3）-1）;U（2*elem（i,3））;U（2*elem（i,4）-1）;U（2*elem（i,4））];
end
fori=1:
m%求单元应力
symsKsEtxyI1I2abB;
e=[1,1;-1,1;-1,-1;1,-1];
forj=1:
4
N（j）=0.25*（1+e（j,1）*Ks）*（1+e（j,2）*Et）;
end
x=0;y=0;
forj=1:
4
x=x+N（j）*X（elem（i,j）,1）;
y=y+N（j）*X（elem（i,j）,2）;
end
J1=jacobian（[x;y],[Ks;Et]）;
J=J1';
forj=1:
4
I1=diff（N（j）,Ks）;
I2=diff（N（j）,Et）;
C=（J^-1）*[I1;I2];
a=C
（1）;
b=C
（2）;
B（1,2*j-1）=a;
B（1,2*j）=0;
B（2,2*j-1）=0;
B（2,2*j）=b;
B（3,2*j-1）=b;
B（3,2*j）=a;%以上同前面部分为得到B阵
end
D=（E/（1-NU*NU））*[1,NU,0;NU,1,0;0,0,（1-NU）/2];
w=D*B*u（:
i）;
w1=subs（w,{Ks,Et},{1,1}）;%求单元上各节点的应力
sigma1（:
i）=double（w1）;
w2=subs（w,{Ks,Et},{-1,1}）;
sigma2（:
i）=double（w2）;
w3=subs（w,{Ks,Et},{-1,-1}）;
sigma3（:
i）=double（w3）;
w4=subs（w,{Ks,Et},{1,-1}）;
sigma4（:
i）=double（w4）;
end
c=[24,29,47,58,78,79,137,149,160,166,186]';%如截图选取圆半径方向的节点号
d=[166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185]';%圆周方向选择的节点号
e=size（c,1）;
f=size（d,1）;
sigmar=zeros（3,e）;%r相同角度不同节点应力矩阵
sigmat=zeros（3,f）;%角度不同r不同节点应力矩阵
msigmar=zeros（1,e）;%半径方向节点处的mises应力
msigmat=zeros（1,f）;%圆周方向节点处的mises应力
fori=1:
e%绕节点平均
g=0;
forj=1:
m%判断节点在单元的那个位置并加上相应的应力值
ifelem（j,1）-c（i）==0
g=g+1;
sigmar（:
i）=sigmar（:
i）+sigma1（:
j）;
end
ifelem（j,2）-c（i）==0
g=g+1;
sigmar（:
i）=sigmar（:
i）+sigma2（:
j）;
end
ifelem（j,3）-c（i）==0
g=g+1;
sigmar（:
i）=sigmar（:
i）+sigma3（:
j）;
end
ifelem（j,4）-c（i）==0
g=g+1;
sigmar（:
i）=sigmar（:
i）+sigma4（:
j）;
end
end
sigmar（:
i）=sigmar（:
i）/g;%求应力绕节点平均
msigmar（:
i）=（0.5*（（sigmar（1,i）-sigmar（2,i））^2+sigmar（1,i）^2+sigmar（2,i）^2+6*（sigmar（3,i））^2））^0.5;%求节点处的mises应力
end
msigmar%mises应力
fori=1:
f%同上
g=0;
forj=1:
m
ifelem（j,1）-d（i）==0
g=g+1;
sigmat（:
i）=sigmat（:
i）+sigma1（:
j）;
end
ifelem（j,2）-d（i）==0
g=g+1;
sigmat（:
i）=sigmat（:
i）+sigma2（:
j）;
end
ifelem（j,3）-d（i）==0
g=g+1;
sigmat（:
i）=sigmat（:
i）+sigma3（:
j）;
end
ifelem（j,4）-d（i）==0
g=g+1;
sigmat（:
i）=sigmat（:
i）+sigma4（:
j）;
end
end
sigmat（:
i）=sigmat（:
i）/g;
msigmat（:
i）=（0.5*（（sigmat（1,i）-sigmat（2,i））^2+sigmat（1,i）^2+sigmat（2,i）^2+6*（sigmat（3,i））^2））^0.5;
end
msigmat
四、计算结果：
1.ANSYS软件计算结果
计算结果分别罗列圆周方向的单元和半径方向的单元位移和应力。
选取半径方向的单元编号为：
c=[24,29,47,58,78,79,137,149,160,166,186]';
选取圆周方向上的单元编号为：
d=[166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185]';
其在图中的位置为图1中红线标注：
图2
在处理ANSYS计算结果时，我导出了节点x，y方向的正应力和剪应力，将其导入excle中并去掉无用项如图2所示，并编写matlab程序将选取的点的mises应力求出来。
图2，求节点mises应力的程序以及计算结果如下：
图3
求选取点的mises应力程序：
clc
A=xlsread（'D:
\data','ansys'）;
c=[24,29,47,58,78,79,137,149,160,166,186]';
d=[166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185]';
e=size（c,1）
f=size（d,
升级会员 