搜档网
搜档网
当前位置：搜档网 › 有限元作业：三角形单元求解Word版

有限元作业：三角形单元求解Word版

有限元作业：三角形单元求解Word版
有限元作业：三角形单元求解Word版

《有限元作业》

年级2015级

学院机电工程学院

专业名称

班级学号

学生姓名

2016年05月

如下图所示为一受集中力P作用的结构，弹性模量E为常量，泊松比V=1/6,厚度为I=1。按平面应力问题计算，运用有限元方法，分别采用三角形及四边形单元求解，求节点位移及单元应力（要求三角形单元数量不少于4个，四边形单元不少于2个）

图(一)

图（二）三角形单元求解

图（三）四边形单元求解

（1）如图划分三角形单元，工分成四个分别为④

（2）如图分别进行编号1、2、3、4、5、6，并建立坐标系

（3）编程进行求解，得出结果，其中假设力P=2000N

调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵

k1 =

1.0e+06 *

7.2857 -3.0000 -2.1429 0.8571 -5.1429 2.1429 -3.0000 7.2857 2.1429 -5.1429 0.8571 -2.1429 -2.1429 2.1429 2.1429 0 0 -2.1429 0.8571 -5.1429 0 5.1429 -0.8571 0 -5.1429 0.8571 0 -0.8571 5.1429 0 2.1429 -2.1429 -2.1429 0 0 2.1429

k2 =

1.0e+06 *

5.1429 0 -5.1429 0.8571 0 -0.8571 0 2.1429 2.1429 -2.1429 -2.1429 0 -5.1429 2.1429 7.2857 -3.0000 -2.1429 0.8571 0.8571 -2.1429 -3.0000 7.2857 2.1429 -5.1429 0 -2.1429 -2.1429 2.1429 2.1429 0 -0.8571 0 0.8571 -5.1429 0 5.1429

k3 =

1.0e+06 *

2.1429 0 -2.1429 -2.1429 0 2.1429 0 5.1429 -0.8571 -5.1429 0.8571 0 -2.1429 -0.8571 7.2857

3.0000 -5.1429 -2.1429 -2.1429 -5.1429 3.0000 7.2857 -0.8571 -2.1429 0 0.8571 -5.1429 -0.8571 5.1429 0 2.1429 0 -2.1429 -2.1429 0 2.1429

k4 =

1.0e+06 *

2.1429 0 -2.1429 -2.1429 0 2.1429

0 5.1429 -0.8571 -5.1429 0.8571 0

-2.1429 -0.8571 7.2857 3.0000 -5.1429 -2.1429

-2.1429 -5.1429 3.0000 7.2857 -0.8571 -2.1429

0 0.8571 -5.1429 -0.8571 5.1429 0

2.1429 0 -2.1429 -2.1429 0 2.1429

调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵

求出的节点位移

U =

-0.0004

0.0008

0.0005

0.0010

0.0007

0.0023

-0.0007

0.0026

调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,Sxy

S1 =

1.0e+03 *

-4.4086

-0.7348

3.5914

S2 =

1.0e+03 *

4.4086

-0.6405

0.4086

S3 =

1.0e+03 *

1.8907

-1.0601

2.1093

S4 =

1.0e+03 *

-1.8907

2.1093

1.8907

二、

（1）如图划分四边形单元，工分成四个分别为

（2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N

调用 Quad2D4Node_Stiffness函数，求出单元刚度矩阵

调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵

求出节点位移

U =

0.0012

0.0017

-0.0012

0.0017

0.0016

0.0049

-0.0017

0.0052

调用Quad2D4Node_Stress函数，求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量

S1 =

1.0e+03 *

0.0000

-0.2478

2.0000

S2 =

1.0e+07 *

0.6856

4.1135

-1.7137

程序附录

一、

1、三角形单元总程序：

E=1e7;

NU=1/6;

t=1;

ID=1;

%调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵

k1=Triangle2D3Node_Stiffness(E,NU,t,0,1,0,0,1,1,ID)

k2=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,1,1,ID)

k3=Triangle2D3Node_Stiffness(E,NU,t,1,1,1,0,2,0,ID)

