有限元分析课程作业Word文档格式.docx
《有限元分析课程作业Word文档格式.docx》由会员分享，可在线阅读，更多相关《有限元分析课程作业Word文档格式.docx（33页珍藏版）》请在冰点文库上搜索。
将式（1-1）代入节点条件式（1-2）中，可求出式（
1-1）
中待定系数，即
a。
a’
a2
bo
b1
b2
U1X1
1
—u2X2
2A
U4X4
U1
U2
U4
—1
X1
X2
X4
GW
y1
y2
y4
a2v2
（ay
a?
U2a3U4）
（dm
b2U2b3U4）
a3V4）（1-6）
药冲1b2V2b3V4）（1-7）
2a（C1V1C2V2C3V4）（1-8）
2A
（1-3）
（1-4）
（1-5）
在式（1-3）~式（1-8）中
1X1y1
1X2y2
1X4y
ai
x4y4
上式中的符号（1,2,3）表示下标轮换，如12,23,31同时更
换。
将单元①各节点的位置坐标&
0,y10必1,y20,&
0,y4
代入得
将系数式（1-3）〜式（1-8）式代入（1-1）中，重写位移函
数，并以节点位移的形式进行表示，有
u（x,y）N^x,y）uN2（x,y）U2v（x,y）N1（x,y）vN2（x,yM
（1-12）
位移函数式（1-11）写成矩阵形式，有
u
（1）（x，y）u!
x，y）
（26）v（x,y）
V1
v
V4
（1-13）
对于单元②，
过程同上，
形状函数矩阵
⑵/、
H（x，y）
（1-14）
位移函数
（2）
u（x,y）
（26）
u（x,y）
v（x,y）
1+x
U3
V3
V
（1-15）
c.单元的应变场描述
对于单元①，
应变函数
xx
⑴（x,y）
（31）
yy
xy
u
x
uv
yx
u（x,y）v（x,y）
Ui
Vi
x0y0u2
V2
-101000
=0-10001
-1-10110
其中几何矩阵
（1）/、
B（x,y）
（3"
6）
0-10001
-1-10110
（1-17）
对于单元②，应变方程
—0
1xy01x01y0
y01xy01x01y
-1
（1-18）
B（x,y）0
（36）
0-1000
1000-1
10-1-10
（1-19）
d.单元的应力场描述
（x,y,Z）yy
"
~2
=D
（33）
（1）
（3
1）
（1-20）
其中，弹性系数矩阵
D）
2
107
~~2
1-1
=3.75
106
（1-21
（1-22）
其中应力函数矩阵
s
（1）为
（1）
（1）
（1）
D）（B）=375
3
-3
1061
1=3.75
106-1
（1-23）
应力方程⑴为
（1）=3.75106-1
q=3.75
（61）
u2v2
v1
Pa
u4
v4
（1-24）
对于单元②，过程同上
弹性系数矩阵D
（2）为
D=3.7510
31
613
00
应力函数矩阵
s
（2）为
-30
S
（2）=13
-10
（36）11
0-1
（1-25）
（1-26）
应力方程
（2）为
u3
（2）6
=3.751061
v3
4Pav4u2
（1-27）
e.单元的势能表达
K
（1）是单元刚度矩阵，即
K
（1）
（1）B
（1）TD
（1）B
（1）d
（1）B
（1）TD
（1）B
（1）dAt（1-28）
其中薄板厚度t0.1m。
将式（1-17）、式（1-21）代入式（1-29）,
得到单元①的刚度阵
K
（1）B
（1）TD
（1）B
（1）tA3.75
0.1
计算得
同理，
u1
u2
v2
7.5
3.75
5.625
1.875
（1）55.625
K
（1）105
-1.875
得到单元②的刚度阵为
（2）55.625
（2）105
K
将两个单元按节点位移所对应的位置进行组装，得到
总体刚度矩阵为
KKK
（1）K
（2）
节点力
F
FPx3FPx3
FRx4
10.110103（FRx101010Frx40）
0.5103（Frx101010FRx40）N
系统的势能
UW=1qTKq-FTq（计算结果在下面呈现）
（4）边界条件的处理及方程求解
边界条件为UiViU4V40。
因此，将针对节点2和节点3的位移求解，节点2和节点3对应总体刚度阵KK中的第3行到第6行、第3列到第6列，贝U需从KK88中提出，置给k，然后生成对应的载荷列阵p,再采用高斯消去法进行求解。
>
k=KK（3:
6,3:
6）;
p=[500;
0;
500;
0];
u=k\p
u=[将列排成了行]
再计算支反力。
在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力；
先将上面得到的位移结果与边界条件的节点位移进行组合，得到整体的位移列阵U（81），再代回原整体刚度方程，计算出所有的节点力，按照位置关系找出对应的支反力。
U=[0;
u;
0]
U=0000
[将列排成了行]
P=KK*U
P=-50050005000-500
[将列排成了行]所以，节点1的支反力为FRx1500N,FRy1-176.4706N，节点2的支反力为FRx2500N,FRy2176.4706N。
根据已求得的位移和支反力计算系统的势能。
A=*U'
*KK*U-P'
*U
A=
（5）结果分析上述支反力计算结果满足静力平衡，验证了以上求解过程及MATLAB算法的正确性。
2、基于四节点四边形单元的有限元建模及分析
（1）结构的离散化与编号
如图3所示一个4节点矩形单元，单元的节点位移共有8个自由度（DOF。
节点编号为123,4，各自的位置坐标为x*，i123,4，各个节点的位移（分别沿x方向和y方向）
为ui,vi,i1,2,3,4。
4o
图3方案二：
使用一个4节点矩形单元
（2）局部坐标系下单元的描述a.单元的几何和节点描述
采用无量纲坐标
=Xy
a'
b
其中a0.5,b0.5。
则单元四个节点的几何位置为
i,i1
1,21
1,31
1,41
将所有节点上的位移组成一个列阵，记作q；
同样，将
所有节点上的各个力也组成一个列阵，记作F，则有
q
（81）
（F）
（U1V1U2V2U3V3U4V4）T
（FX1Fy1Fx2Fy2Fx3Fy3Fx4Fy4）
设位移函数为
u（x,y）aoaxa?
yasxy
v（x,y）bo
b|Xb2yb3xy
由节点条件，
在xx’yyi处，有
u（Xi,yJUi
i123,4
将位移试函数代入节点条件中，求出待定系数
和bo,bl,b2,b3，再代入位移函数中，整理后得
u（x,y）
Ndx,y）u1
n
2（x,y）U2
N3（x,y）U3
N4（x,y）U4
N,x,y）v1
2（x,y）v2
N3（x,y）v3
N4（x,y）v4
其中
M（x,y）
—12x
4
2y
N2（x,y）
112x
Wx,y）
N4（x,y）
如以无量纲坐标来表达，可写成
Ni-1i1i,i123,4
将11,1121,21，31,3141,4
a0,ai,a2，a3
1带入上式
得到形状函数矩阵
N1
—
N2
1-
N3
N4
写成矩阵形式，有
N10N20N30N40v2
0N10N20N30N4u3Nq
C.单元的应变场描述
单元应变为
其中几何矩阵B（x,y）为
应力表达式为
（3i）D）（3i）dm羸儿其中，应力函数矩阵SDB。
e.单元的势能表达
以上已将单元的三大基本变量U,,用基于节点的位移
其中K为4节点矩形单元的刚度矩阵，即
其中，t为薄板的厚度，t0.1m，上式的各个字块矩阵为
检0.1bTdBdr，s1,2,3,4
f.单元刚度阵及刚度方程
单元刚度阵在上面已经列出。
将单元的势能对节点位移
q取一阶极值，可得到单元的刚度方程
处理方法与3节点三角形单元一致，利用上述求解程序具有的可移植性，简化了求解过程。
k=K（3:
u=*
[将列排成了行]
同样注意按照位置关系找出对应的支反
力。
0000
P=K*U
P=
-50050005000-500
[将列排成了行]所以，节点1的支反力为FRx1500N,FRy1-111.1111N，节点
2的支反力为FRx2500N,FRy2111.1111N。
根据已求得的位移和支反力计算系统的势能。
*K*U-P'
（5）结果分析
基于4节点矩形单元计算出的势能小于基于3节点三角
形单元计算出的结果，若将该系统分为更多的单元，计算精度也将提高。
3.两种方案的比较与分析从以上计算可以看出，用三角形单元计算时，由于形函数是完全一次式，因而其应变场在单元内均为常数；
而对于四边形单元，其形函数带有二次式，计算得到的应变场和应力场都是坐标的一次函数，但不是完全的一次函数，对提高精度有一定作用；
根据最小势能原理，势能越小，则整体计算精度越高，比较两种单元计算得到的系统势能，可看出，在相同的节点自由度情况下，矩形单元的计算精度要比三角形单元高。
三、基于MATLAB勺编程实现
1.基于3节点三角形单元的有限元编程实现
（1）程序编写说明
Triangle2D3Node_Stiffness（E,NU,t,xi,yi,xj,yj,xm,ym,ID）
该函数计算单元的刚度矩阵，输入模量E,泊松比NU,
厚度t，三个节点i，j，m平面问题性质指示参数（1为平面应力，2为平面应变），输出单元刚度矩阵k（66）。
Triangle2D3Node_Assemble（KK,k,i,j,m）
该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i，j，m输出整体刚度矩阵KKo
Triangle2D3Node_Stress（E,NU,xi,yi,xj,yj,xm,ym,u
ID）
该函数计算单元的应力，输入弹性模量E，泊松比NU,
厚度t，三个节点i，j，m平面问题性质指示参数（1为平面应力，2为平面应变），单元的位移列阵u（61），输出单元的应力，由于它为常应力单元，则单元的应力分量为Sx，Sy，Sxy。
2）程序清单
%%%%%Triangle2D3Node%%%begin%%%%%%%function
k=Triangle2D3Node_Stiffness（E,NU,t,xi,yi,xj,yj,xm,y
m,ID）
%该函数计算单元的刚度矩阵
%俞入弹性模量E、泊松比NU和厚度t
%俞入3个节点i,j,m的坐标xi,yi,xj,yj,xm,ym
%俞入平面问题性质指示参数ID（1位平面应力，2为平
面应变）
%俞入单元刚度矩阵k（6*6）
%
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）;
ifID==1
D=（E/（1-NU*NU））*[1NU0;
NU10;
00（1-NU）/2];
elseifID==2
D=（E/（1+NU）/（1-2*NU））*[1-NUNU0;
NU1-NU0;
0（1-2*NU）/2];
end
k=t*A*B'
*D*B;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionz=Triangle2D3Node_Assemble（KK,k,i,j,m）%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k
%输入单元的节点编号i，j，m
%输入整体刚度矩阵KKrf
DOF
（1）=2*i-1;
DOF
（2）=2*i;
DOF（3）=2*j-1;
DOF（4）=2*j;
DOF（5）=2*m-1;
DOF（6）=2*m;
forn1=1:
6
forn2=1:
KK（DOF（n1）,DOF（n2））=KK（DOF（n1）,DOF（n2））+k（n1,n2）;
z=KK;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%function
stress=Triangle2D3Node_Stress（E,NU,xi,yi,xj,yj,xm,y
m,u,ID）
%该函数计算单元的应力
%输入平面问题性质指示参数ID（1位平面应力，2为平
面应变），单元的位移列阵u（6*1）
%俞出单元的应力stress（3*1），由于它为常应力单元，则单元的应力分量Sx，Sy，Sxy
A=（xi*xj（ym-yi）+xm*（yi-yj））/2;
betai=yj-ym;
betaj=ym-yi;
betam=yi-yj;
gammai=xm-xj;
gammaj=xi-xm;
B=[batai0bataj0betam0;
gammaibetaigammajbetajgammambetam]/（2*A）;
elseifID==2
endstress=D*B*u;
%%%%%%%%%%%%%%%%%%%%%%%%%%
E=1E7;
NU=1/3;
t=;
c=Triangle2D3Node_Stiffness（E,NU,t,0,0,1,0,0,1,2）
CC=zeros（8,8）;
CC=Triangle2D3Node_Assemble（KK,k1,1,2,4）；
CC=Triangle2D3Node_Assemble（KK,k1,3,4,2）
k1=Triangle2D3Node_Stiffness（E,NU,t,0,0,1,0,0,1,1）
KK=zeros（8,8）;
KK=Triangle2D3Node_Assemble（KK,k1,1,2,4）
KK=Triangle2D3Node_Assemble（KK,k1,3,4,2）
U=[0;
P=KK*U
（3）计算结果
①应变
CC=
+05*
②位移
U=
3节点力
0其中，节点1的支反力为
FRx1500N,FRy1-176.4706N，节点2的支反力为FRx2500N,FRy2176.4706N。
4势能A=
5单元刚度阵
KK=
00
2.基于四节点四边形单元的有限元建模及分析
（1）程序编写说明
Quad2D4Node_Stiffness（E,NU,h,xi,yi,xj,yj,xm,ym,xp,y
p,ID）
该函数计算单元的刚度矩阵，输入模量E,泊松
升级会员 