Calculation stiffness matrix for 8-node hexahedral element.

Hi, I'm trying to solve stifness matrix for standard hexaheral element (100x100x100), I'm doing everything like in a literature but the final result from my solver does not match the result from the program (Midas NFX).The compared value is the nodal forces Qe . Can you help me find the bug?

clear

syms   r s t % r,s,t - local coordinates
E=210000;% Modulus
v=0.28;%Poisson
G=E/(2*(1+v));%Shear modulus
coor_glob=[0 0 100
          100 0 100
          0 100 100
          100 100 100
          0 0 0
          100 0 0 
          0 100 0 
          100 100 0];% global coordinates

x1= coor_glob(1,1); y1= coor_glob(1,2); z1= coor_glob(1,3); 
x2= coor_glob(2,1); y2= coor_glob(2,2); z2= coor_glob(2,3);
x3= coor_glob(3,1); y3= coor_glob(3,2); z3= coor_glob(3,3); 
x4= coor_glob(4,1); y4= coor_glob(4,2); z4= coor_glob(4,3);  
x5= coor_glob(5,1); y5= coor_glob(5,2); z5= coor_glob(5,3); 
x6= coor_glob(6,1); y6= coor_glob(6,2); z6= coor_glob(6,3); 
x7= coor_glob(7,1); y7= coor_glob(7,2); z7= coor_glob(7,3); 
x8= coor_glob(8,1); y8= coor_glob(8,2); z8= coor_glob(8,3);

S=[1/E -v/E -v/E 0 0 0;
  -v/E 1/E -v/E 0 0 0;
  -v/E -v/E 1/E 0 0 0;
  0 0 0 1/(G) 0 0;
  0 0 0 0 1/(G) 0 ;
  0 0 0 0 0 1/(G)];
D=inv(S); %elasticity matrix 

N1=1/8*((1-r)*(1-s)*(1+t));
N2=1/8*((1+r)*(1-s)*(1+t));
N3=1/8*((1-r)*(1+s)*(1+t));
N4=1/8*((1+r)*(1+s)*(1+t));
N5=1/8*((1-r)*(1-s)*(1-t));
N6=1/8*((1+r)*(1-s)*(1-t));
N7=1/8*((1-r)*(1+s)*(1-t));
N8=1/8*((1+r)*(1+s)*(1-t));% shape functions 

x=N1*x1+N2*x2+N3*x3+N4*x4+N5*x5+N6*x6+N7*x7+N8*x8; %interpolate coord.
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,r) diff(y,r) diff(z,r);
    diff(x,s) diff(y,s) diff(z,s);
     diff(x,t) diff(y,t) diff(z,t)];%Jacobian matrix
 
d_N1=inv(J)*[diff(N1,r);diff(N1,s);diff(N1,t)];
d_N2=inv(J)*[diff(N2,r);diff(N2,s);diff(N2,t)];
d_N3=inv(J)*[diff(N3,r);diff(N3,s);diff(N3,t)];
d_N4=inv(J)*[diff(N4,r);diff(N4,s);diff(N4,t)];
d_N5=inv(J)*[diff(N5,r);diff(N5,s);diff(N5,t)];
d_N6=inv(J)*[diff(N6,r);diff(N6,s);diff(N6,t)];
d_N7=inv(J)*[diff(N7,r);diff(N7,s);diff(N7,t)];
d_N8=inv(J)*[diff(N8,r);diff(N8,s);diff(N8,t)];%differentials 

B1=[d_N1(1)	    0	     0
    0	     d_N1(2)	 0
    0	       0	  d_N1(3)
    d_N1(2)	 d_N1(1)	 0
    0	     d_N1(3)  d_N1(2)
    d_N1(3)	   0	  d_N1(1)];%strain matrix for node 1
B2=[d_N2(1)	    0	     0
    0	     d_N2(2)	 0
    0	       0	  d_N2(3)
    d_N2(2)	 d_N2(1)	 0
    0	     d_N2(3)  d_N2(2)
    d_N2(3)	   0	  d_N2(1)];%strain matrix for node 2
B3=[d_N3(1)	    0	     0
    0	     d_N3(2)	 0
    0	       0	  d_N3(3)
    d_N3(2)	 d_N3(1)	 0
    0	     d_N3(3)  d_N3(2)
    d_N3(3)	   0	  d_N3(1)];%strain matrix for node 3
