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:
- Am I applying the boundary conditions in a in-correct way?
- Is the strain displacement matrix B incorrect?
- I don't think I need to shuffle around the stiffness matrix and assemble the global stiffness matrix as I only have 1 elements. Is this 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.
