有限元程序设计大作业DOCWord文档格式.docx
《有限元程序设计大作业DOCWord文档格式.docx》由会员分享，可在线阅读，更多相关《有限元程序设计大作业DOCWord文档格式.docx（36页珍藏版）》请在冰豆网上搜索。
n/2
gNode（（n+1）*m+（n/2+1）+i,:
）=[b/tan（dfi*（i+n/2））,b];
gElement=zeros（m*n,5）;
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
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
6.0286
4
6.1114
3
3.6264
6
4.0315
3.2289
7
3.5248
5
2.9094
8
3.3169
2.7122
9
3.2098
2.5751
10
3.1470
2.8527
11
3.1069
2.7762
3.0170
2.7102
12
3.0604
17
16
2.6707
3.0071
26
25
2.6589
20
2.9911
37
36
2.6545
24
2.9851
50
49
2.6526
28
2.9824
表中，
为径向单元的数目。
分别为与之对应的应力集中系数（
）。
为半板宽b与孔半径r的比值，即
。
解析解中
的取值为3.下图（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
%计算结点坐标
dr=（b-r）/（m+m^2/2+m^3/6）;
dfi=（pi/2）/n;
gNode=zeros（（n+1）*（m+1）,2）;
fori=1:
）=[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
%单元定义
gElement=zeros（m*n,5）;
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];
gElement（:
gMaterial=[E,mu,0];
%定义约束条件
gBC1=zeros（（m+1）*2,3）;
m+1
gBC1（i,1:
2）=[（n+1）*（i-1）+1,（（n+1）*（i-1）+1）*2];
gBC1（m+1+i,1:
2）=[（n+1）*i,（（n+1）*i-1）*2+1];
gBC1（:
3）=0;
%定义分布压力
gDF=zeros（n/2,4）;
n/2
gDF（i,1:
2）=[（n+1）*m+i,（n+1）*m+i+1];
gDF（:
3）=p;
4）=p;
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:
element_number
fprintf（sprintf（'
计算单元刚度矩阵并集成，当前单元:
%d\n'
ie））;
k=StiffnessMatrix（ie）;
AssembleStiffnessMatrix（ie,k）;
%step3.计算分布压力的等效节点力，并集成到整体节点力向量中
[df_number,dummy]=size（gDF）;
foridf=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）;
%step4.处理约束条件，修改刚度矩阵和节点力向量。
采用乘大数法
[bc1_number,dummy]=size（gBC1）;
foribc=1:
bc1_number
g=gBC1（ibc,2）;
f（g）=gBC1（ibc,3）*gK（g,g）*1e15;
gK（g,g）=gK（g,g）*1e15;
%step5.求解方程组，得到节点位移向量
gDelta=gK\f;
%step6.计算每个单元的结点应力
gElementStress=zeros（element_number,5,3）;
delta=zeros（8,1）;
forie=1:
xi=[-111-1];
eta=[-1-111];
forn=1:
B=MatrixB（ie,xi（n）,eta（n））;
D=MatrixD（ie,1）;
delta（1:
2:
7）=gDelta（gElement（ie,1:
4）*2-1）;
delta（2:
8）=gDelta（gElement（ie,1:
4）*2）;
sigma=D*B*delta;
gElementStress（ie,n,:
）=sigma;
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];
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）;
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）;
D=A3*[1A10;
A110;
00A2];
functionB=MatrixB（ie,xi,eta）
%计算单元的应变矩阵B
%xi,eta-----局部坐标
%B---------在局部坐标处的应变矩阵B
[N_x,N_y]=N_xy（ie,xi,eta）;
B=zeros（3,8）;
B（1:
3,（2*i-1）:
2*i）=[N_x（i）0;
0N_y（i）;
N_y（i）N_x（i）];
function[N_x,N_y]=N_xy（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,:
function[N_xi,N_eta]=N_xieta（~,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;
functionJ=Jacobi（ie,xi,eta）
%计算雅克比矩阵
%J-------在局部坐标（xi,eta）处的雅克比矩阵
globalgNodegElement
x=gNode（gElement（ie,1:
4）,1）;
y=gNode（gElement（ie,1:
4）,2）;
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];
functionAssembleStiffnessMatrix（ie,k）
%把单元刚度矩阵集成到整体刚度矩阵
%输入参数:
%ie---单元号
%k---单元刚度矩阵
%返回值:
globalgElementgK
forp=1:
forq=1:
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）;
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;
N1=（
升级会员 