有限元程序的设计大作业.docx
《有限元程序的设计大作业.docx》由会员分享，可在线阅读，更多相关《有限元程序的设计大作业.docx（48页珍藏版）》请在冰豆网上搜索。
有限元程序的设计大作业
1.不同板宽的孔边应力集中问题
姓名：
胡宇
学号：
2120120128
2.摘要
本文采用MATLAB和FOTRAN四节点平面单元，利用有限元数值解法对不同板宽的孔边应力集中问题进行了数值模拟研究。
对于不同的板宽系数（半板宽b/孔半径r），得到了不同的应力集中系数
，并且与解析解进行了比较，验证了有限元解的正确性，并且得出了解析解的适用范围。
3.引言
通常情况下的有限元分析过程是运用可视化分析软件（如ANSYS、ABAQUS、SAP等）进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。
具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。
随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中。
MATLAB是由美国MATHWORKS公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。
它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。
4.MATLAB部分
1，计算模型
本程序采用MATLAB编程，编制平面四边形四节点等参元程序，用以求解近似平面结构问题。
本程序的研究对象为中央开有小孔的长方形板，选取的材料参数为：
板厚h=1、材料强度E=1.0e11Pa、泊松比mu=0.3。
此外，为方便网格的划分和计算，本文所取板的长度与宽度相等。
其孔半径为r=1，板宽为2b待定。
由于本程序的目的在于验证有限元解的正确性和确定解析解的适用范围，因此要求网格足够细密，以满足程序的精度要求。
同时为了减小计算量，我采取网格径向长度递增的网格划分方法。
此种方法特点是，靠近小孔部分的网格细密，在远离小孔的过程中，网格逐渐变得稀疏。
以下为网格节点坐标计算和单元节点排序的MATLAB源程序：
dr=（b-r）/（m+m^2/2+m^3/6）;
dfi=（pi/2）/n;
gNode=zeros（（n+1）*（m+1）,2）;
fori=1:
1:
m
forj=1:
1:
n+1
gNode（（i-1）*（n+1）+j,:
）=[cos（dfi*（j-1））*（r+...dr*（（i-1）+（i-1）^2/2+（i-1）^3/6））,sin（dfi*（j-1））*（r+dr*（（i-1）+（i-1）^2/2+（i-1）^3/6））];
end
end
fori=1:
1:
（n/2+1）
gNode（（n+1）*m+i,:
）=[b,b*tan（dfi*（i-1））];
end
fori=1:
1:
n/2
gNode（（n+1）*m+（n/2+1）+i,:
）=[b/tan（dfi*（i+n/2））,b];
end
gElement=zeros（m*n,5）;
fori=1:
m
forj=1:
n
gElement（（i-1）*n+j,1:
4）=[（i-1）*（n+1）+j,...
i*（n+1）+j,...
i*（n+1）+j+1,...
（i-1）*（n+1）+j+1];
End
end
gElement（:
5）=1;
以上源程序中gNode为节点坐标矩阵，gElement为单元节点排列矩阵。
图
（1）为5x10的网格划分图，图
（2）为结构的微变形图。
图
（1）
图
（2）
假定该板所受的外载为p=1.0e9Pa。
由于是取1/4板进行计算，需要在分界面上加上合适的边界条件。
以小孔圆心为原点，板长方向为x轴，板宽方向为y轴。
我所取得边界条件为：
与x轴平行的分界面上的节点的y向固定，x向可以移动；与y轴平行的分界面上的节点的x向固定，y向可以移动。
划分好网格，约束和载荷加好后。
就可以计算了。
取b=10，m=5，n=10时的应力图。
图（3）为x向应力分布图，图（4）为y向应力分布图，图（5）为剪切力分布图。
图（3）
图（4）
图（5）
2，孔边应力研究
无限大板宽的孔边应力集中问题的弹性力学的解析解为：
在孔边的y轴上
有分布:
为了验证有限元解的正确性，我取10x20的网格进行了计算，得到了y轴上的x向应力分布曲线，并且与解析解进行了比较。
如图（6）。
y轴上的x向应力解析解为：
图（6）
同时，小孔边缘的剪切力分布如图（7）所示。
剪切力的解析解为
图（7）
图中，红色*点为有限元节点解，黑色x点为解析解点。
由图可知，有限元解和解析解基本相同，有限元解是有效的。
同时，为了研究不同板宽对应力集中系数的影响，我选取了多组数据进行对比。
如表
（1）所示。
表
（1）
bor（
）
b-r
2
1
1
6.0286
4
6.1114
3
2
1
3.6264
6
4.0315
4
3
2
3.2289
7
3.5248
5
4
2
2.9094
8
3.3169
6
5
2
2.7122
9
3.2098
7
6
2
2.5751
10
3.1470
8
7
3
2.8527
11
3.1069
9
8
3
2.7762
11
3.0170
10
9
3
2.7102
12
3.0604
17
16
4
2.6707
16
3.0071
26
25
5
2.6589
20
2.9911
37
36
6
2.6545
24
2.9851
50
49
7
2.6526
28
2.9824
表中，
为径向单元的数目。
分别为与之对应的应力集中系数（
）。
为半板宽b与孔半径r的比值，即
。
解析解中
的取值为3.下图（7）为解析解和有限元解的比较。
图（7）
图中黑色o点为理论应力集中系数，红色x点为较低网格密度下的应力集中系数变化曲线，蓝色*点为较高网格密度下的应力集中系数变化曲线。
由图可知，同为有限元解，网格密度越高，精确度就越高。
同时，我们还可以知道，解析解并不是适用于所有的情况。
当板宽与孔径的比值小于5时，解析解与有限元解的差别变大，即解析解不再适用。
5.FORTRAN部分
为了进一步验证结果的正确性，我同时运用FORTRAN对相同的网格进行了计算。
计算结果表明，有限元法所得的结果是可信的。
6.结论
有限元解和解析解各有优劣，有限元解的优点是基本上对任意情况都适用，但是一般情况下计算量都很大，结果获取速度慢。
解析解则求解迅速，但并不适用于所有情况。
本文所研究的孔边应力集中问题就很好的体现了两种方法的特点。
解析解并不适用于半板宽和孔半径之比小于5的情况，但是对于板宽远大于孔径的情况具有很高的精确度。
而有限元解的精确度依赖于网格的密度和算法，在算法相同的情况下，为增加精度而加密网格会使计算量增加，在精度要求不高的情况下适用。
附录一MATLAB源程序
functionhuyu2120120128
globalrbmnpEmu
%定义板的尺寸
r=1;%小孔半径
b=10;%半板宽
formatshorte
%定义网格密度
m=1;%径向划分的数目
n=2;%周向划分的数目（必须为偶数）n=2*m为最优
%定义材料性质
E=1e11;%弹性模量
mu=0.3;%泊松比
%定义均布荷载大小
p=1.0e9;%单位Pa
%
FemModel;%定义有限元模型
SolveModel;%求解有限元模型
DisplayResults;%把计算结果显示出来
return
functionFemModel
%定义有限元模型
%全局变量及名称
%r-----圆孔半径
%b-----板宽
%m---------半径方向单元数目
%n---------圆周方向单元数目
%E---------弹性模量
%mu---------泊松比
%p----------均布载荷大小
%gNode------节点坐标
%gElement---单元定义
%gMaterial--材料性质
%gBC1-------第一类约束条件
%gBC3-------斜约束
%gDF--------分布力
%返回值：
%无
globalgNodegElementgMaterialgBC1gDF
globalrbmnpEmu
%计算结点坐标
dr=（b-r）/（m+m^2/2+m^3/6）;
dfi=（pi/2）/n;
gNode=zeros（（n+1）*（m+1）,2）;
fori=1:
1:
m
forj=1:
1:
n+1
gNode（（i-1）*（n+1）+j,:
）=[cos（dfi*（j-1））*（r+...
dr*（（i-1）+（i-1）^2/2+（i-1）^3/6））,sin（dfi*（j-1））*（r+dr*（（i-1）+（i-1）^2/2+（i-1）^3/6））];
end
end
fori=1:
1:
（n/2+1）
gNode（（n+1）*m+i,:
）=[b,b*tan（dfi*（i-1））];
end
fori=1:
1:
n/2
gNode（（n+1）*m+（n/2+1）+i,:
）=[b/tan（dfi*（i+n/2））,b];
end
%单元定义
gElement=zeros（m*n,5）;
fori=1:
m
forj=1:
n
gElement（（i-1）*n+j,1:
4）=[i*（n+1）+j,...
i*（n+1）+j+1,...
（i-1）*（n+1）+j+1,...
（i-1）*（n+1）+j];
end
end
gElement（:
5）=1;
%定义材料性质
gMaterial=[E,mu,0];
%定义约束条件
gBC1=zeros（（m+1）*2,3）;
fori=1:
1:
m+1
gBC1（i,1:
2）=[（n+1）*（i-1）+1,（（n+1）*（i-1）+1）*2];
end
fori=1:
1:
m+1
gBC1（m+1+i,1:
2）=[（n+1）*i,（（n+1）*i-1）*2+1];
end
gBC1（:
3）=0;
%定义分布压力
gDF=zeros（n/2,4）;
fori=1:
1:
n/2
gDF（i,1:
2）=[（n+1）*m+i,（n+1）*m+i+1];
end
gDF（:
3）=p;
gDF（:
4）=p;
return
functionSolveModel
%求解有限元模型
%说明：
%该函数求解有限元模型，过程如下
%1.计算单元刚度矩阵，集成整体刚度矩阵
%2.计算单元的等效节点力，集成整体节点力向量
%3.处理约束条件，修改整体刚度矩阵和节点力向量
%4.求解方程组，得到整体节点位移向量gMaterial
globalgNodegElementgBC1gDFgKgDeltagElementStress
%step1.定义整体刚度矩阵和节点力向量
[node_number,dummy]=size（gNode）;
gK=sparse（node_number*2,node_number*2）;
f=sparse（node_number*2,1）;
%step2.计算单元刚度矩阵，并集成到整体刚度矩阵中
[element_number,dummy]=size（gElement）;
forie=1:
1:
element_number
fprintf（sprintf（'计算单元刚度矩阵并集成，当前单元:
%d\n',ie））;
k=StiffnessMatrix（ie）;
AssembleStiffnessMatrix（ie,k）;
end
%step3.计算分布压力的等效节点力，并集成到整体节点力向量中
[df_number,dummy]=size（gDF）;
foridf=1:
1:
df_number
enf=EquivalentDistPressure（gDF（idf,1）,gDF（idf,2）,gDF（idf,3）,gDF（idf,4））;
i=gDF（idf,1）;
j=gDF（idf,2）;
f（（i-1）*2+1:
（i-1）*2+2,1）=f（（i-1）*2+1:
（i-1）*2+2,1）+enf（1:
2）;
f（（j-1）*2+1:
（j-1）*2+2,1）=f（（j-1）*2+1:
（j-1）*2+2,1）+enf（3:
4）;
end
%step4.处理约束条件，修改刚度矩阵和节点力向量。
采用乘大数法
[bc1_number,dummy]=size（gBC1）;
foribc=1:
1:
bc1_number
g=gBC1（ibc,2）;
f（g）=gBC1（ibc,3）*gK（g,g）*1e15;
gK（g,g）=gK（g,g）*1e15;
end
%step5.求解方程组，得到节点位移向量
gDelta=gK\f;
%step6.计算每个单元的结点应力
gElementStress=zeros（element_number,5,3）;
delta=zeros（8,1）;
forie=1:
element_number
xi=[-111-1];
eta=[-1-111];
forn=1:
4
B=MatrixB（ie,xi（n）,eta（n））;
D=MatrixD（ie,1）;
delta（1:
2:
7）=gDelta（gElement（ie,1:
4）*2-1）;
delta（2:
2:
8）=gDelta（gElement（ie,1:
4）*2）;
sigma=D*B*delta;
gElementStress（ie,n,:
）=sigma;
end
end
return
functionk=StiffnessMatrix（ie）
%计算平面应力等参数单元的刚度矩阵
%输入参数：
%ie--单元号
%返回值：
%k--单元刚度矩阵
%说明：
%用高斯积分法求解平面等参数单元的刚度矩阵
k=zeros（8,8）;
D=MatrixD（ie,1）;
%3x3高斯积分点和权系数
%x=[-0.774596669241483,0.0,0.774596669241483];
%w=[0.555555555555556,0.888888888888889,0.555555555555556];
%2x2高斯积分点和权系数
x=[-0.577350269189626,0.577350269189626];
w=[1,1];
fori=1:
1:
length（x）
forj=1:
1:
length（x）
B=MatrixB（ie,x（i）,x（j））;
J=Jacobi（ie,x（i）,x（j））;
k=k+w（i）*w（j）*transpose（B）*D*B*det（J）;
end
end
return
functionD=MatrixD（ie,opt）
%计算单元的弹性矩阵D
%输入参数：
%ie---------单元号
%opt---------问题选项
%1==>平面应力
%2==>平面应变
%返回值：
%D---------弹性矩阵D
globalgElementgMaterial
E=gMaterial（gElement（ie,5）,1）;%弹性模量
mu=gMaterial（gElement（ie,5）,2）;%泊松比
ifopt==1%平面应力的弹性常数
A1=mu;
A2=（1-mu）/2;
A3=E/（1-mu^2）;
else%平面应变的弹性常数
A1=mu/（1-mu）;
A2=（1-2*mu）/2/（1-mu）;
A3=E*（1-mu）/（1+mu）/（1-2*mu）;
end
D=A3*[1A10;
A110;
00A2];
return
functionB=MatrixB（ie,xi,eta）
%计算单元的应变矩阵B
%输入参数：
%ie---------单元号
%xi,eta-----局部坐标
%返回值：
%B---------在局部坐标处的应变矩阵B
[N_x,N_y]=N_xy（ie,xi,eta）;
B=zeros（3,8）;
fori=1:
1:
4
B（1:
3,（2*i-1）:
2*i）=[N_x（i）0;
0N_y（i）;
N_y（i）N_x（i）];
end
return
function[N_x,N_y]=N_xy（ie,xi,eta）
%计算形函数对整体坐标的导数
%输入参数：
%ie---------单元号
%xi,eta-----局部坐标
%返回值：
%N_x-------在局部坐标处的形函数对x坐标的导数
%N_y-------在局部坐标处的形函数对y坐标的导数
J=Jacobi（ie,xi,eta）;
[N_xi,N_eta]=N_xieta（ie,xi,eta）;
A=J\[N_xi;N_eta];
N_x=A（1,:
）;
N_y=A（2,:
）;
return
function[N_xi,N_eta]=N_xieta（~,xi,eta）
%计算形函数对局部坐标的导数
%输入参数：
%ie---------单元号
%xi,eta-----局部坐标
%返回值：
%N_xi-------在局部坐标处的形函数对xi坐标的导数
%N_eta-------在局部坐标处的形函数对eta坐标的导数
x=[-1,1,1,-1];
e=[-1,-1,1,1];
N_xi=zeros（1,4）;
N_eta=zeros（1,4）;
N_xi
（1）=x
（1）*（1+e
（1）*eta）/4;
N_eta
（1）=e
（1）*（1+x
（1）*xi）/4;
N_xi
（2）=x
（2）*（1+e
（2）*eta）/4;
N_eta
（2）=e
（2）*（1+x
（2）*xi）/4;
N_xi（3）=x（3）*（1+e（3）*eta）/4;
N_eta（3）=e（3）*（1+x（3）*xi）/4;
N_xi（4）=x（4）*（1+e（4）*eta）/4;
N_eta（4）=e（4）*（1+x（4）*xi）/4;
return
functionJ=Jacobi（ie,xi,eta）
%计算雅克比矩阵
%输入参数：
%ie---------单元号
%xi,eta-----局部坐标
%返回值：
%J-------在局部坐标（xi,eta）处的雅克比矩阵
globalgNodegElement
x=gNode（gElement（ie,1:
4）,1）;
y=gNode（gElement（ie,1:
4）,2）;
[N_xi,N_eta]=N_xieta（ie,xi,eta）;
x_xi=N_xi*x;
x_eta=N_eta*x;
y_xi=N_xi*y;
y_eta=N_eta*y;
J=[x_xi,y_xi;x_eta,y_eta];
return
functionAssembleStiffnessMatrix（ie,k）
%把单元刚度矩阵集成到整体刚度矩阵
%输入参数:
%ie---单元号
%k---单元刚度矩阵
%返回值:
%无
globalgElementgK
fori=1:
1:
4
forj=1:
1:
4
forp=1:
1:
2
forq=1:
1:
2
s=（i-1）*2+p;
t=（j-1）*2+q;
S=（gElement（ie,i）-1）*2+p;
T=（gElement（ie,j）-1）*2+q;
gK（S,T）=gK（S,T）+k（s,t）;
end
end
end
end
return
functionenf=EquivalentDistPressure（i,j,pi,pj）
%计算线性分布压力的等效节点力
%输入参数:
%i,j----------结点号
%pi,pj--------跟结点号对应的压力值
%返回值:
%enf-----------等效节点力向量
globalgNode
%获取结点坐标
xi=gNode（i,1）;
xj=gNode（j,1）;
yi=gNode（i,2）;
yj=gNode（j,2）;
sig（1:
4,1）=0;
%2x2高斯积分点和权系数
x=[-0.577350269189626,0.577350269189626];
w=[1,1];
fori=1:
1:
length（x）
N1=（1-x（i）
升级会员 