TransWikia.com

Uniaxial stretching solution not uniform in FEM code

Computational Science Asked on August 22, 2021

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

2 Answers

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.

enter image description here

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.

enter image description here

Correct answer by Bill Greene on August 22, 2021

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

enter image description here

  • result

enter image description here

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.

Answered by Xu Hui on August 22, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP