有限元分析简单实例之平面矩形薄板(matlab)
问题描述
对于如图所示的一个平面矩形薄板结构,施加如右图所示的几个方向力,对其进行有限元分析,计算各个节点的位移及支座反力。(其中F是合力,E是弹性模量,μ是泊松比,t是厚度)
要用到的函数
(1)计算单元的刚度矩阵
function k = Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)
% 计算单元的刚度矩阵
% 输入弹性模量E,泊松比NU,厚度t
% 输入三个节点的坐标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 = [betai 0 betaj 0 betam 0;
0 gammai 0 gammaj 0 gammam;
gammai betai gammaj betaj gammam betam]/(2*A);
if ID==1
D=(E/(1-NU*NU))*[1 NU 0;NU 1 0;0 0 (1-NU)/2];
elseif ID==2
D =(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];
end
k = t*A*B'*D*B;
end
(2)进行单元刚度矩阵的组装
function z = Triangle2D3Node_Assembly(KK,k,i,j,m)
% 该函数进行单元刚度矩阵的组装
% 输入单元刚度矩阵k,单元节点编号i,j,m
% 输出整体刚度矩阵KK
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;
for n1=1:6
for n2=1:6
KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
end
end
z = KK;
end
(3)计算单元的应力
function stress = Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)
% 计算单元的应力
% 输入弹性模量E,泊松比NU,厚度t
% 输入三个节点的坐标xi,yi,xj,yj,xm,ym
% 输入平面问题性质指标参数ID(1为平面应力,2为平面应变),单元的位移矩阵u(6*1)
% 输出单元的应力stress(3*1),由于它为常应力单元,单元的应力分量为Sx,Sy,Sxy
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 = [betai 0 betaj 0 betam 0;
0 gammai 0 gammaj 0 gammam;
gammai betai gammaj betaj gammgm betam]/(2*A);
if ID==1
D=(E/(1-NU*NU))*[1 NU 0;NU 1 0;0 0 (1-NU)/2];
elseif ID==2
D =(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];
end
stress = D*B*u;
end
解答步骤
(1)变量输入及单元刚度矩阵的求解
输入弹性模量E,泊松比NU,薄板厚度t,指示参数ID,针对单元1和2,分辨求其刚度矩阵。
>> E= 1e7;
>> NU = 1/3;
>> t =0.1;
>> ID = 1;
>>> k1 = Triangle2D3Node_Stiffness(E,NU,t,2,0,0,1,0,0,ID)
k1 =
1.0e+06 *
0.2813 0 0 0.1875 -0.2813 -0.1875
0 0.0938 0.1875 0 -0.1875 -0.0938
0 0.1875 0.3750 0 -0.3750 -0.1875
0.1875 0 0 1.1250 -0.1875 -1.1250
-0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750
-0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188
>> k2 = Triangle2D3Node_Stiffness(E,NU,t,0,1,2,0,2,1,ID)
k2 =
1.0e+06 *
0.2813 0 0 0.1875 -0.2813 -0.1875
0 0.0938 0.1875 0 -0.1875 -0.0938
0 0.1875 0.3750 0 -0.3750 -0.1875
0.1875 0 0 1.1250 -0.1875 -1.1250
-0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750
-0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188
(2)刚度矩阵的组装
一共四个节点,每个节点有x,y方向上2个*度,总*度为8,因此,结构的总的刚度矩阵为8*8。运用Triangle2D3Node_Assembly函数进行刚度矩阵的组装。
>> KK = zeros(8,8);
>> KK=Triangle2D3Node_Assembly(KK,k1,2,3,4);
>> KK=Triangle2D3Node_Assembly(KK,k1,3,2,1);
>> KK
KK =
1.0e+06 *
0.6563 0.3750 -0.3750 -0.1875 -0.2813 -0.1875 0 0
0.3750 1.2188 -0.1875 -1.1250 -0.1875 -0.0938 0 0
-0.3750 -0.1875 0.6563 0 0 0.3750 -0.2813 -0.1875
-0.1875 -1.1250 0 1.2188 0.3750 0 -0.1875 -0.0938
-0.2813 -0.1875 0 0.3750 0.6563 0 -0.3750 -0.1875
-0.1875 -0.0938 0.3750 0 0 1.2188 -0.1875 -1.1250
0 0 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750
0 0 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188
(3) 边界条件处理及位移矩阵的求解
节点3,4不动,因此针对节点1和节点2的位移进行求解,节点1和节点2的位移对应KK矩阵中前四行和前四列,用k表示。然后定义对应的载荷列阵p,采用高斯消去法进行求解
>> k = KK(1:4,1:4)
k =
1.0e+06 *
0.6563 0.3750 -0.3750 -0.1875
0.3750 1.2188 -0.1875 -1.1250
-0.3750 -0.1875 0.6563 0
-0.1875 -1.1250 0 1.2188
>> p =[0;-50000;0;-50000]
p =
0
-50000
0
-50000
>> u=k\p
u =
0.1877
-0.8992
-0.1497
-0.8422
(4)计算各个节点的节点力(包含支反力和所施加的外力)
>> U =[u;0;0;0;0]
U =
0.1877
-0.8992
-0.1497
-0.8422
0
0
0
0
>> P=KK*U
P =
1.0e+05 *
0.0000
-0.5000
-0.0000
-0.5000
-2.0000
-0.0702
2.0000
1.0702
可得到支反力
以上内容源自《有限元分析及应用》清华大学 曾攀老师主讲