有限元分析简单实例之平面矩形薄板（matlab）

2024-01-18 22:43| 来源: 网络整理| 查看: 265

有限元分析简单实例之平面矩形薄板（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

可得到支反力在这里插入图片描述以上内容源自《有限元分析及应用》清华大学曾攀老师主讲

【本文地址】

公司简介

联系我们

今日新闻

推荐新闻

专题文章