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

2021-03-11 19:04:13 阅读：643 来源： 互联网

标签：yi 有限元 0.3750 0.1875 KK 1.1250 薄板 matlab NU

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

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

标签：yi,有限元,0.3750,0.1875,KK,1.1250,薄板,matlab,NU
来源： https://blog.csdn.net/w1029126938/article/details/114676929

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

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

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

问题描述

要用到的函数

解答步骤