弹力有限元作业.docx
《弹力有限元作业.docx》由会员分享，可在线阅读，更多相关《弹力有限元作业.docx（15页珍藏版）》请在冰豆网上搜索。
弹力有限元作业
弹力有限元作业
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=
0.1667
>>t=1
t=
1
>>k1=LinearTriangleElementStiffness（E,NU,t,0,1,1,2,0,2,1）
k1=
0.214300-0.2143-0.21430.2143
00.5143-0.085700.0857-0.5143
0-0.08570.51430-0.51430.0857
-0.2143000.21430.2143-0.2143
-0.21430.0857-0.51430.21430.7286-0.3000
0.2143-0.51430.0857-0.2143-0.30000.7286
>>k2=LinearTriangleElementStiffness（E,NU,t,1,2,0,1,1,1,1）
k2=
0.214300-0.2143-0.21430.2143
00.5143-0.085700.0857-0.5143
0-0.08570.51430-0.51430.0857
-0.2143000.21430.2143-0.2143
-0.21430.0857-0.51430.21430.7286-0.3000
0.2143-0.51430.0857-0.2143-0.30000.7286
>>k3=LinearTriangleElementStiffness（E,NU,t,0,0,1,1,0,1,1）
k3=
0.214300-0.2143-0.21430.2143
00.5143-0.085700.0857-0.5143
0-0.08570.51430-0.51430.0857
-0.2143000.21430.2143-0.2143
-0.21430.0857-0.51430.21430.7286-0.3000
0.2143-0.51430.0857-0.2143-0.30000.7286
（三）集成整体刚度矩阵
由于结构有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=
0.5143000-0.51430.08570-0.085700
00.2143000.2143-0.2143-0.2143000
0000000000
0000000000
-0.51430.2143000.7286-0.3000-0.21430.085700
0.0857-0.214300-0.30000.72860.2143-0.514300
0-0.214300-0.21430.21430.2143000
-0.08570000.0857-0.514300.514300
0000000000
0000000000
>>K=LinearTriangleAssemble（K,k2,1,4,2）
K=
0.72860-0.21430.2143-0.51430.08570-0.300000
00.72860.0857-0.51430.2143-0.2143-0.3000000
-0.21430.08570.7286-0.300000-0.51430.214300
0.2143-0.5143-0.30000.7286000.0857-0.214300
-0.51430.2143000.7286-0.3000-0.21430.085700
0.0857-0.214300-0.30000.72860.2143-0.514300
0-0.3000-0.51430.0857-0.21430.21430.7286000
-0.300000.2143-0.21430.0857-0.514300.728600
0000000000
0000000000
>>K=LinearTriangleAssemble（K,k3,5,2,4）
K=
0.72860-0.21430.2143-0.51430.08570-0.300000
00.72860.0857-0.51430.2143-0.2143-0.3000000
-0.21430.08571.2429-0.300000-1.02860.30000-0.0857
0.2143-0.5143-0.30000.9429000.3000-0.4286-0.21430
-0.51430.2143000.7286-0.3000-0.21430.085700
0.0857-0.214300-0.30000.72860.2143-0.514300
0-0.3000-1.02860.3000-0.21430.21431.4571-0.3000-0.21430.0857
-0.300000.3000-0.42860.0857-0.5143-0.30001.45710.2143-0.5143
000-0.214300-0.21430.21430.21430
00-0.08570000.0857-0.514300.5143
（四）引入边界条件
本题的边界条件如下：
U3x=U3y=U4x=U4y=U5x=U5y=U1x=U2x=0
F1y=-0.5F1x=0
（五）求解
>>k=K（2:
2:
4,2:
2:
4）
k=
0.7286-0.5143
-0.51430.9429
>>f=[-0.5;0]
f=
-0.5000
0
>>u=k\f
u=
-1.1159
-0.6087
所以：
v1=-1.1159F/Et,v2=-0.6087F/Et
（六）求支座反力
>>U=[0;-1.1159;0;-0.6087;0;0;0;0;0;0]
U=
0
-1.1159
0
-0.6087
0
0
0
0
0
0
>>F=K*U
F=
-0.1304
-0.5000
0.0870
-0.0000
-0.2391
0.2391
0.1522
0.2609
0.1304
0
所以，F3x=-0.2391F,F3y=0.2391F
F4x=0.1522F,F4y=0.2609F
F5x=0.1304F,F5y=0
3-11挡水墙的密度为
，厚度为
，如下图1示，水的密度为
，试求应力分量。
已知：
。
图1
采用ansys软件进行解答
解答步骤：
1，采用国际单位制
2，定义一个四节点平面，并定义其属性，混凝土坝的属性为密度2.8e3，弹性模量3.2e10，
泊松比为0.125，厚度为1.
3，网格划分
建立四个关键点并对该单元平面作进一步划分：
5*25个小单元格（40cmX40cm），并对各小单元进行节点编号、单元编号。
4，在底边施加约束，选择线上dof。
5，施加荷载。
如下图所示：
对本题共有两荷载作用于挡水墙，一是挡水墙自重（输入重力加速度9.8），一是水体对其压力。
本题将自重和水体压力组合为荷载组共同作用于挡水墙上，其中挡水墙以１ｍ的厚度计算。
模型求解
运行ansys求解。
求得应力应变图结果如下：
、
、
的应力图示如下图示：
观察应力图可大致得到如下规律：
图：
从上往下颜色越来越深，
越来越大；
图：
从中间往两边颜色越来越深，
越来越大在两边缘处取得各自的取大值；
图：
关于
轴对称，由上往下
越来越大。
所计算的结果输出另用txt结果输出。
理论值检验计算
建立如图1所示的直角坐标系，经理论推导得到其应力分量理论解如下：
；
；
；
采用matlab进行检验计算结果：
程序为：
%symsxysxsysxy;
x=input（'输入x坐标'）
y=input（'输入y坐标'）
%!
P2=1000;P1=2500;G=9.8;B=2.5;
sx=2*1000*9.8*x^3*y/8+3*1000*9.8*x*y/10-4*1000*9.8*x*y^3/8-2800*10*x;
sy=1000*9.8*x*（2*y^3/8-3*y/4-1/2）;
sxy=-1000*9.8*x^2*（3*y^2/8-3/8）-1000*9.8*y*（-y^3/8+3*y/20-2/80/y）;
sx
sy
sxy
计算结果与理论计算结果相差偏大，可能由于单元格划分过大引起。
升级会员 