For my own interest I'm developing a finite element method (FEM) code for solving small problem. I already have working codes for 2D elements, trias and quadriliterals.
But I can't make any sense of the stiffness matrix for a hexahedron element. I certain I have the shape functions correct, and the derivates as well.
Below is a minimal working example in Matlab of my progress so far. I use Gauss integration.
It's only 1 elements, a 20-by-20-by-20 cube, with a point load of 100e3 N applied in one of the top corners, and fixed on the bottom. This corresponds to DOF's 1-12 being zero. The material is linear elastic with E=210e3 and a Poisson's ratio nu=0.3.
Things I could see being incorrect:
My logic for the element stiffness calculation is, in pseudo code:
for each gauss point:
Find values for gauss points xhi, eta, my
Find values for shape functions at current Gauss point
Find values for shape function derivates wrt xhi, eta and my at current Gauss point
Calculate Jacobian
Calculate shape functions derivates wrt x,y and z
Calculate strain displacement matrix B
Find element stiffness matrix Ke
clc
clear all
% 8 gauss points in total
xhi = [ -0.57735,-0.57735,-0.57735,-0.57735, 0.57735, 0.57735,0.57735,0.57735];
eta = [ -0.57735,-0.57735, 0.57735,0.57735,-0.57735,-0.57735,0.57735,0.57735];
my = [ -0.57735, 0.57735,-0.57735,0.57735,-0.57735, 0.57735,-0.57735,0.57735];
ngp=8; % 2 gauss points in each direction
coord = [0. 0. 0.;
20. 0. 0.;
20. 20. 0.;
0. 20. 0.;
0. 0. 20.;
20. 0. 20.;
20. 20. 20.;
0. 20. 20];
E=210e3;
poisson = 0.3;
Dcoeff = E/((1 + poisson) * (1 - 2 * poisson));
D = Dcoeff * [ (1 - poisson), poisson, poisson, 0, 0, 0;
poisson, (1 - poisson), poisson, 0, 0, 0;
poisson, poisson, (1 - poisson), 0, 0, 0;
0, 0, 0, ((1 - 2 * poisson)/2), 0, 0;
0, 0, 0, 0, ((1 - 2 * poisson)/2), 0;
0, 0, 0, 0, 0, ((1 - 2 * poisson)/2)];
zero = zeros(1,8);
K_cpp = zeros(24,24);
% if this was a real problem i would have to loop over each element
% but i only have 1 element, so only need to loop over the Gauss points
for i=1:8
% // shape functions, 1/8=0.125
N(1) = (1-xhi(i))*(1-eta(i))*(1-my(i))*0.125;
N(2) = (1+xhi(i))*(1-eta(i))*(1-my(i))*0.125;
N(3) = (1+xhi(i))*(1+eta(i))*(1-my(i))*0.125;
N(4) = (1-xhi(i))*(1+eta(i))*(1-my(i))*0.125;
N(5) = (1-xhi(i))*(1-eta(i))*(1+my(i))*0.125;
N(6) = (1+xhi(i))*(1-eta(i))*(1+my(i))*0.125;
N(7) = (1+xhi(i))*(1+eta(i))*(1+my(i))*0.125;
N(8) = (1-xhi(i))*(1+eta(i))*(1+my(i))*0.125;
% // derive shape functions wrt xhi
dNdXhi(1) = -(1-eta(i))*(1-my(i))*0.125;
dNdXhi(2) = (1-eta(i))*(1-my(i))*0.125;
dNdXhi(3) = (1+eta(i))*(1-my(i))*0.125;
dNdXhi(4) = -(1+eta(i))*(1-my(i))*0.125;
dNdXhi(5) = -(1-eta(i))*(1+my(i))*0.125;
dNdXhi(6) = (1-eta(i))*(1+my(i))*0.125;
dNdXhi(7) = (1+eta(i))*(1+my(i))*0.125;
dNdXhi(8) = -(1+eta(i))*(1+my(i))*0.125;
% // derive shape functions wrt eta
dNdEta(1) = -(1-xhi(i))*(1-my(i))*0.125;
dNdEta(2) = -(1+xhi(i))*(1-my(i))*0.125;
dNdEta(3) = (1+xhi(i))*(1-my(i))*0.125;
dNdEta(4) = (1-xhi(i))*(1-my(i))*0.125;
dNdEta(5) = -(1-xhi(i))*(1+my(i))*0.125;
dNdEta(6) = -(1+xhi(i))*(1+my(i))*0.125;
dNdEta(7) = (1+xhi(i))*(1+my(i))*0.125;
dNdEta(8) = (1-xhi(i))*(1+my(i))*0.125;
% // derive shape functions wrt my
dNdMy(1) = -(1-xhi(i))*(1-eta(i))*0.125;
dNdMy(2) = -(1+xhi(i))*(1-eta(i))*0.125;
dNdMy(3) = -(1+xhi(i))*(1+eta(i))*0.125;
dNdMy(4) = -(1-xhi(i))*(1+eta(i))*0.125;
dNdMy(5) = (1-xhi(i))*(1-eta(i))*0.125;
dNdMy(6) = (1+xhi(i))*(1-eta(i))*0.125;
dNdMy(7) = (1+xhi(i))*(1+eta(i))*0.125;
dNdMy(8) = (1-xhi(i))*(1+eta(i))*0.125;
% // find Jacobian for each Gauss point
dNdXhidEtadMy(1,:) = dNdXhi;
dNdXhidEtadMy(2,:) = dNdEta;
dNdXhidEtadMy(3,:) = dNdMy;
% J = [dNdXhi*coord(:,1),dNdXhi*coord(:,1),dNdXhi*coord(:,1);
J = dNdXhidEtadMy*coord;
% // find shape functions derivate matrix wrt x, y & z
dNdxdydz = inv(J)*dNdXhidEtadMy;
% // construct B matrix
q_x = dNdxdydz(1,:);
q_y = dNdxdydz(2,:);
q_z = dNdxdydz(3,:);
B_cpp = [ q_x, zero, zero;
zero, q_y, zero;
zero, zero, q_z;
q_y, q_x, zero;
zero, q_z, q_y;
q_z, zero, q_x];
% // compute Gauss point contribution to element Ke matrix
K_cpp = K_cpp + (B_cpp'*D*B_cpp*det(J));
end
% nod 1-4 dof 1-3 = 0:
% apply boundary conditions on nodes 0-4 --> dofs 1-12
for i=1:12
K_cpp(:,i)=zeros(1,24);
K_cpp(i,:)=zeros(24,1);
K_cpp(i,i)=1;
end
f=[zeros(1,23) 100e3]';
u=inv(K_cpp)*f;
This yields the solution u:
u =
1.0e+14*[0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.0293, -1.0531, -1.0511, -1.0503, -1.3835, -1.3835, -0.3304, -0.3304, -1.3791, -1.3895, -0.3304, -0.3304,
So the boundary conditions of prescibed zero displacement on the 12 first DOF's seems to work, but obviously the displacement is completely out of order. It should be something like 0.25mm, according to a commercial FE code I've double checked against.
I managed to solve the problem, with the help of the fine people at reddit.com/r/fea.
The issue was the strain-displacement matrix B.
Here's the incorrect B matrix, that I had:
[N1x N2x ... zeros(1,2*8)]
[zeros(1,2*8) N1y, N2y,... zeros(1,1*8)]
[B] = [zeros(1,2*8) N1z N2z ... ]
[N1y N1x ... N1x N2x ... zeros(1,1*8)]
[0 0 ... N1z N2z ... N1y N2y ...]
[N1z N2z ... zeros(1,1*8) N1x N2x ...]
Here is the correct defintion.
[N1x 0 0 N2x 0 0 ... ]
[0 N1y 0 0 N2y 0 ... ]
[B] = [0 0 N1z 0 0 N2z ... ]
[N1y N1x 0 N2y N2x 0 ... ]
[0 N1z N1y 0 N2z N2y ... ]
[N1z 0 N1x N2z 0 N2x ... ]