hex/brick/六面体的单元刚度矩阵

%-----------------------------------------------------------
function k=Hexahedral3D8Node_Stiffness（ENUx1y1z1x2y2z2x3y3z3x4y4z4x5y5z5x6y6z6x7y7z7x8y8z8）
%该函数计算单元的刚度矩阵
%输入弹性模量E，泊松比NU
%输入8个节点的坐标x1、y1、z1、x2、y2、z2、x3、y3、z3、x4、y4、z4、
%输入x5、y5、z5、x6、y6、z6、x7、y7、z7、x8、y8、z8
%输出单元刚度矩阵k（24X24） 
%---------------------------------------------------------------
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（xs）diff（ys）diff（zs）;diff（xt）diff（yt）diff（zt）;diff（xn）diff（yn）diff（zn）];
Jdet=det（J）;
J;
Jdet;
%   定义 B 矩阵的系数
a=diff（yt）*diff（zn）-diff（zt）*diff（yn）;
b=diff（ys）*diff（zn）-diff（zs）*diff（yn）;
c=diff（ys）*diff（zt）-diff（zs）*diff（yt）;
d=diff（xt）*diff（zn）-diff（zt）*diff（xn）;
e=diff（xs）*diff（zn）-diff（zs）*diff（xn