k4=Triangle2D3Node_Stiffness(E,NU,t,2,0,2,1,1,1,ID)

%调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵

KK = zeros(12,12);

KK=Triangle2D3Node_Assembly(KK,k1,1,2,3);

KK=Triangle2D3Node_Assembly(KK,k2,2,4,3);

KK=Triangle2D3Node_Assembly(KK,k3,3,4,5);

KK=Triangle2D3Node_Assembly(KK,k4,5,6,3)

% 边界条件的处理及刚度方程求解

k=KK(5:12,5:12)

p=[0;0;0;0;0;0;0;2000]

u=k\p

%支反力的计算

U=[0;0;0;0;u] %为节点位移

P=KK*U

%调用Triangle2D3Node_Strain函数，求出应变SN1、SN2、SN3中求出的分别为SNx,SNy,SNxy

u1=[U(1);U(2);U(3);U(4);U(5);U(6)];

u2=[U(3);U(4);U(7);U(8);U(5);U(6)];

u3=[U(5);U(6);U(7);U(8);U(9);U(10)];

u4=[U(9);U(10);U(11);U(12);U(5);U(6)];

SN1=Triangle2D3Node_Strain(0,1,0,0,1,1,u1)

SN2=Triangle2D3Node_Strain(0,0,1,0,1,1,u2)

SN3=Triangle2D3Node_Strain(1,1,1,0,2,0,u3)

SN4=Triangle2D3Node_Strain(2,0,2,1,1,1,u4)

%调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,Sxy

u1=[U(1);U(2);U(3);U(4);U(5);U(6)];

u2=[U(3);U(4);U(7);U(8);U(5);U(6)];

u3=[U(5);U(6);U(7);U(8);U(9);U(10)];

u4=[U(9);U(10);U(11);U(12);U(5);U(6)];

S1=Triangle2D3Node_Stress(E,NU,0,1,0,0,1,1,u1,ID)

S2=Triangle2D3Node_Stress(E,NU,0,0,1,0,1,1,u2,ID)

S3=Triangle2D3Node_Stress(E,NU,1,1,1,0,2,0,u3,ID)

S4=Triangle2D3Node_Stress(E,NU,2,0,2,1,1,1,u4,ID)

2、求刚度矩阵程序

function k=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)

%该函数计算单元的刚度矩阵

%输入弹性模量E，泊松比NU，厚度t

%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym

%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)

%输出单元刚度矩阵k(6X6)

%---------------------------------------------------------------

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;

3、求整体刚度矩阵

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

4、求应变程序

function strain=Triangle2D3Node_Strain(xi,yi,xj,yj,xm,ym,u)

%该函数计算单元的应变

%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym

%输入单元的位移列阵u(6X1)

%输出单元的应力strain(3X1)，由于它为常应变单元，则单元的应变分量为SNx,SNy,SNz

%---------------------------------------------------------------

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);

strain = B*u;

5、求应力程序

function stress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID) %该函数计算单元的应力

%输入弹性模量E，泊松比NU，厚度t

%输入三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym

%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)，单元的位移列阵u(6X1)

%输出单元的应力stress(3X1)，由于它为常应力单元，则单元的应力分量为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 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];

stress = D*B*u;

二、

1、四边形单元总程序：

E=1e7;

NU=1/6;

h=1;

ID=1;

%调用 Quad2D4Node_Stiffness函数，求出单元刚度矩阵

k1= Quad2D4Node_Stiffness(E,NU,h,0,1,0,0,1,0,1,1,ID)

k2= Quad2D4Node_Stiffness(E,NU,h,1,0,2,0,2,1,1,1,ID)

%调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵

KK=zeros(12,12);

KK= Quad2D4Node_Assembly(KK,k1,1,2,3,4);

KK= Quad2D4Node_Assembly(KK,k2,3,5,6,4)

% 边界条件的处理及刚度方程求解

k=KK(5:12,5:12)

p=[0;0;0;0;0;0;0;2000]

u=k\p

%支反力的计算

U=[0;0;0;0;u] %为节点位移

P=KK*U

%调用Quad2D4Node_Stress函数，求出单元应力中的的S1、S2、S3分别为Sx,Sy,Sxy应力分量

