八节点六面体单元有限元程序

admin

%%%%%%%%%%Hexahedral3D8Node%%%%%%%%%%%
function k=Hexahedral3D8Node_Stiffness(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8)
%用于计算单元的刚度矩阵阵
%输入量为弹性模量E，泊松比NU，输入的8个节点的坐标x1y1,z1,....x8,y8,z8
%输出量为单元刚度矩阵k（24*24）
syms s t n;%定义局部坐标系
%定义型函数N
N1=(1+s)*(1-t)*(1-n)/8;
N2=(1+s)*(1+t)*(1-n)/8;
N3=(1-s)*(1+t)*(1-n)/8;
N4=(1-s)*(1-t)*(1-n)/8;
N5=(1+s)*(1-t)*(1+n)/8;
N6=(1+s)*(1+t)*(1-n)/8;
N7=(1-s)*(1+t)*(1+n)/8;
N8=(1-s)*(1-t)*(1+n)/8;
%%定义坐标变换
x=N1*x1+N2*x2+N3*x3+N4*x4+N5*x5+N6*x6+N7*x7+N8*x8;
y=N1*y1+N2*y2+N3*y3+N4*y4+N5*y5+N6*y6+N7*y7+N8*y8;
z=N1*z1+N2*z2+N3*z3+N4*z4+N5*z5+N6*z6+N7*z7+N8*z8;
%定义雅克比矩阵
J=[diff(x,s),diff(y,s),diff(z,s);diff(x,t),diff(y,t),diff(z,t);diff(x,n),diff(y,n),diff(z,n)];
Jdet=det(J);
%%%定义B矩阵的系数
a=diff(y,t)*diff(z,n)-diff(z,t)*diff(y,n);
b=diff(y,s)*diff(z,n)-diff(z,s)*diff(y,n);
c=diff(y,s)*diff(z,t)-diff(z,s)*diff(y,t);
d=diff(x,t)*diff(z,n)-diff(z,t)*diff(x,n);
e=diff(x,s)*diff(z,n)-diff(z,s)*diff(x,n);
f=diff(x,s)*diff(z,t)-diff(z,s)*diff(x,n);
g=diff(x,t)*diff(y,n)-diff(y,t)*diff(x,n);
h=diff(x,s)*diff(y,n)-diff(y,s)*diff(x,n);
l=diff(x,s)*diff(y,t)-diff(y,s)*diff(x,t);
%通过循环计算各个矩阵
Ns=[N1,N2,N3,N4,N5,N6,N7,N8];
Bs=sym(zeros(6,3,8));
for i=1:8
    Bs(:,:,i)=[a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n),0,0；
        0,-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(i),n),0;0,0,g*diff(Ns(i),s)-h*diff(Ns(i),t)+l*diff(Ns(i),n);
        0,g*diff(Ns(i),s)-h*diff(Ns(i),t)+l*diff(Ns(i),n),-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(i),n);
        g*diff(Ns(i),s)-h*diff(Ns(i),t)+l*diff(Ns(i),n),0,a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n)]/Jdet;
end
%%计算B矩阵
B=[Bs(:,:,1),Bs(:,:,2),Bs(:,:,3),Bs(:,:,4),Bs(:,:,5),Bs(:,:,6),Bs(:,:,7),Bs(:,:,8)];
B
%%计算D矩阵
D=(E/((1+NU)*(1-2*NU)))*[1-NU,NU,NU,0,0,0;NU,1-NU,NU,0,0,0;NU,NU,1-NU,0,0,0;0,0,0,0.5-NU,0,0;0,0,0,0,0.5-NU,0;
    0,0,0,0,0,0.5-NU];
