有限元基础编程 | 梁结构铰节点如何处理？

在梁结构中，常常包含一些铰节点，根据大学时所学的力学常识，这时铰节点处的弯矩为0。如下图所示，在进行有限元分析时，需要将结构离散为两个梁单元去分析，这时需要注意的是，铰节点只能被考虑一次。

 

在上图中，铰节点可以归属到单元1的右端，或者单元2的左端，若单元1的右端考虑铰节点，且单元2的左端考虑铰节点，最终求得的刚度矩阵就会奇异。

现考虑纯弯Euler-Bernoulli梁单元右端铰节点的问题，单元刚度方程可表示为：

对于单元的铰接端，只有位移自由度参与刚度矩阵集成，转动自由度不参与刚度矩阵集成，该自由度属于内部自由度，应在单元层次上进行自由度释放，也称静力凝聚，在平面非协调元章节中着重介绍了该数值技术。

将纯弯梁单元刚度矩阵可分块写为：

进一步写为：

将上式展开：

基于上式的第二个方程解出    ：

代入第一个方程：

进一步表示为：

其中    就是凝聚矩阵，可变为：

在铰接端弯矩为0，但是转角不等于0，可将刚度矩阵转角自由度的行和列的元素置为0即可，单元1的凝聚矩阵可表示为：

单元2的凝聚矩阵可表示为：

在对包含铰节点的梁单元实际编程时，可参考《Introduction to finite elementanalysis using MATLAB® andabaqus》的处理方式。

首先，创建矩阵Hinge(nel,2)并将所有元素初始化为1，假设单元     的铰节点在其左端，就将该单元的第一个节点设置为0：Hinge(i,1)=0；假设单元     的的铰节点在其右端，就将该单元的第二个节点设置为0：Hinge(j,2)=0。

单元刚度矩阵可参考如下代码：

if Hinge(i, 1) == 0
    disp('hinge(1,1)=0')
    kl=[ 3*EI/L^3      0     -3*EI/L^3     3*EI/L^2  ; ...
        0         0          0             0    ; ...
        -3*EI/L^3      0      3*EI/L^3     -3*EI/L^2 ; ...
        3*EI/L^2      0     -3*EI/L^2      3*EI/L   ];    

elseif Hinge(i, 2) == 0
    disp('hinge(1,2)=0')
    kl=[ 3*EI/L^3     3*EI/L^2     -3*EI/L^3    0 ; ...
        3*EI/L^2     3*EI/L       -3*EI/L^2    0 ; ...
        -3*EI/L^3    -3*EI/L^2      3*EI/L^3    0 ; ...
        0          0              0        0 ] ;   
else
    kl=[ 12*EI/L^3     6*EI/L^2     -12*EI/L^3     6*EI/L^2 ; ...
        6*EI/L^2     4*EI/L        -6*EI/L^2     2*EI/L   ; ...
        -12*EI/L^3    -6*EI/L^2      12*EI/L^3     -6*EI/L^2 ; ...
        6*EI/L^2      2*EI/L        -6*EI/L^2      4*EI/L   ];
end