u1=[U(1);U(2);U(3);U(4);U(5);U(6);U(7);U(8)];

u2=[U(5);U(6);U(9);U(10);U(11);U(12);U(7);(8)];

S1= Quad2D4Node_Stress(E,NU,0,1,0,0,1,0,1,1,u1,ID)

S2= Quad2D4Node_Stress(E,NU,1,0,2,0,2,1,1,1,u2,ID)

2、求刚度矩阵程序

function k= Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID) %该函数计算单元的刚度矩阵

%输入弹性模量E，泊松比NU，厚度h

%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp

%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)

%输出单元刚度矩阵k(8X8)

%---------------------------------------------------------------syms s t;

a = (yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;

b = (yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;

c = (xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;

d = (xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;

B1 = [a*(t-1)/4-b*(s-1)/4 0 ; 0 c*(s-1)/4-d*(t-1)/4 ;

c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4];

B2 = [a*(1-t)/4-b*(-1-s)/4 0 ; 0 c*(-1-s)/4-d*(1-t)/4 ;

c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4];

B3 = [a*(t+1)/4-b*(s+1)/4 0 ; 0 c*(s+1)/4-d*(t+1)/4 ;

c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4];

B4 = [a*(-1-t)/4-b*(1-s)/4 0 ; 0 c*(1-s)/4-d*(-1-t)/4 ;

c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4];

Bfirst = [B1 B2 B3 B4];

Jfirst = [0 1-t t-s s-1 ; t-1 0 s+1 -s-t ;

s-t -s-1 0 t+1 ; 1-s s+t -t-1 0];

J = [xi xj xm xp]*Jfirst*[yi ; yj ; ym ; yp]/8;

B = Bfirst/J;

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

BD = J*transpose(B)*D*B;

r = int(int(BD, t, -1, 1), s, -1, 1);

z = h*r;

k = double(z);

3、求总体刚度矩阵程序

function z = Quad2D4Node_Assembly(KK,k,i,j,m,p)

%该函数进行单元刚度矩阵的组装

%输入单元刚度矩阵k，单元的节点编号i、j、m、p

%输出整体刚度矩阵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;

DOF(7)=2*p-1;

DOF(8)=2*p;

for n1=1:8

for n2=1:8

KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);

end

end

z=KK;

4、求应力程序

function stress= Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID) %该函数计算单元的应力

%输入弹性模量E，泊松比NU，厚度h，

%输入4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp,

%输入平面问题性质指示参数ID(1为平面应力，2为平面应变)

%输入单元的位移列阵u(8X1)

%输出单元的应力stress(3X1)

%由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy

%---------------------------------------------------------------syms s t;

a = (yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;

b = (yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;

c = (xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;

d = (xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;

B1 = [a*(t-1)/4-b*(s-1)/4 0 ; 0 c*(s-1)/4-d*(t-1)/4 ;

c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4];

B2 = [a*(1-t)/4-b*(-1-s)/4 0 ; 0 c*(-1-s)/4-d*(1-t)/4 ;

c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4];

B3 = [a*(t+1)/4-b*(s+1)/4 0 ; 0 c*(s+1)/4-d*(t+1)/4 ;

c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4];

B4 = [a*(-1-t)/4-b*(1-s)/4 0 ; 0 c*(1-s)/4-d*(-1-t)/4 ;

c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4];

Bfirst = [B1 B2 B3 B4];

Jfirst = [0 1-t t-s s-1 ; t-1 0 s+1 -s-t ;

s-t -s-1 0 t+1 ; 1-s s+t -t-1 0];

J = [xi xj xm xp]*Jfirst*[yi ; yj ; ym ; yp]/8;

B = Bfirst/J;

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

str1 = D*B*u;

str2 = subs(str1, {s,t}, {0,0});

stress = double(str2);

（注：可编辑下载，若有不当之处，请指正，谢谢!）

相关主题

相关文档

最新文档

© 2013-2022 www.sodocs.net  站点地图 | 侵权投诉

闽ICP备11023808号-8  本站资源均为网友上传分享，本站仅负责收集和整理，有任何问题请在对应网页下方投诉通道反馈