D
%%计算单元刚度矩阵
BD=Jdet*transpose(B)*D*B;
z=int(int(int(BD,n,-1,1),t,-1,1)s,-1,1);
z
k=double(z);
function KK=Hexahedral3D8Node_Assemble(KK,k,i,j,l,m,n,o,p,q)
%%该函数进行单元刚度矩阵的组装
%输入量为单元刚度矩阵k，节点编号i,j,l,m,n,o,p,q
%输出量为整体的刚度矩阵KK
DOF=[3*i-2:3*i,3*j-2:3*j,3*l-2:3*l,3*m-2:3*m,3*n-2:3*n,3*o-2:3*o,3*p-2:3*p,3*q-2:3*q];
for n1=1:24
    for n2=1:24
        KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2)
    end
end
function stress=Hexahedral3D8Node_Stress(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,u)
%%该函数计算单元中心点的应力
%输入量为弹性模量E，泊松比NU，输入八个节点的坐标x1,y1,z1,......x2,y2,z2
%输入为单元位移矩阵u(6*1)
%输出为中心点的应力stress（6*1），Sx,Sy,Sz,Sxy,Syz,Szx
syms s t n;
%%定义型函数
N1=(1+s)*(1-t)*(1-n)/8;
N2=(1+s)*(1+t)*(1-m)/8;
N3=(1-s)*(1+t)*(1-n)/8;
N4=(1-s)*(1-t)*(1-n)/8;
N5=(1+s)*(1-t)*(1+n)/8;
N6=(1+s)*(1+t)*(1-n)/8;
N7=(1-s)*(1+t)*(1+n)/8;
N8=(1-s)*(1-t)*(1+n)/8;
%%定义坐标变换
x=N1*x1+N2*x2+N3*x3+N4*x4+N5*x5+N6*x6+N7*x7+N8*x8;
y=N1*y1+N2*y2+N3*y3+N4*y4+N5*y5+N6*y6+N7*y7+N8*y8;
z=N1*z1+N2*z2+N3*z3+N4*z4+N5*z5+N6*z6+N7*z7+N8*z8;
%定义雅克比矩阵
J=[diff(x,s),diff(y,s),diff(z,s);diff(x,t),diff(y,t),diff(z,t);diff(x,m),diff(y,n),diff(z,n)];
Jdet=det(J);
%%%定义B矩阵的系数
a=diff(y,t)*diff(z,n)-diff(z,t)*diff(y,n);
b=diff(y,s)*diff(z,n)-diff(z,s)*diff(y,n);
c=diff(y,s)*diff(z,t)-diff(z,s)*diff(y,t);
d=diff(x,t)*diff(z,n)-diff(z,t)*diff(x,n);
e=diff(x,s)*diff(z,n)-diff(z,s)*diff(x,n);
f=diff(x,s)*diff(z,t)-diff(z,s)*diff(x,f);
g=diff(x,t)*diff(y,n)-diff(y,t)*diff(x,n);
h=diff(x,s)*diff(y,n)-diff(y,s)*diff(x,n);
l=diff(x,s)*diff(y,t)-diff(y,s)*diff(x,t);
%通过循环计算各个矩阵
Ns=[N1,N2,N3,N4,N5,N6,N7,N8];
Bs=sym(zeros(6,3,8));
for i=1:8
    Bs(:,:,i)=[a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n),0,;
        0,-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(i),n),-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(i),n);
        g*diff(Ns(i),s)-h*diff(Ns(i),t)+l*diff(Ns(i),n),0,a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n)]/Jdet;
end
%%计算B矩阵
B=[Bs(:,:,1),Bs(:,:,2),Bs(:,:,3),Bs(:,:,4),Bs(:,:,5),Bs(:,:,6),Bs(:,:,7),Bs(:,:,8)];
B
%计算D矩阵
D=(E/((1+NU)*(1-2*NU)))*[1-NU,NU,NU,0,0,0;NU,1-NU,NU,0,0,0;NU,NU,1-NU,0,0,0;0,0,0,0.5-NU,0,0;0,0,0,0,0.5-NU,0;
    0,0,0,0,0,0.5-NU];
D
%%计算应力向量
w=D*B*u;
%在单元的中心处计算应力
wcent=subs(w,{s,t,n},{0,0,0})
stress=double(wcent)