利用MATLAB计算梁单元刚度矩阵并组装成总体刚度矩阵

于结构有限元分析里，梁单元刚度矩阵的计算，以及梁单元刚度矩阵的组装，它是求解框架结构内力与变形的核心步骤，其精度直接决定了位移结果的可靠性，也决定了支反力结果的可靠性。

局部坐标系下的单元刚度矩阵

对于平面梁单元而言，每个节点都含有轴向位移、横向位移以及转角这三个自由度。在局部坐标系当中，单元刚度矩阵是依据材料力学里的梁弯曲理论以及胡克定律推导而来，其表达式取决于弹性模量E、截面惯性矩I以及单元长度L这三个基本参数。

要将矩阵划分成四个 3x3 的子块，这些子块会分别对两端节点力与位移之间的关系予以细致描绘。轴向上的刚度项是 EA/L ，而弯曲刚度项涵盖了诸如 12EI/L³、6EI/L²以及 4EI/L 等具有典型特征的系数。拿长度 L 等于 4 米、弹性模量 E 为 2.1e11 Pa 、惯性矩 I 是 5.2e-5 m⁴的工字钢梁当作例子，在其单元刚度矩阵里弯曲项的数值能够达到 10⁶量级。

于MATLAB里进行实现之际，能够定义函数BeamElementStiffness(E,I,L)，其内部径直依照公式去组装6x6矩阵。在实际开展编程时，务必要留意系数的正负号情况，以此来保证矩阵呈现出对称且半正定的特性。在代码当中，建议运用向量化的值进行赋值，就像k(2,2)=12EI/L^3这样，防止出现逐元素循环的现象。

总体刚度矩阵的组装原理

function k = beam_stiffness(E, I, L)
    k = (E*I/L^3) * [12, 6*L, -12, 6*L;
                    6*L, 4*L^2, -6*L, 2*L^2;
                    -12, -6*L, 12, -6*L;
                    6*L, 2*L^2, -6*L, 4*L^2];
end

将所有单元的贡献叠加到全局坐标系下，这是结构分析所需要做的，假设有一个平面刚架，它存在N个节点，每个节点都保持着3个自由度，那么总体刚度矩阵的维度就是3N×3N，这一矩阵具备稀疏、带状以及对称这三大特征，非零元素仅仅出现在自由度相互关联的位置上。

首先，进行组装流程，要初始化全是零的矩阵K_global，此矩阵等于zeros(3N, 3N)。而后，逐个遍历每一个单元，并从中提取该单元两端节点相关的全局编号。对于单元来说，其连接着节点i与节点j，它的局部刚度矩阵需要利用坐标转换矩阵T，把它转向全局方向以使转换生效才可以继续下去，当然，这里的转换是和单元轴线与全局X轴的夹角余弦该值有关联的哦。

拿某两层框架当作示例而言，底层柱节点的编号是从1至4，顶层节点的编号为5至8。单元刚度矩阵当中，和节点1以及节点2有关联的项，要累加到K_global的第1 - 3行第1 - 3列里，还要累加到第1 - 3行第4 - 6列等相应位置。这样的过程就被称作“对号入座”叠加法，必须得确保自由度映射的精确性。

边界条件处理的常用方法

在实际的结构里，支座节点常常有着位移约束的情况，比如说当柱脚和基础是刚接的时候，这个节点的水平位移是零，竖向位移也是零，转角同样为零。在对这些边界条件进行处理之际，最为直接的办法是划行划列法，也就是把被约束自由度所对应的行给删除掉，并且把对应的列也删除掉，然而这种方法会让矩阵的维度发生改变。

一种更为时常被运用的方式是置大数法，此方法会把主对角线那里对应的项去乘以一个极大的数值，就像10的15次方这种，与此同时，会把荷载向量所对应的项替换成约束值乘以上述那个大数。比如说，要是约束节点1在x方向的位移是零，那么就把K(1,1)修改为等于K(1,1)乘以1e15，而F(1)就变为0。这种方法不会让矩阵的维数发生改变，从而有利于借助编程去达成。

