原创力文档- 1、本文档共15页，可阅读全部内容。
- 2、原创力文档（book118）网站文档一经付费（服务费），不意味着购买了该文档的版权，仅供个人/单位学习、研究之用，不得用于商业用途，未经授权，严禁复制、发行、汇编、翻译或者网络传播等，侵权必究。
- 3、本站所有内容均由合作方或网友上传，本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺！文档内容仅供研究参考，付费前请自行鉴别。如您付费，意味着您自己接受本站规则且自行承担风险，本站不退款、不进行额外附加服务；查看《如何避免下载的几个坑》。如果您已付费下载过本站文档，您可以点击 这里二次下载。
- 4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等，请点击“版权申诉”（推荐），也可以打举报电话：400-050-0827(电话支持时间：9:00-18:30)。
查看更多
平面四节点等参单元matlab实现[精选]
计算力学报告 平面四节点等参单元
学生姓名： 朱 彬
学号：
问题描述及分析
在无限大平面内有一个小圆孔。孔内有一集中力p，试求用有限元法编程和用ANSYS软件求出各点应力分量和位移分量，并比较二者结果。
根据圣维南原理建立半径为10mm的大圆，设小圆孔的半径a=0.5mm，在远离大圆边界的地方模型是比较精确的。由于作用在小圆孔上的力引起的位移随距离的衰减非常快，所以可以把大圆边界条件设为位移为零。
二、有限元划分描述
在划分单元时，单元数量比较多，于是我采取了使用ansys软件建模自动划分单元网格的方法。具体操作如下：
打开ansys，在单元类型中选择solid-Quad 4 node 182单元；建立类半径为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); %初始总体刚度矩阵
for i=1:m
syms Ks Et x y I1 I2 a b B; %定义可能存在的变量
e=[1,1;-1,1;-1,-1;1,-1];
for j=1:4 %形函数
N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et);
end
x=0;y=0;
for j=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;
for j=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];%高斯积分权系数
for j=1:5
文档评论（0）