4
$\begingroup$

I write a simple MATLAB code for solving solid FEM problem.

The problem looks like that

(1)      (2)
x-------x
|     / |
|   /   |
| /     |
x-------x
(3)     (4)

Each node has two degrees of freedom(DOF), which stand for $x$ and $y$ displacements.

Boundary condition

I want only constraint node 1's $x,y$ DOF and node 2's $x DOF, that is constraining DOFs 1,2,3 to zero.

Load condition

All nodes have a node force equal to 1.

Result

I found that the global stiffness matrix has a rank 7 but not the rank 8, which leads to the problem can not be solved.

Problem

In the 2D solid FEM case, the 3 DOFs are needed to be constrained for limiting rigid body translation along $x,y$ axis, and the rotaion. But after I constrain 3 DOFs, the problem still can not be solved. Why?

Appendix

The matlab code looks like

clc
clear
close all
%% model
x_range=10;
y_range=10;
num_x_node=2;
num_y_node=2;
num_node=num_x_node*num_y_node;
DOF=2;
%% generate grid
node_coordinate_table=calculateNodeCoordinateTableForRectangleDomain(num_x_node,num_y_node,x_range,y_range);
element_node_table=calculateElementNodeTableForRectangleDomainTriangleGrid(num_x_node,num_y_node);
%% calculate element matrix
K=calcGlobalMatrix(node_coordinate_table,element_node_table,DOF);
%% load
P=ones(num_node*DOF,1);
%% constrain
is_fix_dof=zeros(num_node*DOF,1);
is_fix_dof(1)=1;
is_fix_dof(2)=1;
is_fix_dof(3)=1;
%% modify discrete equation by constrain (cross zero center one method)
for i=1:1:num_node*DOF
    if is_fix_dof(i)==1
        K(i,:)=0;
        K(:,i)=0;
        K(i,i)=1;
        P(i)=0;
    end
end
%% solution
x=K\P
%% plot grid
figure
set(gcf,'position',[50,50,1000,800])
hold on
plotNode(node_coordinate_table);
plotElement(element_node_table,node_coordinate_table);

%% function
function K=calcGlobalMatrix(node_coordinate_table,element_node_table,DOF)
num_node=size(node_coordinate_table,1);
K=zeros(DOF*num_node);
element_matrix_array=calculateElementMatrixArray(node_coordinate_table,element_node_table);
num_element=size(element_node_table,1);
for i=1:1:num_element
    K0=element_matrix_array{i};
    node_global_index_array=element_node_table(i,:);
    for node_i=1:1:3
        for node_j=1:1:3
            for i1=1:1:DOF
                for j1=1:1:DOF
                    node_global_index_i=node_global_index_array(node_i);
                    node_global_index_j=node_global_index_array(node_j);
                    row=(node_i-1)*DOF+i1;
                    col=(node_j-1)*DOF+j1;
                    ROW=(node_global_index_i-1)*DOF+i1;
                    COL=(node_global_index_j-1)*DOF+j1;
                    K(ROW,COL)=K(ROW,COL)+K0(row,col);
                end
            end
        end
    end
end

end

function [element_matrix_array]=calculateElementMatrixArray(node_coordinate_table,element_node_table)
num_element=size(element_node_table,1);
element_matrix_array=cell(num_element,1);
for i=1:1:num_element
    node_a_id=element_node_table(i,1);
    node_b_id=element_node_table(i,2);
    node_c_id=element_node_table(i,3);
    x=zeros(3,1);
    x(1)=node_coordinate_table(node_a_id,1);
    x(2)=node_coordinate_table(node_b_id,1);
    x(3)=node_coordinate_table(node_c_id,1);
    y=zeros(3,1);
    y(1)=node_coordinate_table(node_a_id,2);
    y(2)=node_coordinate_table(node_b_id,2);
    y(3)=node_coordinate_table(node_c_id,2);
    element_matrix_array{i}=calculateElementMatrixByNodeCoordinate(x,y);
end

end

function plotNode(node_coordinate_table)
x=node_coordinate_table(:,1);
y=node_coordinate_table(:,2);
plot(x,y,'o');
num_node=size(node_coordinate_table,1);
for i=1:1:num_node
    node_x=x(i);
    node_y=y(i);
    text(node_x,node_y,num2str(i));
end

end

function plotElement(element_node_table,node_coordinate)
num_element=size(element_node_table,1);
for i=1:1:num_element
    node_a_id=element_node_table(i,1);
    node_b_id=element_node_table(i,2);
    node_c_id=element_node_table(i,3);
    x_a=node_coordinate(node_a_id,1);
    y_a=node_coordinate(node_a_id,2);
    x_b=node_coordinate(node_b_id,1);
    y_b=node_coordinate(node_b_id,2);
    x_c=node_coordinate(node_c_id,1);
    y_c=node_coordinate(node_c_id,2);
    plot([x_a,x_b,x_c,x_a],[y_a,y_b,y_c,y_a],'k-');
    x_o=(x_a+x_b+x_c)/3;
    y_o=(y_a+y_b+y_c)/3;
    text(x_o,y_o,[num2str(i)]);% -1 for plot element index
