# Uniaxial stretching solution not uniform in FEM code

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);

% 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;

% 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;
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;
end
clear node_1 node_2 node_3 node_4
end
end

%% apply constrain to global stiffness matrix
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=K\P;

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');
$$$$

• I corrected the writing in your question, please check that it preserves what you wanted to say. – nicoguaro Jul 16 '20 at 16:43
• Yes, the z-displacement of all the top nodes should be the same. I suggest testing the 8-node, 5-element case. If you can't get the correct solution for that case, there is no point in including more elements in your model. – Bill Greene Jul 16 '20 at 17:03
• Thanks for your modify, it is much clear now @nicoguaro. – Xu Hui Jul 17 '20 at 0:29
• @BillGreene I have already check the 8-node, 5 element, that is nx=2,ny=2,nz=2. At this case, the displacement is uniform. however of nx=2,ny=2,nz=3, the displacement is not uniform. – Xu Hui Jul 17 '20 at 0:35

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.
**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.

• Sir, thanks for your point out the "non-conforming". So in fem, the edge should never cross? And futhermore, I can not understand why the displacement is not match when the edge is "non-conforming". could you please give me a explation? The dof is on the node, and we can solve the problem, so we can obtain the displacement on node. Ath this stage, what is "displacement not match"? I can not understand your words "As a consequence the displacements in the elements on one side of that face don't match the displacements in the adjacent elements on the other side of the face." – Xu Hui Jul 19 '20 at 9:37

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.
**
**
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.

• It seems that you are basing your conclusion on the contours, but if you check the values of the displacement for the top plane you get: 0.209003, 0.215859`. The values are not that different. – nicoguaro Jul 17 '20 at 17:40
• Sir, in theory, the value in theory should be equal but not similar, as the solution obtained by Bill Greene. – Xu Hui Jul 19 '20 at 9:38