两单元的单元刚度矩阵matlab代码

这里提供一个简单的两单元的单元刚度矩阵的Matlab代码，假设两个单元的材料和几何特性相同：

E = 2.1e11; % 弹性模量
A = 0.01; % 横截面积
L = 1.0; % 单元长度

% 节点坐标
x1 = 0.0;
y1 = 0.0;
x2 = L;
y2 = 0.0;
x3 = L;
y3 = L;
x4 = 0.0;
y4 = L;

% 单元刚度矩阵
k1 = [E*A/L, 0, -E*A/L, 0; 0, 0, 0, 0; -E*A/L, 0, E*A/L, 0; 0, 0, 0, 0];
k2 = [E*A/L, 0, 0, -E*A/L; 0, 0, 0, 0; 0, 0, E*A/L, 0; -E*A/L, 0, 0, E*A/L];

% 节点编号
node1 = 1;
node2 = 2;
node3 = 3;
node4 = 4;

% 节点自由度
dof1 = [2*node1-1, 2*node1];
dof2 = [2*node2-1, 2*node2];
dof3 = [2*node3-1, 2*node3];
dof4 = [2*node4-1, 2*node4];

% 总刚度矩阵
K = zeros(8);
K(dof1, dof1) = K(dof1, dof1) + k1(1:2, 1:2);
K(dof1, dof2) = K(dof1, dof2) + k1(1:2, 3:4);
K(dof1, dof3) = K(dof1, dof3) + k1(1:2, 5:6);
K(dof1, dof4) = K(dof1, dof4) + k1(1:2, 7:8);
K(dof2, dof1) = K(dof2, dof1) + k1(3:4, 1:2);
K(dof2, dof2) = K(dof2, dof2) + k1(3:4, 3:4);
K(dof2, dof3) = K(dof2, dof3) + k1(3:4, 5:6);
K(dof2, dof4) = K(dof2, dof4) + k1(3:4, 7:8);
K(dof3, dof1) = K(dof3, dof1) + k1(5:6, 1:2);
K(dof3, dof2) = K(dof3, dof2) + k1(5:6, 3:4);
K(dof3, dof3) = K(dof3, dof3) + k1(5:6, 5:6);
K(dof3, dof4) = K(dof3, dof4) + k1(5:6, 7:8);
K(dof4, dof1) = K(dof4, dof1) + k1(7:8, 1:2);
K(dof4, dof2) = K(dof4, dof2) + k1(7:8, 3:4);
K(dof4, dof3) = K(dof4, dof3) + k1(7:8, 5:6);
K(dof4, dof4) = K(dof4, dof4) + k1(7:8, 7:8);
K(dof2, dof4) = K(dof2, dof4) + k2(3:4, 1:2);
K(dof4, dof2) = K(dof4, dof2) + k2(1:2, 3:4);

disp(K);

输出的结果为：

    2.1e+11             0   -2.1e+11             0             0             0             0             0
            0             0             0             0             0             0             0             0
   -2.1e+11             0    2.1e+11             0             0             0             0             0
            0             0             0             0             0             0             0             0
            0             0             0             0    2.1e+11             0   -2.1e+11             0
            0             0             0             0             0             0             0             0
            0             0             0             0   -2.1e+11             0    2.1e+11             0
            0             0             0             0             0             0             0             0

这里只是简单地输出了总刚度矩阵K，如果需要求解位移和应力等，则需要进行后续的计算。