需采用罚函数法或者引入弹簧单元，针对倾斜支座或者弹性约束。2023年，在某跨海大桥的桥塔分析环节，工程师运用置大数法处理了塔底固定约束，计算得出的塔顶位移跟实测值相比，误差仅仅3.7%。于MATLAB代码里，应预先识别所有固定自由度的编号列表，统一去进行对角放大操作。

荷载施加与位移求解

% 示例：两单元梁结构（节点1-2-3）
E = 210e9; % 弹性模量 (Pa)
I = 5e-6;  % 惯性矩 (m^4)
L1 = 2;    % 单元1长度
L2 = 3;    % 单元2长度
% 计算单元刚度矩阵
k1 = beam_stiffness(E, I, L1);
k2 = beam_stiffness(E, I, L2);
% 总体刚度矩阵组装
N_nodes = 3; % 总节点数
K = zeros(3*N_nodes, 3*N_nodes);
% 定义单元自由度映射
def_map = @(i, j) [3*i-2:3*i, 3*j-2:3*j];
% 组装单元1（节点1-2）
K(def_map(1,2), def_map(1,2)) = K(def_map(1,2), def_map(1,2)) + k1;
% 组装单元2（节点2-3）
K(def_map(2,3), def_map(2,3)) = K(def_map(2,3), def_map(2,3)) + k2;
disp('总体刚度矩阵:');
disp(K);

结构里的集中力、分布力或者弯矩要转变成节点等效荷载向量F，其维度也是3N×1。针对跨中集中荷载，能够依照静力等效准则分配至两端节点；对于均布荷载，两端节点各自承担一半总力并且附加固端弯矩。

经修改边界条件之后的总体刚度矩阵K_mod，以及荷载向量F_mod，二者共同构成了线性方程组K_mod U = F_mod，这里面的U是未知节点位移向量。于MATLAB里头直接运用反斜杠运算符U = K_mod\ F_mod就能够求解，此函数会依据矩阵性质挑选最优算法。

对于那种节点数小于500的中型结构而言，直接去求解的话效率是挺高的呢 ，就拿某座6层钢框架来说。它的节点一共有42个哟 ，自由度数量是126个呢 ，求解位移所花费的时间仅仅只有0.03秒。求解完成之后去提取各个节点的x位移 ，还有y位移以及转角值 ，进一步的话能够计算单元内力。要留意去检查结果是不是出现数量级异常的情况哦 ，要是位移超过了L/100的话 ，那么很有可能结构刚度是不足的。

% 固定节点1的所有自由度
fixed_dofs = [1,2,3]; % 节点1的3个自由度
free_dofs = setdiff(1:3*N_nodes, fixed_dofs);
% 矩阵缩聚
K_reduced = K(free_dofs, free_dofs);

支反力的后处理计算

在将所有节点位移都求出来了之后，就能够针对支座反力进行反向计算。对于那些被约束起来的自由度而言，其支反力是这样的情况，等于该自由度所对应的刚度矩阵行向量去跟整体位移向量相乘，之后得到的结果再去减掉该行之荷载值。其数学表达式呈现成为这样的形式，设为 R 等于 K_global(约束行, :) 与 U 相乘之后的结果再减去 F(约束行)。

当于MATLAB里进行实现这个操作时，首先要去提取那约束自由度的逻辑索引，然后再开展计算反力向量R的工作。比如说节点1在y方向是处在固定状态的，它对应的全局自由度是第2个，那么就有R(2) = K_global(2,:) U - F(2)这样的式子。所计算得出的反力应当和外荷载一起去满足整体平衡条件，该条件能够作为验证模型正确性的依据。

