8节点六面体单元刚度矩阵的计算

xinyuxf

找了很多资料也没找到8节点六面体单元的单元刚度矩阵的求法，无奈自己编了一个，分享一下，也希望大家多多指教其中的不足和错误。
%%%%X、Y节点坐标；E、UM弹性模量和泊松比；lmd单元连接的各节点坐标
function EK = RectESM(X,Y,Z,E,NU,lmd)
x1 = X(lmd(1)); x2 = X(lmd(2)); x3 = X(lmd(3)); x4 = X(lmd(4));   %%%%8个节点的坐标
y1 = Y(lmd(1)); y2 = Y(lmd(2)); y3 = Y(lmd(3)); y4 = Y(lmd(4));
z1 = Z(lmd(1)); z2 = Z(lmd(2)); z3 = Z(lmd(3)); z4 = Z(lmd(4));
x5 = X(lmd(5)); x6 = X(lmd(6)); x7 = X(lmd(7)); x8 = X(lmd(8));   
y5 = Y(lmd(5)); y6 = Y(lmd(6)); y7 = Y(lmd(7)); y8 = Y(lmd(8));
z5 = Z(lmd(5)); z6 = Z(lmd(6)); z7 = Z(lmd(7)); z8 = Z(lmd(8));
D=[1-NU  NU  NU  0  0  0
       NU  1-NU NU  0  0  0      
       NU   NU 1-NU 0  0  0
       0 0 0 (1-2*NU)/2 0 0
       0 0 0 0 (1-2*NU)/2 0
      0 0 0 0 0 (1-2*NU)/2];
