有限元课程设计

相关主题

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

一．问题描述

如图所示的平面矩形结构，设E=1，NU=0.25，h=1，考虑以下约束和外载：

位移边界条件BC（u）：U A=0，V A=0，U D=0，

力边界条件BC（p）：在CD边上有均布载荷q=1，

建模情形：使用四个四节点矩形单元，

试在该建模情形下，求各节点的位移以及各个单元的应力分布。

二．Matlab程序

(1).函数定义：

function k= Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID) 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 xjxmxp]*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);

end

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

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

z=KK;

end

function stress= Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID) syms s t;