MATLAB 三维梁单元的坐标转换矩阵

在三维梁单元的有限元计算中，我们需要进行坐标转换，将局部坐标系下的应力、应变、刚度矩阵等转换到全局坐标系下。坐标转换矩阵（Coordinate Transformation Matrix）在这个过程中起到了关键作用。以下是 MATLAB 中三维梁单元的坐标转换矩阵示例代码：

%% 定义局部坐标系下的梁单元的节点坐标（三维空间）
node = [
    0, 0, 0;
    0, 0, 1;
    0, 1, 0;
    1, 0, 0;
    1, 1, 0;
    1, 0, 1;
    0, 1, 1;
    1, 1, 1
];

%% 定义局部坐标系下的梁单元的单元拓扑（即节点编号）
elem = [
    1, 2;
    1, 3;
    1, 4;
    2, 7;
    2, 6;
    3, 7;
    3, 5;
    4, 6;
    4, 5;
    5, 8;
    6, 8;
    7, 8
];

%% 定义梁的长度和杨氏模量
L = 1; % 梁的长度
E = 2e11; % 杨氏模量

%% 计算局部坐标系下的梁单元的刚度矩阵
k_local = zeros(6,6);
for i = 1:size(elem,1)
    x1 = node(elem(i,1),1); y1 = node(elem(i,1),2); z1 = node(elem(i,1),3);
    x2 = node(elem(i,2),1); y2 = node(elem(i,2),2); z2 = node(elem(i,2),3);
    L_e = sqrt((x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2); % 计算单元长度
    k_e = E * [1,-1;-1,1] * [A/L_e, -A/L_e; -A/L_e, A/L_e]; % 计算单元刚度矩阵
    k_local([2*i-1,2*i],[2*i-1,2*i]) = k_local([2*i-1,2*i],[2*i-1,2*i]) + k_e; % 将单元刚度矩阵组装到总刚度矩阵中
end

%% 计算局部坐标系到全局坐标系的转换矩阵
theta = pi/4; % 定义旋转角度（弧度）
Rz = [cos(theta), -sin(theta), 0; sin(theta), cos(theta), 0; 0, 0, 1]; % 绕 z 轴旋转的旋转矩阵
Ry = [cos(theta), 0, sin(theta); 0, 1, 0; -sin(theta), 0, cos(theta)]; % 绕 y 轴旋转的旋转矩阵
Rx = [1, 0, 0; 0, cos(theta), -sin(theta); 0, sin(theta), cos(theta)]; % 绕 x 轴旋转的旋转矩阵
R = Rz * Ry * Rx; % 组合三个旋转矩阵得到总的旋转矩阵
T = [R, zeros(3,3); zeros(3,3), R]; % 构造局部坐标系到全局坐标系的转换矩阵

%% 计算全局坐标系下的梁单元的刚度矩阵
k_global = T' * k_local * T;

%% 打印全局坐标系下的梁单元的刚度矩阵
disp('Global stiffness matrix: ');
disp(k_global);

上述代码中，我们首先定义局部坐标系下的梁单元的节点坐标和单元拓扑，然后计算局部坐标系下的梁单元的刚度矩阵。接着，我们定义旋转角度并计算局部坐标系到全局坐标系的转换矩阵 T。最后，将局部坐标系下的刚度矩阵通过转换矩阵 T 转换到全局坐标系下，得到全局坐标系下的梁单元的刚度矩阵 k_global。