我编写了一个简单的 MATLAB 代码来解决固体 FEM 问题。

问题看起来像这样

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

每个节点有两个自由度（DOF），代表$x$和$y$位移。

边界条件

node 1我只想要约束$x,y$DOF 和node 2's $x DOF，即将 DOF 约束1,2,3为零。

负载条件

所有节点的节点力都等于 1。

结果

我发现全局刚度矩阵的秩是 7 而不是 8，导致问题无法解决。

问题

在 2D 实体 FEM 情况下，需要约束 3 个自由度以限制刚体沿$x,y$轴和旋转。但是我约束了3个自由度后，问题仍然无法解决。为什么？

附录

matlab代码看起来像

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
```