平面问题的有限元方法86作业集中力平面问题.docx
《平面问题的有限元方法86作业集中力平面问题.docx》由会员分享，可在线阅读，更多相关《平面问题的有限元方法86作业集中力平面问题.docx（10页珍藏版）》请在冰点文库上搜索。
平面问题的有限元方法86作业集中力平面问题
题目：
为一个受到集中力P作用的结构，t=1,试按平面应力问题计算，采用三角形单元，求出节点位移。
假设E为弹性模量，泊松比
=
，
t=1cm,平面应力问题
将结构划成单元，并对其编号。
（1）定义单元
单元定义和有关数据见下表。
编程为：
betai=yj-ym;
betaj=ym-yi;
betam=yi-yj;
gammai=xm-xj;
gammaj=xi-xm;
gammam=xj-xi;
依次代入节点的坐标求得下表：
i
j
m
△
单
元
号
1
（0,1）
2
（0，0）
3
（1，1）
-101
1-10
2
（0，0）
4
（1，0）
3
（1，1）
-110
0-11
3
（1，1）
4
（1，0）
5
（2，0）
0-11
1-10
3
（1，1）
5
（2，0）
6
（2，1）
-101
0-11
（2）求各单元的刚度矩阵
从上表可以看出，4个单元刚度不同
单元
的刚度矩阵为：
其中子矩阵表达式为：
[y1]=LinearTriangleElementStiffness（1,0.1667,1,0,1,0,0,1,1,1）
[y2]=LinearTriangleElementStiffness（1,0.1667,1,0,0,1,0,1,1,1）
[y3]=LinearTriangleElementStiffness（1,0.1667,1,1,1,1,0,2,0,1）
[y4]=LinearTriangleElementStiffness（1,0.1667,1,1,1,2,0,2,1,1）
（3）总体刚度矩阵组装
单元编号
①
②
③
④
整体编码
1,2,3
2,4,3
3,4,5
3,5,6
局部编码
以整编码体表示的单元刚度矩阵
总体刚度矩阵为
运行如下程序：
得到：
总体刚度矩阵
（4）求总体载荷列阵
式子中
限制了结构的刚体位移。
（5）引入几何边界条件为：
根据划0置1法,经过处理的总体矩阵为：
k（:
[1234]）=[]；k（[1234],:
）=[]
求解上述方程为：
假设P=[11111111]'
C=[u3v3u4v4u5v5u6v6]'
c=inv（k）*A
%%%%%%%%%%%%%%%求整体刚度矩阵
clc
clearall
symsEuta%%定义变量
%%%%%%%%%%%%%%%%%%%%%%%%写出单元刚度矩阵,单元面积为a^2/2
%%%%%%%%%%%%单元1
u=0.1667;E=1;t=1;a=1;
bi=-1;ci=a;
bj=0;cj=-a;
bm=1;cm=0;
mianji=a^2/2;
B1=1/2/mianji*[bi0bj0bm0
0ci0cj0cm
cibicjbjcmbm];
D=E/（1-u^2）*[1u0
u10
00（1-u）/2];%%%弹性矩阵D
k1=transpose（B1）*D*B1*t*mianji
%%%%%%%%%%%%单元2
bi=-1;ci=0;
bj=1;cj=-1;
bm=0;cm=1;
mianji=a^2/2;
B2=1/2/mianji*[bi0bj0bm0
0ci0cj0cm
cibicjbjcmbm];
D=E/（1-u^2）*[1u0
u10
00（1-u）/2];%%%弹性矩阵D
k2=transpose（B2）*D*B2*t*mianji
%%%%%%%%%%%%单元3
bi=0;ci=1;
bj=-1;cj=-1;
bm=1;cm=0;
mianji=a^2/2;
B3=1/2/mianji*[bi0bj0bm0
0ci0cj0cm
cibicjbjcmbm];
D=E/（1-u^2）*[1u0
u10
00（1-u）/2];%%%弹性矩阵D
k3=transpose（B3）*D*B3*t*mianji
%%%%%%%%%%%%单元4
bi=-1;ci=0;
bj=0;cj=-1;
bm=1;cm=1;
mianji=a^2/2;
B4=1/2/mianji*[bi0bj0bm0
0ci0cj0cm
cibicjbjcmbm];
D=E/（1-u^2）*[1u0
u10
00（1-u）/2];%%%弹性矩阵D
k4=transpose（B4）*D*B4*t*mianji
%%%%%%%%%%%%%%%%%%%%%%%%%单元刚度矩阵再分解
k1_11=k1（1:
2,1:
2）;
k1_12=k1（1:
2,3:
4）;
k1_13=k1（1:
2,5:
6）;
k1_21=k1（3:
4,1:
2）;
k1_22=k1（3:
4,3:
4）;
k1_23=k1（3:
4,5:
6）;
k1_31=k1（5:
6,1:
2）;
k1_32=k1（5:
6,3:
4）;
k1_33=k1（5:
6,5:
6）;
%%%%%%%%%%%%%
k2_22=k2（1:
2,1:
2）;
k2_24=k2（1:
2,3:
4）;
k2_23=k2（1:
2,5:
6）;
k2_42=k2（3:
4,1:
2）;
k2_44=k2（3:
4,3:
4）;
k2_43=k2（3:
4,5:
6）;
k2_32=k2（5:
6,1:
2）;
k2_34=k2（5:
6,3:
4）;
k2_33=k2（5:
6,5:
6）;
%%%%%%%%%%%%%
k3_33=k2（1:
2,1:
2）;
k3_34=k2（1:
2,3:
4）;
k3_35=k2（1:
2,5:
6）;
k3_43=k2（3:
4,1:
2）;
k3_44=k2（3:
4,3:
4）;
k3_45=k2（3:
4,5:
6）;
k3_53=k2（5:
6,1:
2）;
k3_54=k2（5:
6,3:
4）;
k3_55=k2（5:
6,5:
6）;
%%%%%%%%%%%%%
k4_33=k2（1:
2,1:
2）;
k4_35=k2（1:
2,3:
4）;
k4_36=k2（1:
2,5:
6）;
k4_53=k2（3:
4,1:
2）;
k4_55=k2（3:
4,3:
4）;
k4_56=k2（3:
4,5:
6）;
k4_63=k2（5:
6,1:
2）;
k4_65=k2（5:
6,3:
4）;
k4_66=k2（5:
6,5:
6）;
%%%%%%%%%%%单元刚度矩阵组装，形成整体刚度矩阵%%%%%%%
k_11=k1_11;k_12=k1_12;k_13=k1_13;k_14=zeros（2,2）;
k_15=zeros（2,2）;k_16=zeros（2,2）;k_21=k1_21;k_22=k1_22+k2_22;
k_23=k1_23+k2_23;
k_24=k2_24;k_25=zeros（2,2）;k_26=zeros（2,2）;k_31=k1_31;
k_32=k1_32+k2_32;k_33=k1_33+k2_33+k3_33+k4_33;
k_34=k2_34+k3_33;
k_35=k3_35+k4_35;
k_36=k4_36;k_41=zeros（2,2）;k_42=k2_42;k_43=k2_43+k3_43;k_44=k2_44+k3_44;
k_45=k3_45;k_46=zeros（2,2）;k_51=zeros（2,2）;k_52=zeros（2,2）;
k_53=k3_53+k4_53;k_54=k3_54;k_55=k3_55+k3_55+k4_55;k_56=k4_56;
k_61=zeros（2,2）;k_62=zeros（2,2）;k_63=k4_63;k_64=zeros（2,2）;
k_65=k4_65;
k_66=k4_66;
%%%%
k=[k_11k_12k_13k_14k_15k_16;
k_21k_22k_23k_24k_25k_26;
k_31k_32k_33k_34k_35k_36;
k_41k_42k_43k_44k_45k_46;
k_51k_52k_53k_54k_55k_56;
k_61k_62k_63k_64k_65k_66]
矩阵压缩：
k（:
[1234]）=[]
k（[1234],:
）=[]
升级会员 