我正在尝试用 FEA 数值求解悬臂梁的位移。梁被建模为由一组 8 节点六面体单元组成的 3D 实体，这些单元处于未变形状态，构造为立方体。每个元素在求解位移之前具有相同的体积，其中每边为 1 厘米。见下图：
我在正 z 方向（全局坐标）上的梁顶部的 4 个最外侧节点处施加均匀分布的载荷。我得到的结果是这个可怕的变形结构！
显然我的解决方案有问题，梁应该只是弯曲而没有这些可怕的伪影。我现在将尝试解释到目前为止我所考虑的事情，对我来说是裸露的：
节点编号
对于每个元素，我从左下角开始逆时针对 8 个节点（全局节点编号）进行编号，对“底面”进行完整旋转，然后对“顶面”进行相同的旋转，（见下图） . 通过这个编号模式，我设法为每个单元刚度矩阵获得了一个正雅可比行列式。它确实会变成负数，正如预期的那样，如果编号是顺时针进行的。然而，雅可比非常小，实际上它对于每个单元刚度矩阵几乎是恒定的，例如 1.250000000000003e-07，仅与最后一个小数点略有不同。
元素矩阵的特征值
在测试刚度矩阵的特征值时，许多单元刚度矩阵共有 6 个为零的特征值，对应于此处类似帖子的答案中所述的刚体位移。在那篇文章中，回答者指出：
如果你得到超过六个，你的刚度矩阵没有必要的等级，因此是有缺陷的。
我的元素矩阵是 24x24，因为每个节点都有 3 个自由度。如果我只用 3 个自由度对所有内容进行建模，是否应该总共有 6 个特征值为零？我认为额外的 3 个值也与旋转自由度有关，即使我只用 3 个平移自由度模拟了我的问题。
我正在创建在 3 个方向中的每个方向上总共有2 个高斯点的元素矩阵：因为形状函数是线性的，所以高斯点的数量应该足够了。
for i = 1:numgp
for j = 1:numgp
for k = 1:numgp
xi = gp(i);
eta = gp(j);
zeta = gp(k);
%In here I compute Jacobian, strain matrix
mat = mat + w(i)*w(j)*w(k)*B'*c_mat*B*J;
end
end
end
并计算 8 个形状函数中每一个在自然坐标中的导数：例如第一个形状函数的导数：
dN1dXI = -(1/8)*(eta-1)*(zeta-1);
dN1dETA = -(1/8)*(xi-1)*(zeta-1);
dN1dZETA = -(1/8)*(xi-1)*(eta-1);
并构造梯度矩阵如下：
%Gradient matrix of shape functions N1-N8
GN = [
dN1dXI, dN2dXI, dN3dXI, dN4dXI, dN5dXI, dN6dXI, dN7dXI, dN8dXI;
dN1dETA, dN2dETA, dN3dETA, dN4dETA, dN5dETA, dN6dETA, dN7dETA, dN8dETA;
dN1dZETA, dN2dZETA, dN3dZETA, dN4dZETA, dN5dZETA, dN6dZETA, dN7dZETA, dN8dZETA;
];
%Jacobian matrix
%C is a 8x3 matrix with the global (cartesian) coordinates
%the columns x y and z where the i:th row is the i:th node in the element.
Jmat = GN*C;
%Jacobian
J = det(Jmat)
Bmat = Jmat\GN;
B1x = Bmat(1, 1); B2x = Bmat(1, 2); B3x = Bmat(1, 3); B4x = Bmat(1, 4);
B5x = Bmat(1, 5); B6x = Bmat(1, 6); B7x = Bmat(1, 7); B8x = Bmat(1, 8);
B1y = Bmat(2, 1); B2y = Bmat(2,2); B3y = Bmat(2,3); B4y = Bmat(2, 4);
B5y = Bmat(2, 5); B6y = Bmat(2, 6); B7y = Bmat(2,7); B8y = Bmat(2,8);
B1z = Bmat(3, 1); B2z = Bmat(3, 2); B3z = Bmat(3, 3); B4z = Bmat(3, 4);
B5z = Bmat(3, 5); B6z = Bmat(3, 6); B7z = Bmat(3, 7); B8z = Bmat(3,8);
B1 = [B1x 0 0; 0 B1y 0; 0 0 B1z; 0 B1z B1y; B1z 0 B1x; B1y B1x 0];
B2 = [B2x 0 0; 0 B2y 0; 0 0 B2z; 0 B2z B2y; B2z 0 B2x; B2y B2x 0];
%..and so on for all the 8 nods
B = [B1 B2 B3 B4 B5 B6 B7 B8];
然而，这不会导致数值对称的单元刚度矩阵。它们几乎是对称的，但值略有不同，导致它们不完全对称。
边界条件
在求解位移之前，我简单地删除了我组装的（全局刚度矩阵）中的前 48 行和列，它们对应于悬臂梁的狄利克雷边界条件，也就是说，总共有 16 个节点（每个节点有 3平移自由度）应该具有零位移。
我很感激我能得到的任何帮助，我可以包含更多的代码片段，但正如你可能理解的那样，有很多代码，所以我会像现在一样保持它的精简。