《有限元基础编程百科全书》更新：考虑剪切变形的Euler-Bernoulli梁单元

Euler-Bernoulli梁的假设中，忽略了剪切变形的影响，为了考虑该影响，单元类型可选取Timoshenko梁，也可以在Euler-Bernoulli梁假设的基础上进行修正。

《有限单元法—编程与软件应用》和《MATLAB有限元结构动力学分析与工程应用》均有提及考虑剪切变形的Euler-Bernoulli梁单元。

考虑剪切变形后，    方向上的位移变为：

其中，    表示弯曲变形引起的位移，    表示剪切变形引起的位移。单元节点自由度可表示为：

由此得到弯曲部分的刚度矩阵和    和剪切部分的刚度矩阵    ，具体的推导过程可见上面提及的参考文献，写程序时只需要将已有的刚度矩阵写进程序里即可。

在上述描述中，每个节点有三个自由度：    、    、    ，实际计算中，利用单元的平衡方程和几何关系可以消去一个变量，得到两个独立变量：    、    的方程，此时再加入轴向刚度，最终的单元刚度矩阵可表示为：

其中，    ，    为剪应力不均匀系数，矩形截面梁    取0.85，圆形截面梁    取0.89，环形截面梁    取0.53，其他截面类型可参考Abaqus帮助文档相关参数设定。具体的单元程序可见MFEAOOP单元库中的B23M.m。

这里贴上B23M.m的程序，仅供参考！

function k = B23M(prop, elemNodeCoordinate)
    % B23M - 剪切修正的平面 Euler-Bernoulli 梁单元(cubic beam)
    E = prop.E;  % 弹性模量
    v = prop.v;  % 泊松比
    
    Profile = prop.Profile; % 截面类型
     if strcmpi(Profile, 'Rectangular')
        % 矩形截面惯性矩计算
        b = prop.a;  % 宽度
        h = prop.b;  % 高度
        A = b * h;      % 截面面积
        I = (b * h^3) / 12;  % 关于 y 轴的惯性矩 (沿高度方向)
        kappa = 0.85;  % 截面剪切修正系数
    elseif strcmpi(Profile, 'Circular')
        % 圆形截面惯性矩计算
        r = prop.r;  % 半径
        A = pi * r^2;  % 截面面积
        I = (pi * r^4) / 4;  % 关于 y 轴的惯性矩 (与半径平方成比例)
        kappa = 0.89;  % 截面剪切修正系数
    elseif strcmpi(Profile, 'Pipe')
        % 环形截面惯性矩计算
        D = prop.D * 2;  % 外径
        d = prop.d * 2;  % 内径
        A = pi * (D^2 - d^2) / 4;  % 截面面积
        I = pi * (D^4 - d^4) / 64;  % 环形截面惯性矩
        kappa = 0.53;  % 截面剪切修正系数
    else
        error('Unknown section type. Please choose "Rectangular" or "Circular".');
    end

    % 计算杆单元长度
    x = elemNodeCoordinate(:, 1);
    y = elemNodeCoordinate(:, 2);
    L = sqrt((x(2) - x(1))^2 + (y(2) - y(1))^2);

    % 坐标变换矩阵
    l = (x(2) - x(1)) / L;
    m = (y(2) - y(1)) / L;
    R=[l m 0 0 0 0
      -m l 0 0 0 0
       0 0 1 0 0 0
       0 0 0 l m 0
       0 0 0 -m l 0
       0 0 0 0 0 1];

    G = E / (2 * (1 + v));
    b = (12 * E * I * kappa) / (G * A * L^2);
    % 计算局部刚度矩阵
    kLocal = [
        E * A / L,   0,                               0,                                  -E * A / L,    0,                                  0;
        0,           12 * E * I / ((1 + b) * L^3),    6 * E * I / ((1 + b) * L^2),         0,           -12 * E * I / ((1 + b) * L^3),       6 * E * I / ((1 + b) * L^2);
        0,           6 * E * I / ((1 + b) * L^2),     (4 + b) * E * I / ((1 + b) * L),     0,           -6 * E * I / ((1 + b) * L^2),        (2 - b) * E * I / ((1 + b) * L);
       -E * A / L,   0,                               0,                                   E * A / L,    0,                                  0;
        0,          -12 * E * I / ((1 + b) * L^3),   -6 * E * I / ((1 + b) * L^2),         0,            12 * E * I / ((1 + b) * L^3),       -6 * E * I / ((1 + b) * L^2);
        0,           6 * E * I / ((1 + b) * L^2),     (2 - b) * E * I / ((1 + b) * L),     0,            -6 * E * I / ((1 + b) * L^2),        (4 + b) * E * I / ((1 + b) * L)
    ]; 
    % 计算全局刚度矩阵
    k = R' * kLocal * R;
end