B4=[d_N4(1)	    0	     0
    0	     d_N4(2)	 0
    0	       0	  d_N4(3)
    d_N4(2)	 d_N4(1)	 0
    0	     d_N4(3)  d_N4(2)
    d_N4(3)	   0	  d_N4(1)];%strain matrix for node 4
B5=[d_N5(1)	    0	     0
    0	     d_N5(2)	 0
    0	       0	  d_N5(3)
    d_N5(2)	 d_N5(1)	 0
    0	     d_N5(3)  d_N5(2)
    d_N5(3)	   0	  d_N5(1)];%strain matrix for node 5
B6=[d_N6(1)	    0	     0
    0	     d_N6(2)	 0
    0	       0	  d_N6(3)
    d_N6(2)	 d_N6(1)	 0
    0	     d_N6(3)  d_N6(2)
    d_N6(3)	   0	  d_N6(1)];%strain matrix for node 6
B7=[d_N7(1)	    0	     0
    0	     d_N7(2)	 0
    0	       0	  d_N7(3)
    d_N7(2)	 d_N7(1)	 0
    0	     d_N7(3)  d_N7(2)
    d_N7(3)	   0	  d_N7(1)];%strain matrix for node 7
B8=[d_N8(1)	    0	     0
    0	     d_N8(2)	 0
    0	       0	  d_N8(3)
    d_N8(2)	 d_N8(1)	 0
    0	     d_N8(3)  d_N8(2)
    d_N8(3)	   0	  d_N8(1)];%strain matrix for node 8

B=[B1 B2 B3 B4 B5 B6 B7 B8];%main strin matrix for hexahedral

f(r,s,t)=B'*D*B*det(J);% function to integration

N=[-sqrt(1/3) sqrt(1/3)];% Gauss-Legendre quadrature points fot N=2
w=[1 1];% Gauss-Legendre quadrature weights fot N=2
wij=[w(1)*w(1)*w(2);
     w(2)*w(1)*w(2);
     w(1)*w(2)*w(2);
     w(2)*w(2)*w(2);
     w(1)*w(2)*w(2);
     w(2)*w(1)*w(1);
     w(1)*w(2)*w(1);
     w(2)*w(2)*w(1)];%matrix of weights for future multiplikation

%% Gauss-Legendre integration:
k=f(N(1),N(1),N(2))*wij(1)+f(N(2),N(1),N(2))*wij(2)+f(N(1),N(2),N(2))*wij(3)+f(N(2),N(2),N(2))*wij(4)+f(N(1),N(1),N(1))*wij(5)+f(N(2),N(1),N(1))*wij(6)+f(N(1),N(2),N(1))*wij(7)+f(N(2),N(2),N(1))*wij(8);
     
K=eval(k); %elemennt stiffnes matrix
u=[0 0 0 0 0 0 0 0 0 -2.9665e-003 7.7690e-004 7.7690e-004 0 0 0 0 0 0 0 0 0 0 0 0]'; %translations from FEM program
Qe=K*u %nodal forces 

%% Nodal forces from FEM program
Qe_fem=[4279.5049 2764.0085 -241.1585 -2782.678 -869.3088 -1551.515 4012.6929 1291.8174 1291.8174 -10000 0 0 3283.1487 675.8129 675.8129 -289.4951 -2069.6565 -2069.6565 4279.5049 -241.1585 2764.0085 -2782.678 -1551.515 -869.3088]';

Answer 1

[Boris Kotsyubuk]

It seems like your 4th row in your strain matrices (B1, B2, etc.) is incorrect. I believe the 4th row should be your last row in the matrix, while the 5th and 6th rows move up one row. I think the strain matrix for B1 should look like this:

B1 = [ d_N1(1), 0, 0;

0, d_N1(2), 0;

0, 0, d_N1(3);

0, d_N1(3), d_N1(2);

d_N1(3), 0, d_N1(1);

d_N1(2), d_N1(1), 0 ];

Hope this helps!

[Alehegn Tesfaye]

I did B matrix on octave exactly as you did but it still won't work, any help.

% Coordinates parameters

N=32;

x1=0; y1=2; z1=0;

x2=0; y2=0; z2=0;

x3=2; y3=0; z3=0;

x4=2; y4=2; z4=0;

x5=0; y5=2; z5=2;

x6=0; y6=0; z6=2;

x7=2; y7=0; z7=2;

x8=2; y8=2; z8=2;

% Element parameters

E=(170+N);

u=0.3;

% Plane Stress or Strain parameters