D=E/(1+NU)/(1-2*NU)*D;
A=[-0.774597 0 0.774597];            
H=[0.555556 0.888889 0.555556];       %%%            下面采用3节点高斯积分法求解三次积分         %%%
z=zeros(24);                        
for p=1:3                           
    r=A(p);                           
    for m=1:3
        s=A(m);
        for n=1:3
            t=A(n);
            a=((s-1)*(1-t)*x1+(1-s)*(1-t)*x2+(1+s)*(1-t)*x3+(1+s)*(t-1)*x4+(s-1)*(1+t)*x5+(1-s)*(1+t)*x6+(1+s)*(1+t)*x7-(1+s)*(1+t)*x8)/8;
            b=((r-1)*(1-t)*x1+(1+r)*(t-1)*x2+(1+r)*(1-t)*x3+(1-r)*(1-t)*x4+(r-1)*(1+t)*x5-(1+r)*(1+t)*x6+(1+r)*(1+t)*x7+(1-r)*(1+t)*x8)/8;
            c=((r-1)*(1-s)*x1+(1+r)*(s-1)*x2-(1+r)*(1+s)*x3+(r-1)*(1+s)*x4+(1-r)*(1-s)*x5+(1+r)*(1-s)*x6+(1+r)*(1+s)*x7+(1-r)*(1+s)*x8)/8;
            d=((s-1)*(1-t)*y1+(1-s)*(1-t)*y2+(1+s)*(1-t)*y3+(1+s)*(t-1)*y4+(s-1)*(1+t)*y5+(1-s)*(1+t)*y6+(1+s)*(1+t)*y7-(1+s)*(1+t)*y8)/8;
            e=((r-1)*(1-t)*y1+(1+r)*(t-1)*y2+(1+r)*(1-t)*y3+(1-r)*(1-t)*y4+(r-1)*(1+t)*y5-(1+r)*(1+t)*y6+(1+r)*(1+t)*y7+(1-r)*(1+t)*y8)/8;
            f=((r-1)*(1-s)*y1+(1+r)*(s-1)*y2-(1+r)*(1+s)*y3+(r-1)*(1+s)*y4+(1-r)*(1-s)*y5+(1+r)*(1-s)*y6+(1+r)*(1+s)*y7+(1-r)*(1+s)*y8)/8;
            g=((s-1)*(1-t)*z1+(1-s)*(1-t)*z2+(1+s)*(1-t)*z3+(1+s)*(t-1)*z4+(s-1)*(1+t)*z5+(1-s)*(1+t)*z6+(1+s)*(1+t)*z7-(1+s)*(1+t)*z8)/8;
            h=((r-1)*(1-t)*z1+(1+r)*(t-1)*z2+(1+r)*(1-t)*z3+(1-r)*(1-t)*z4+(r-1)*(1+t)*z5-(1+r)*(1+t)*z6+(1+r)*(1+t)*z7+(1-r)*(1+t)*z8)/8;
            i=((r-1)*(1-s)*z1+(1+r)*(s-1)*z2-(1+r)*(1+s)*z3+(r-1)*(1+s)*z4+(1-r)*(1-s)*z5+(1+r)*(1-s)*z6+(1+r)*(1+s)*z7+(1-r)*(1+s)*z8)/8;
            Jfirst=[a d g;b e h;c f i];
            J=det(Jfirst);
            Nr(1)=(s-1)*(1-t)/8; Ns(1)=(r-1)*(1-t)/8;  Nt(1)=(r-1)*(1-s)/8; Nr(2)=(1-s)*(1-t)/8; Ns(2)=(1+r)*(t-1)/8;  Nt(2)=(1+r)*(s-1)/8;
            Nr(3)=(1+s)*(1-t)/8; Ns(3)=(1+r)*(1-t)/8;  Nt(3)=(-1-r)*(1+s)/8; Nr(4)=(1+s)*(t-1)/8; Ns(4)=(1-r)*(1-t)/8;  Nt(4)=(r-1)*(1+s)/8;
            Nr(5)=(s-1)*(1+t)/8; Ns(5)=(r-1)*(1+t)/8;  Nt(5)=(1-r)*(1-s)/8; Nr(6)=(1-s)*(1+t)/8; Ns(6)=(-1-r)*(1+t)/8;  Nt(6)=(1+r)*(1-s)/8;
            Nr(7)=(1+s)*(1+t)/8; Ns(7)=(1+r)*(1+t)/8;  Nt(7)=(1+r)*(1+s)/8; Nr(8)=(-1-s)*(1+t)/8; Ns(8)=(1-r)*(1+t)/8;  Nt(8)=(1-r)*(1+s)/8;
            for k=1:1:8
               Nrst=[Nr(k);Ns(k);Nt(k)];
               Nxyz=inv(Jfirst)*Nrst;
               B(:,k*3-2:k*3)=[Nxyz(1) 0 0; 0 Nxyz(2) 0;0 0 Nxyz(3);Nxyz(2) Nxyz(1) 0;0  Nxyz(3) Nxyz(2);Nxyz(3) 0 Nxyz(1)];
            end
            BD=transpose(B)*D*B*J;
            z=z+H(p)*H(m)*H(n)*BD;
        end
    end
end   
EK=double(z);
return;

请大家指指错，呵

[ 本帖最后由 xinyuxf 于 2006-9-14 10:37 编辑 ]

fx2003521

先看看先

xinyuxf

恩，以后再做个实例，验证一下，肯定有不足的地方，而且程序应该还能简化

skyfing12

好东西,强烈顶上!!
老兄能否仔细给兄弟介绍一下有限元计算的步骤!
小弟初出江湖,不知如何下手!!
skyfing12@yahoo.com.cn

mxlzhenzhu

兄弟，我也是用你种Gauss积分方案，但是，你觉得精度如何？我的试验结果表明，不好啊。

我试过2点，3点一直到6点的Gauss积分方案，结果都差不多，和有限元计算结果对比，非常不好啊，这问题困扰我很久了。

我对比过，用标准的六面体，长宽高为2，用有限元计算的和自己计算的不一样。
我去下载了一个matlab femtoolbox，发现他的计算结果和我的一样。

我用四面体计算结果和有限元计算结果，一模一样，那没的说。

排除一切的一切原因，我认为通用有限元里面不是用Gauss积分；不知道我又是不是在信口开河了？求指导。

mxlzhenzhu

什么是Hourglassing,控制原理是什么？