八节点六面体等参数元计算程序

八节点六面体等参数元，单元刚度矩阵计算函数和整体矩阵拼装函数

其中：

Stiffnesske 为单元刚度计算函数；

BDcalc 为BD矩阵计算函数，被 Stiffnesske 函数所调用；

StiffnessAssemble 为整体矩阵拼装函数；

编写这几个函数的初衷是：最初使用了《有限元基础教程》(高等

教育出版社)所自带的计算函数，书中的计算函数采用符号变量编写，

计算速度较慢，且若是非正方体单元，计算时间非常非常长。因此编

写了这几个非符号变量形式的函数，计算速度极大的提高了。这几个

函数均已经过校验，计算结果与书中程序计算结果误差非常小。

自编各函数如下：

1.Stiffnesske函数：

-----------------------------------------------------------------------

function Ke=Stiffnesske(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为弹性模量和泊松比，x1-z8为单元八节点的坐标，输出Ke为单元刚度矩阵

%**************编写：杨玉磊 西安交通大学 机械学院 硕1005***************

Loc=[x1 y1 z1;x2 y2 z2;x3 y3 z3;x4 y4 z4;x5 y5 z5;x6 y6 z6;

x7 y7 z7;x8 y8 z8;];

gsx=[-0.7745966692 0 0.7745966692]; %高斯求积坐标点与系数

gsw=[0.55555555556 0.888888888889 0.55555555556];

Ke=zeros(24,24);

for ii=1:3 %三维高斯求积

sx=gsx(ii);

sw=gsw(ii);

for jj=1:3

nx=gsx(jj);

nw=gsw(jj);

for kk=1:3

tx=gsx(kk);

tw=gsw(kk);

Ke=Ke+sw*nw*tw*BDcalc(sx,nx,tx,Loc,E,NU);

end

end

end

-----------------------------------------------------------------------

2.BDcalc函数：

-----------------------------------------------------------------------

function BD=BDcalc(s,n,t,Loc,E,NU)

%根据给定的s,n,t数值计算函数值BD，用于外围的积分

%**************编写：杨玉磊 西安交通大学 机械学院 硕1005***************

N1=(1-s)*(1-n)*(1+t)/8;

N2=(1+s)*(1-n)*(1+t)/8;

N3=(1+s)*(1-n)*(1-t)/8;

N4=(1-s)*(1-n)*(1-t)/8;

N5=(1-s)*(1+n)*(1+t)/8;

N6=(1+s)*(1+n)*(1+t)/8;

N7=(1+s)*(1+n)*(1-t)/8;

N8=(1-s)*(1+n)*(1-t)/8;

dNsnt=[-(1-n)*(1+t)/8, -(1-s)*(1+t)/8, (1-s)*(1-n)/8; %N1-8对s,n,t的导数矩阵

(1-n)*(1+t)/8, -(1+s)*(1+t)/8, (1+s)*(1-n)/8;

(1-n)*(1-t)/8, -(1+s)*(1-t)/8, -(1+s)*(1-n)/8;

-(1-n)*(1-t)/8, -(1-s)*(1-t)/8, -(1-s)*(1-n)/8;

-(1+n)*(1+t)/8, (1-s)*(1+t)/8, (1-s)*(1+n)/8;

(1+n)*(1+t)/8, (1+s)*(1+t)/8, (1+s)*(1+n)/8;

(1+n)*(1-t)/8, (1+s)*(1-t)/8, -(1+s)*(1+n)/8;

-(1+n)*(1-t)/8, (1-s)*(1-t)/8, -(1-s)*(1+n)/8;];

dNsnt=dNsnt';

J=dNsnt*

Loc;

detJ=det(J);

dNxyz=J\dNsnt;

B=zeros(6,24);

for ii=1:8 %计算B矩阵

Bii=[dNxyz(1,ii) 0 0;0 dNxyz(2,ii) 0;0 0 dNxyz(3,ii);

dNxyz(2,ii) dNxy