弹力有限元作业.docx
《弹力有限元作业.docx》由会员分享，可在线阅读，更多相关《弹力有限元作业.docx（12页珍藏版）》请在冰豆网上搜索。
弹力有限元作业
6-4对下图所示的离散结构，试求结点1，2的位移及铰支座3，4，5的反力（按平面应力问题计算）
采用MATLAB进行运算程序：
一、计算弹性模量E、泊松比NU、厚度t、节点坐标为（xi,yi）、（xj,yj）、（xm,ym）的单元刚度矩阵。
p=1表明函数用于平面应力情况。
p=2表明函数用于平面应变情况。
functiony=LinearTriangleElementStiffness（E,NU,t,xi,yi,xj,yj,xm,ym,p）
A=（xi*（yj-ym）+xj*（ym-yi）+xm*（yi-yj））/2;
betai=yj-ym;
betaj=ym-yi;
betam=yi-yj;
gammai=xm-xj;
gammaj=xi-xm;
gammam=xj-xi;
B=[betai0betaj0betam0;
0gammai0gammaj0gammam;
gammaibetaigammajbetajgammambetam]/（2*A）;
ifp==1
D=（E/（1-NU*NU））*[1NU0;NU10;00（1-NU）/2];
elseifp==2
D=（E/（1+NU）/（1-2*NU））*[1-NUNU0;NU1-NU0;00（1-2*NU）/2];
end
y=t*A*B'*D*B;
二、计算整体刚度矩阵K。
functiony=LinearTriangleAssemble（K,k,i,j,m）
K（2*i-1,2*i-1）=K（2*i-1,2*i-1）+k（1,1）;
K（2*i-1,2*i）=K（2*i-1,2*i）+k（1,2）;
K（2*i-1,2*j-1）=K（2*i-1,2*j-1）+k（1,3）;
K（2*i-1,2*j）=K（2*i-1,2*j）+k（1,4）;
K（2*i-1,2*m-1）=K（2*i-1,2*m-1）+k（1,5）;
K（2*i-1,2*m）=K（2*i-1,2*m）+k（1,6）;
K（2*i,2*i-1）=K（2*i,2*i-1）+k（2,1）;
K（2*i,2*i）=K（2*i,2*i）+k（2,2）;
K（2*i,2*j-1）=K（2*i,2*j-1）+k（2,3）;
K（2*i,2*j）=K（2*i,2*j）+k（2,4）;
K（2*i,2*m-1）=K（2*i,2*m-1）+k（2,5）;
K（2*i,2*m）=K（2*i,2*m）+k（2,6）;
K（2*j-1,2*i-1）=K（2*j-1,2*i-1）+k（3,1）;
K（2*j-1,2*i）=K（2*j-1,2*i）+k（3,2）;
K（2*j-1,2*j-1）=K（2*j-1,2*j-1）+k（3,3）;
K（2*j-1,2*j）=K（2*j-1,2*j）+k（3,4）;
K（2*j-1,2*m-1）=K（2*j-1,2*m-1）+k（3,5）;
K（2*j-1,2*m）=K（2*j-1,2*m）+k（3,6）;
K（2*j,2*i-1）=K（2*j,2*i-1）+k（4,1）;
K（2*j,2*i）=K（2*j,2*i）+k（4,2）;
K（2*j,2*j-1）=K（2*j,2*j-1）+k（4,3）;
K（2*j,2*j）=K（2*j,2*j）+k（4,4）;
K（2*j,2*m-1）=K（2*j,2*m-1）+k（4,5）;
K（2*j,2*m）=K（2*j,2*m）+k（4,6）;
K（2*m-1,2*i-1）=K（2*m-1,2*i-1）+k（5,1）;
K（2*m-1,2*i）=K（2*m-1,2*i）+k（5,2）;
K（2*m-1,2*j-1）=K（2*m-1,2*j-1）+k（5,3）;
K（2*m-1,2*j）=K（2*m-1,2*j）+k（5,4）;
K（2*m-1,2*m-1）=K（2*m-1,2*m-1）+k（5,5）;
K（2*m-1,2*m）=K（2*m-1,2*m）+k（5,6）;
K（2*m,2*i-1）=K（2*m,2*i-1）+k（6,1）;
K（2*m,2*i）=K（2*m,2*i）+k（6,2）;
K（2*m,2*j-1）=K（2*m,2*j-1）+k（6,3）;
K（2*m,2*j）=K（2*m,2*j）+k（6,4）;
K（2*m,2*m-1）=K（2*m,2*m-1）+k（6,5）;
K（2*m,2*m）=K（2*m,2*m）+k（6,6）;
y=K;
三、计算单元位移矢量为u时的单元应力。
p=1表明函数用于平面应力情况。
p=2表明函数用于平面应变情况。
该函数返回单元应力矢量。
functiony=LinearTriangleElementStresses（E,NU,xi,yi,yj,xm,ym,p,u）
A=（xi*（yj-ym）+xj*（ym-yi）+xm*（yi-yj））/2;
betai=yj-ym;
betaj=ym-yi;
betam=yi-yj;
gammai=xm-xj;
gammaj=xi-xm;
gammam=xj-xi;
B=[betai0betaj0betam0;
0gammai0gammaj0gammam;
gammaibetaigammajbetajgammambetam]/（2*A）;
ifp==1
D=（E/（1-NU*NU））*[1NU0;NU10;00（1-NU）/2];
elseifp==2
D=（E/（1+NU）/（1-2*NU））*[1-NUNU0;NU1-NU0;00（1-2*NU）/2];
end
y=D*B*u
假设E=1，F=1，厚度t=1，杆件长度为1，最后算出的位移乘以F/Et，力乘以F，以课本5结点为坐标原点
（一）离散化
单元编号
结点i
结点j
结点m
1
2
3
4
1
5
1
4
2
3
2
4
（二）写出单元刚度矩阵
通过调用MATLAB的LinearTriangleElementStiffness函数，得到两个单元刚度矩阵k1,k2和k3，每个矩阵都是6×6的。
程序调用输出结果为：
>>E=1
E=
1
>>NU=1/6
NU=
>>t=1
t=
1
>>k1=LinearTriangleElementStiffness（E,NU,t,0,1,1,2,0,2,1）
k1=
00
00
00
00
>>k2=LinearTriangleElementStiffness（E,NU,t,1,2,0,1,1,1,1）
k2=
00
00
00
00
>>k3=LinearTriangleElementStiffness（E,NU,t,0,0,1,1,0,1,1）
k3=
00
00
00
00
（三）集成整体刚度矩阵
由于结构有5个结点，所以整体刚度矩阵式10×10的。
程序调用输出结果为：
>>K=zeros（10,10）
K=
0000000000
0000000000
0000000000
0000000000
0000000000
0000000000
0000000000
0000000000
0000000000
0000000000
>>K=LinearTriangleAssemble（K,k1,4,1,3）
K=
000000
000000
0000000000
0000000000
0000
0000
000000
000000
0000000000
0000000000
>>K=LinearTriangleAssemble（K,k2,1,4,2）
K=
0000
0000
0000
0000
0000
0000
0000
0000
0000000000
0000000000
>>K=LinearTriangleAssemble（K,k3,5,2,4）
K=
0000
0000
000
000
0000
0000
0
0
000000
000000
（四）引入边界条件
本题的边界条件如下：
U3x=U3y=U4x=U4y=U5x=U5y=U1x=U2x=0
F1y=F1x=0
（五）求解
>>k=K（2:
2:
4,2:
2:
4）
k=
>>f=[;0]
f=
0
>>u=k\f
u=
所以：
v1=Et,v2=Et
（六）求支座反力
>>U=[0;;0;;0;0;0;0;0;0]
U=
0
0
0
0
0
0
0
0
>>F=K*U
F=
0
所以，F3x=,F3y=
F4x=,F4y=
F5x=,F5y=0
3-11挡水墙的密度为
，厚度为
，如下图1示，水的密度为
，试求应力分量。
已知：
。
图1
采用ansys软件进行解答
解答步骤：
1，采用国际单位制
2，定义一个四节点平面，并定义其属性，混凝土坝的属性为密度，弹性模量，
泊松比为，厚度为1.
3，网格划分
建立四个关键点并对该单元平面作进一步划分：
5*25个小单元格（40cmX40cm），并对各小单元进行节点编号、单元编号。
4，在底边施加约束，选择线上dof。
5，施加荷载。
如下图所示：
对本题共有两荷载作用于挡水墙，一是挡水墙自重（输入重力加速度），一是水体对其压力。
本题将自重和水体压力组合为荷载组共同作用于挡水墙上，其中挡水墙以１ｍ的厚度计算。
模型求解
运行ansys求解。
求得应力应变图结果如下：
、
、
的应力图示如下图示：
观察应力图可大致得到如下规律：
图：
从上往下颜色越来越深，
越来越大；
图：
从中间往两边颜色越来越深，
越来越大在两边缘处取得各自的取大值；
图：
关于
轴对称，由上往下
越来越大。
所计算的结果输出另用txt结果输出。
理论值检验计算
建立如图1所示的直角坐标系，经理论推导得到其应力分量理论解如下：
；
；
；
采用matlab进行检验计算结果：
程序为：
%symsxysxsysxy;
x=input（'输入x坐标'）
y=input（'输入y坐标'）
%!
P2=1000;P1=2500;G=;B=;
sx=2*1000**x^3*y/8+3*1000**x*y/10-4*1000**x*y^3/8-2800*10*x;
sy=1000**x*（2*y^3/8-3*y/4-1/2）;
sxy=-1000**x^2*（3*y^2/8-3/8）-1000**y*（-y^3/8+3*y/20-2/80/y）;
sx
sy
sxy
计算结果与理论计算结果相差偏大，可能由于单元格划分过大引起。
精心搜集整理，只为你的需要
升级会员 