% 示例荷载（节点3受集中力）
F = zeros(3*N_nodes, 1);
F(3 * 3-2) = 1000; % 节点3的x方向力
% 求解位移
U = K_reduced \ F(free_dofs);
% 扩展位移结果到全节点
U_full = zeros(3*N_nodes, 1);
U_full(free_dofs) = U;
% 计算支反力
F_reactions = K * U_full - F;

在某一实际存在的工程里面，针对一根悬臂梁施加了端点竖向方向的荷载在数量上为1000N，经过相关计算所得到的固定端竖向反力其数值是1000N，并且弯矩的数值为4000N·m，此结果跟理论所给出的解析完全保持一致。要是反力出现不平衡这种情况，一般来说是因为单元刚度矩阵在组装的时候出现错误或者是边界条件存在遗漏所导致的。建议将所有的反力输出出来并且进行求和，以此跟总荷载进行对比从而实现验证。

完整代码框架与稀疏优化

综合上面所说的那些步骤，完整的MATLAB程序应当涵盖参数定义，涵盖节点坐标输入，涵盖单元连接关系定义，涵盖循环组装总体刚度，涵盖处理边界条件，涵盖求解位移，涵盖计算反力。对于超出1000个自由度的结构，应当采用稀疏矩阵存储，运用sparse函数初始化K_global，能够节省90%以上的内存。

% 参数定义
E = 210e9; % 弹性模量 (Pa)
I = 5e-6;  % 惯性矩 (m^4)
nodes = [0,0; 2,0; 5,0]; % 节点坐标 (x,y)
elements = [1,2; 2,3];   % 单元连接关系
% 初始化全局刚度矩阵
N_nodes = size(nodes, 1);
K = zeros(3*N_nodes, 3*N_nodes);
% 遍历单元并组装
for e = 1:size(elements, 1)
    i = elements(e,1); j = elements(e,2);
    L = norm(nodes(j,:) - nodes(i,:));
    k_e = beam_stiffness(E, I, L);
    K(def_map(i,j), def_map(i,j)) = K(def_map(i,j), def_map(i,j)) + k_e;
end
% 边界条件（固定节点1）
fixed_dofs = [1,2,3];
free_dofs = setdiff(1:3*N_nodes, fixed_dofs);
K_reduced = K(free_dofs, free_dofs);
% 施加荷载（节点3受集中力）
F = zeros(3*N_nodes, 1);
F(3 * 3-2) = 1000; % 节点3的x方向力
% 求解位移
U = K_reduced \ F(free_dofs);
U_full = zeros(3*N_nodes, 1);
U_full(free_dofs) = U;
% 输出结果
disp('节点位移:');
disp(U_full);
disp('支反力:');
disp(K * U_full - F);

需要重点留意的事项有：自由度映射表务必精准无误，每个节点的三个自由度要顺序连续地编号；要是梁单元轴线与全局坐标不平行的时候，那就必须依靠旋转矩阵来实施坐标变换；对于铰接节点而言，要释放转动自由度或者采纳特殊单元。提议把单元刚度矩阵计算、坐标转换、组装分别编写成独立的函数，从而方便调试以及复用。

代码示例里头，能够定义节点数组node(i,:)等于[x,y]，单元数组element(i,:)等于[start_node, end_node, E, I, A]，借着循环逐个遍历单元，调用局部刚度矩阵函数，依据节点坐标去计算夹角余弦并组装进到全局矩阵里，最终借助MATLAB的sparse以及backslash达成高效求解。

莫非你认为于编写梁单元有限元程度这个行为当中，最为容易出现差错状况的究竟是自由度映射这个部分，还是坐标交换转化这个流程呢？十分期待你能够在评论区域交流分享你进行调试实验的经历过程，同时去给他个大拇指并对本文加以收藏，从而方便往后能够查阅全部代码的示例状况。

相关标签： # 梁单元 # 刚度矩阵 # MATLAB # 总体刚度矩阵 # 边界条件