D=(E/((1+u)+(1-(2*u))))*[1-u u u 0 0 0;u 1-u u 0 0 0;u u 1-u 0 0 0;0 0 0 (1-(2*u))/2 0 0;0 0 0 0 (1-(2*u))/2 0; 0 0 0 0 0 (1-(2*u))/2];

% Load parameters

% Shape Function

syms s t w real;

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

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

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

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

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

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

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

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

% Coordinate Mapping

x=simplify(N1*x1+N2*x2+N3*x3+N4*x4+N5*x5+N6*x6+N7*x7+N8*x8);

y=simplify(N1*y1+N2*y2+N3*y3+N4*y4+N5*y5+N6*y6+N7*y7+N8*y8);

z=simplify(N1*z1+N2*z2+N3*z3+N4*z4+N5*z5+N6*z6+N7*z7+N8*z8);

% N Matrix

% Jacobian Matrix

syms J J0 real;

dxds=diff(x,s);

dxdt=diff(x,t);

dxdw=diff(x,w);

dyds=diff(y,s);

dydt=diff(y,t);

dydw=diff(y,w);

dzds=diff(z,s);

dzdt=diff(z,t);

dzdw=diff(z,w);

J=[dxds dyds dzds;dxdt dydt dzdt;dxdw dydw dzdw];

Jdet=simplify(det(J));

Jinv=simplify(inv(J));

% B Matrix

syms B real;

d_N1=Jinv*[diff(N1,s);diff(N1,t);diff(N1,w)];

d_N2=Jinv*[diff(N2,s);diff(N2,t);diff(N2,w)];

d_N3=Jinv*[diff(N3,s);diff(N3,t);diff(N3,w)];

d_N4=Jinv*[diff(N4,s);diff(N4,t);diff(N4,w)];

d_N5=Jinv*[diff(N5,s);diff(N5,t);diff(N5,w)];

d_N6=Jinv*[diff(N6,s);diff(N6,t);diff(N6,w)];

d_N7=Jinv*[diff(N7,s);diff(N7,t);diff(N7,w)];

d_N8=Jinv*[diff(N8,s);diff(N8,t);diff(N8,w)];

B1=[d_N1(1) 0 0

0 d_N1(2) 0

0 0 d_N1(3)

0 d_N1(3) d_N1(2)

d_N1(3) 0 d_N1(1)

d_N1(2) d_N1(1) 0];

B2=[d_N2(1) 0 0

0 d_N2(2) 0

0 0 d_N2(3)

0 d_N2(3) d_N2(2)

d_N2(3) 0 d_N2(1)

d_N2(2) d_N2(1) 0];

B3=[d_N3(1) 0 0

0 d_N3(2) 0

0 0 d_N3(3)

0 d_N3(3) d_N3(2)

d_N3(3) 0 d_N3(1)

d_N3(2) d_N3(1) 0];

B4=[d_N4(1) 0 0

0 d_N4(2) 0

0 0 d_N4(3)

0 d_N4(3) d_N4(2)

d_N4(3) 0 d_N4(1)

d_N4(2) d_N4(1) 0];

B5=[d_N5(1) 0 0

0 d_N5(2) 0

0 0 d_N5(3)

0 d_N5(3) d_N5(2)

d_N5(3) 0 d_N5(1)

d_N5(2) d_N5(1) 0];

B6=[d_N6(1) 0 0

0 d_N6(2) 0

0 0 d_N6(3)

0 d_N6(3) d_N6(2)

d_N6(3) 0 d_N6(1)

d_N6(2) d_N6(1) 0];

B7=[d_N7(1) 0 0

0 d_N7(2) 0

0 0 d_N7(3)

0 d_N7(3) d_N7(2)

d_N7(3) 0 d_N7(1)

d_N7(2) d_N7(1) 0];

B8=[d_N8(1) 0 0

0 d_N8(2) 0

0 0 d_N8(3)

0 d_N8(3) d_N8(2)

d_N8(3) 0 d_N8(1)

d_N8(2) d_N8(1) 0];

B=simplify([B1 B2 B3 B4 B5 B6 B7 B8]);

BT=B';

% Integration function

syms F real;

F=BT*D*B*Jdet;

% Stiffness Matrix

p=(1/sqrt(3));

o=function_handle(F,'vars',[s,t,w]);

K=(o(-p,-p,-p)+o(p,-p,-p)+o(p,p,-p)+o(-p,p,-p)+o(-p,-p,p)+o(p,-p,p)+o(p,p,p)+o(-p,p,p))