end
end

function node_coordinate_table=calculateNodeCoordinateTableForRectangleDomain(num_x_node,num_y_node,x_range,y_range)
num_node=num_x_node*num_y_node;
num_x_element=num_x_node-1;
num_y_element=num_y_node-1;
node_coordinate_table=zeros(num_node,2);
n=1;
dx=x_range/num_x_element;
dy=y_range/num_y_element;
for j=1:1:num_y_node
    for i=1:1:num_x_node
        x = (i-1) * dx;
        y = y_range - (j-1) * dy;
        node_coordinate_table(n,1)=x;
        node_coordinate_table(n,2)=y;
        n=n+1;
    end
end
end

function element_node_table=calculateElementNodeTableForRectangleDomainTriangleGrid(num_x_node,num_y_node)
% element node table generate is order in matlab format, that is, id start from 1
num_x_element=num_x_node-1;
num_y_element=num_y_node-1;
num_element=num_x_element*num_y_element*2;
element_node_table = zeros(num_element,3);
n=1;
for j=1:1:num_y_element
    for i=1:1:num_x_element
        %  a--c
        %  | /|
        %  |/ |
        %  b--d
        node_a_id=i+(j-1)*num_x_node;
        node_b_id=node_a_id+num_x_node;
        node_c_id=node_a_id+1;
        node_d_id=node_b_id+1;
        % upper element abc
        element_node_table(n,1)=node_a_id;
        element_node_table(n,2)=node_b_id;
        element_node_table(n,3)=node_c_id;
        n=n+1;
        % lower element cbd
        element_node_table(n,1)=node_c_id;
        element_node_table(n,2)=node_b_id;
        element_node_table(n,3)=node_d_id;
        n=n+1;
    end
end
end

function K0=calculateElementMatrixByNodeCoordinate(x,y)
% reference
% finite element method, wang xucheng chapter 2

% model
E_0=100;
nu_0=0;
Thickness=1;
% calculate element matrix
a=zeros(3,1);
b=zeros(3,1);
c=zeros(3,1);
for i=1:1:3
    if(i==1);j=2;m=3;end
    if(i==2);j=3;m=1;end
    if(i==3);j=1;m=2;end
    a(i)=x(j)*y(m)-x(m)*y(j);
    b(i)=y(j)-y(m);
    c(i)=-x(j)+x(m);
end
A=1/2*(b(i)*c(j)-b(j)*c(i));
K0=zeros(6,6);
for r=1:1:3
    for s=1:1:3
        K0(2*r-1,2*s-1)=b(r)*b(s)+(1-nu_0)/2*c(r)*c(s);
        K0(2*r-1,2*s)=nu_0*b(r)*c(s)+(1-nu_0)/2*c(r)*b(s);
        K0(2*r,2*s-1)=nu_0*c(r)*b(s)+(1-nu_0)/2*b(r)*c(s);
        K0(2*r,2*s)=c(r)*c(s)+(1-nu_0)/2*b(r)*b(s);
    end
end
K0=K0*E_0*Thickness/(4*A*(1-nu_0^2));
end
```
$\endgroup$
5
$\begingroup$

The particular set of constraints you have chosen does not prevent a rigid body rotation about node 1. Thus the stiffness matrix is singular, as you have noted. One way to prevent this rigid body rotation is to set the y-displacement at node 2 to zero. You could also constrain the x-displacement at either node 3 or node 4 to prevent the rotation.

One way to see this is to replace the constraint by a reaction force acting in the same direction as the constraint. Preventing the rotation of the body about node 1 requires that the reaction force due to the constraint create a torque about node 1. A reaction force in the x-direction at node 2 produces no torque.

$\endgroup$
  • $\begingroup$ thanks for your reply. I can not understand why it can not prevent rigid body rotation. In my mind, after I fix 2 dof on node 1, so the system can not have translation. After I fix x dof on node 2, the system can not have rotation. Could you please tell me why? $\endgroup$ – Xu Hui Jan 16 at 13:50
  • 2
    $\begingroup$ If you fix both x and y-dofs on node 1, and the x-dof on node 2, all you fix is the horizontal distance between the two nodes; don't you want to fix the y-dof on node 2 to prevent the rotation? $\endgroup$ – Wolfgang Bangerth Jan 16 at 17:15
  • $\begingroup$ I now can totally understand the issue! :) thanks for your time. $\endgroup$ – Xu Hui Jan 17 at 4:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.