平面3节点三角形单元刚度矩阵matlab程序

%%%%%%%%%%%% Triangle2D3Node %%% begin %%%%%%%%%%%%%%%%%

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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%