Uniaxial stretching solution not uniform in FEM code

Question

I am trapped here for a long time. I wrote a toy Matlab FEM code. I want to run the follow simulation.
Mesh
Suppose we have a cube, and we divide it into subcube along $x,y,z$ axis, then each subcube be divide into 5 tetrahedral elements. The cube along $x,y,z$ axis is divide into $n_x,n_y,n_z$ node, respectively. So there are $n_x times n_y times n_z$ nodes in the grid. Each node in grid can be index with $node(i,j,k)$
Boundary conditions

The cube's $xoy$ plane has fixed the $z$ DOF, cube's $xoz$ plane has fixed $y$ DOF, and cube's $yoz$ plane has fixed $x$ DOF.

The cube's top plane has a uniform pressure applied (tensile).

So the cube is under uniaxial (tensile) stress.
Problems
I think under uniaxial stretching, the displacement along the z-axis of nodes having the same $z$ coordinate should be equal. That is $u_z(x, y, const) = const$.
However, my code doesn't give me this solution. I do not know it is my code's bug or it is the truth at they are not equal.
Main code
Actually, I already check the code with examples from some books. I already checked the element matrix and matrix assemble. My code can obtain the same solution with a standard book. Maybe there is something wrong with my loading. But I can not find a clue.
The full code and checking can be found in this repository
clc
clear
close all

num_x_node=4;
num_y_node=4;
num_z_node=4;

x_range=10;
y_range=10;
z_range=10;

grid= CubeDomainTetGrid(num_x_node, num_y_node, num_z_node, x_range, y_range, z_range);
node_coordinate_table=grid.nodeCoordinateTable();
element_node_table=grid.elementNodeTable();

dof_per_node=3;
num_node = grid.numNode();
num_dof = num_node * dof_per_node;

is_constrain = zeros(num_dof,1);

load_value = zeros(num_dof,1);
is_load = zeros(num_dof,1);

% assemble global stiffness matrix
element_stiffness_matrix_array = calculateElementStiffnessMatrixArray(node_coordinate_table, element_node_table);
K = assembleGlobalMatrix(node_coordinate_table, element_node_table, dof_per_node, element_stiffness_matrix_array);

%% constrain
% all node at x-y plane is fix z dof to zero
for i=1:1:num_x_node
    for j=1:1:num_y_node
        for k=1
            node_index = grid.nodeIndex(i,j,k);
            for dof_index = 3
                dof_global_index = (node_index - 1) * dof_per_node + dof_index;
                is_constrain(dof_global_index)=1;
            end
        end
    end
end
% all node at y-z plane is fix x dof to zero
for i=1:1:1
    for j=1:1:num_y_node
        for k=1:1:num_z_node
            node_index = grid.nodeIndex(i,j,k);
            for dof_index = 1
                dof_global_index = (node_index - 1) * dof_per_node + dof_index;
                is_constrain(dof_global_index)=1;
            end
        end
    end
end
% all node at x-z plane is fix y dof to zero
for i=1:1:num_x_node
    for j=1
        for k=1:1:num_z_node
            node_index = grid.nodeIndex(i,j,k);
            for dof_index = 2
                dof_global_index = (node_index - 1) * dof_per_node + dof_index;
                is_constrain(dof_global_index)=1;
            end
        end
    end
end
clear node_index;
clear dof_index;
clear dof_global_index;

%% load
% node(:,:,num_z_node) is apply pressure
pressure = 6;

% tranvers cell on surface and accumulate pressure
% the face see from top like
% the node force can be calulate by accumulate each small triangle
%  A y
%  |
%  |---> x
%  
%  -------------
%  |  |  |  |
%  |  |  |  |
%  |  |  |  |     
%  -------------
%  |  |  |  |
%  |  |  |  |
%  |  |  |  |
%  -------------
%
%      1   2
%      -----
%      |  |
%      |  |
%      |  |
%      -----
%      3   4

