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

xinyuxf · 发表于 2006-9-14 08:13

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

找了很多资料也没找到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 · 发表于 2006-9-14 14:47

先看看先

xinyuxf · 发表于 2006-9-14 16:10

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

skyfing12 · 发表于 2006-12-19 22:52

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

mxlzhenzhu · 发表于 2013-8-13 16:57

本帖最后由 mxlzhenzhu 于 2013-8-13 23:41 编辑

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

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

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

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

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

mxlzhenzhu · 发表于 2013-8-13 17:37

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

账号		自动登录	找回密码
密码			我要加入

[计算力学] 8节点六面体单元刚度矩阵的计算

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分