for i=1:1:num_x_node-1
    for j=1:1:num_y_node-1
        % calculate force on each sub cube's top surface
        dx= x_range / (num_x_node - 1);
        dy= y_range / (num_y_node - 1);
        area_triangle= dx*dy/2;
        force_on_triangle=area_triangle*pressure;
        num_node_per_triangle = 3;
        f_on_node = force_on_triangle/num_node_per_triangle;
        node_1=grid.nodeIndex(i,j+1,num_z_node);
        node_2=grid.nodeIndex(i+1,j+1,num_z_node);
        node_3=grid.nodeIndex(i,j,num_z_node);
        node_4=grid.nodeIndex(i+1,j,num_z_node);
        % lower triangle
        for node_index = [node_1 node_3 node_4]
            dof_index = 3;
            dof_global_index = (node_index - 1) * dof_per_node + dof_index;
            load_value(dof_global_index)=load_value(dof_global_index)+f_on_node;
        end
        % upper triangle
        for node_index = [node_1 node_2 node_4]
            dof_index = 3;
            dof_global_index = (node_index - 1) * dof_per_node + dof_index;
            load_value(dof_global_index)=load_value(dof_global_index)+f_on_node;
        end
        clear node_1 node_2 node_3 node_4
    end
end

%% apply constrain to global stiffness matrix
P=load_value;
for i=1:1:num_dof
    if is_constrain(i)
        K(i,:)=0;
        K(:,i)=0;
        K(i,i)=1;
        P(i)=0;
    end
end

%% solution
U=KP;

displacement=(reshape(U,dof_per_node,num_node))';
new_node_coordinate_table=node_coordinate_table+displacement*1e8;

figure
hold on
plotTetGrid(node_coordinate_table,element_node_table);
scatter3(new_node_coordinate_table(:,1),new_node_coordinate_table(:,2),new_node_coordinate_table(:,3));
axis equal

x_displacement=zeros(num_node,1);
y_displacement=zeros(num_node,1);
z_displacement=zeros(num_node,1);
for i=1:1:num_node
    x_displacement(i) = U((i - 1)*dof_per_node + 1);
    y_displacement(i) = U((i - 1)*dof_per_node + 2);
    z_displacement(i) = U((i - 1)*dof_per_node + 3);
end
figure
plot(z_displacement);
xlabel('node index');
ylabel('z displacement');
```

Bill Greene · Answer

The reason that this particular mesh does not give the correct, uniform
displacement solution to this problem is that it is "non-conforming."
Specifically, at the intersection of the two cubes in the model, the
two element edges that cross that face don't align with each other but
instead cross each other.
A typical face in a non-conforming mesh are shown in the figure below.
Elements 1 and 2 are on one side of this face and elements 3 and
4 are on the other. Consider an arbitrary location in the interior
of face (element) 1. The displacements at this point are computed
by interpolating the values at nodes 1, 3, and 4. Assume that
face (element) 3 also contains this location. In this face, the
displacements at that location are interpolated from nodes 1, 2, and 3.
Since the displacements at the element nodes are, in general, arbitrary,
the displacements interpolated at the selected location will not be the
same. This is why the mesh is referred to as non-conforming.

I have attached an abaqus model of this part where the mesh is conforming
at this face. Note that this model contains six tetrahedra in each cube
for a total of twelve overall. The nodal displacements are exactly uniform
for this model and the values of $sigma_{zz}$ in all elements are equal to
the applied pressure.
*node, nset=all_nodes
  1,                0,                0,                0
  2,                0,               10,                0
  3,               10,                0,                0
  4,               10,               10,                0
  5,                0,                0,               10
  6,                0,               10,               10
  7,               10,                0,               10
  8,               10,               10,               10
  9,                0,                0,               20
 10,                0,               10,               20
 11,               10,                0,               20
 12,               10,               10,               20
*element, type=c3d4, elset=all_elements
  1,   7,  10,  11,   8
  2,  12,  11,  10,   8
  3,   6,   3,   2,   4
  4,   6,   7,   3,   4
  5,   7,  10,   9,  11
  6,   6,  10,   7,   8
  7,   6,   8,   7,   4
  8,   2,   5,   1,   3
  9,   6,   5,   2,   3
 10,   6,   7,   5,   3
 11,   6,   9,   5,   7
 12,   6,  10,   9,   7
*material, name=the_mat
*elastic
         100.,  .3
*solid section, elset=all_elements, material=the_mat
1.0
*nset,nset=x0
1,2,5,6,9,10,
*nset,nset=y0
1,3,5,7,9,11,
*nset,nset=z0
1,2,3,4,
*nset,nset=top 
9,10,11,12,
*surface, TYPE=CUTTING SURFACE, DEFINITION=COORDINATES, name=top_surf
0.,0.,20.,   0.,0.,1.
all_elements
*Boundary
x0, 1, 1
y0, 2,2
z0, 3, 3
*Step, name=Step-1, nlgeom=NO
*Static
1., 1., 1e-05, 1.
*Dsload
**Surf-1, P, -1.
top_surf,P,-1.
*node print
u
*el print
S
*End Step

The figure below shows a contour plot of the z-displacement which, as can be seen,
is uniform in the x and y directions.

Xu Hui · Answer

Well, there is no bug in my code related to FEM. It seems that differences can appear if the mesh is not good. I checked my code with the software Abaqus. The element matrices and element vectors are all same.
I can obtain the same solution using Abaqus and my code. The displacement is not constant.

Work condition

result

One can see the displacement in $z$ is not uniform in the result obtained by this mesh.
The following input file for Abaqus can be used to reproduce the solution.
*Heading
** Job name: Job-test Model name: Job_my
** Generated by: Abaqus/CAE 2016
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=PART-1
*Node
      1,           0.,           0.,           0.
      2,           0.,          10.,           0.
      3,          10.,           0.,           0.
      4,          10.,          10.,           0.
      5,           0.,           0.,          10.
      6,           0.,          10.,          10.
      7,          10.,           0.,          10.
      8,          10.,          10.,          10.
      9,           0.,           0.,          20.
     10,           0.,          10.,          20.
     11,          10.,           0.,          20.
     12,          10.,          10.,          20.
*Element, type=C3D4
 1,  1,  3,  4,  7
 2,  1,  4,  2,  6
 3,  7,  5,  6,  1
 4,  7,  6,  8,  4
 5,  1,  7,  4,  6
 6,  5,  7,  8, 11
 7,  5,  8,  6, 10
 8, 11,  9, 10,  5
 9, 11, 10, 12,  8
10,  5, 11,  8, 10
*Elset, elset=Set-1, generate
  1,  10,   1
** Section: Section-1
*Solid Section, elset=Set-1, material=MATERIAL-1
,
*End Part
**  
**
** ASSEMBLY
**
*Assembly, name=Assembly
**  
*Instance, name=PART-1-1, part=PART-1
*End Instance
**  
*Nset, nset=Set-1, instance=PART-1-1, generate
 1,  4,  1
*Nset, nset=Set-2, instance=PART-1-1, generate
  1,  11,   2
*Nset, nset=Set-3, instance=PART-1-1
  1,  2,  5,  6,  9, 10
*Elset, elset=_Surf-1_S1, internal, instance=PART-1-1
 8, 9
*Surface, type=ELEMENT, name=Surf-1
_Surf-1_S1, S1
*End Assembly
** 
** MATERIALS
** 
*Material, name=MATERIAL-1
*Elastic
100.,0.
** 
** BOUNDARY CONDITIONS
** 
** Name: BC-1 Type: Displacement/Rotation
*Boundary
Set-1, 3, 3
** Name: BC-2 Type: Displacement/Rotation
*Boundary
Set-2, 2, 2
** Name: BC-3 Type: Displacement/Rotation
*Boundary
Set-3, 1, 1
** ----------------------------------------------------------------
** 
** STEP: Step-1
** 
*Step, name=Step-1, nlgeom=NO
*Static
1., 1., 1e-05, 1.
** 
** LOADS
** 
** Name: Load-1   Type: Pressure
*Dsload
Surf-1, P, -1.
** 
** OUTPUT REQUESTS
** 
*Restart, write, frequency=0
** 
** FIELD OUTPUT: F-Output-1
** 
*Output, field, variable=PRESELECT
** 
** HISTORY OUTPUT: H-Output-1
** 
*Output, history, variable=PRESELECT
*End Step

What I learned from this course is that, if a strange bug can happen, we can using commercial software as for comparison purposes.

Uniaxial stretching solution not uniform in FEM code

Mesh

Boundary conditions

Problems

Main code

2 Answers

Add your own answers!

Ask a Question