我正在为 18 节点 (3x3x2) 3D 元素 FEM 编写代码。然而,即使我(非常)确定所有形状函数都是正确的等等,每当我尝试反转刚度矩阵以求解位移时,我都会收到来自 Matlab 的警告消息,告诉我我的矩阵接近奇异矩阵。
是否有我应该测试以解决此问题的清单?什么样的错误通常会导致这些错误?或者它是否过于宽泛以至于无法做出有意义的回答?
编辑:感谢您的所有意见。很抱歉这么晚才回复,但这是我的代码的链接。本来不想发的,这么烦人的debug,让陌生人给自己做好像有点过分了。
我正在为 18 节点 (3x3x2) 3D 元素 FEM 编写代码。然而,即使我(非常)确定所有形状函数都是正确的等等,每当我尝试反转刚度矩阵以求解位移时,我都会收到来自 Matlab 的警告消息,告诉我我的矩阵接近奇异矩阵。
是否有我应该测试以解决此问题的清单?什么样的错误通常会导致这些错误?或者它是否过于宽泛以至于无法做出有意义的回答?
编辑:感谢您的所有意见。很抱歉这么晚才回复,但这是我的代码的链接。本来不想发的,这么烦人的debug,让陌生人给自己做好像有点过分了。
我的评论特别适用于 3D 结构有限元,但这些原则适用于更一般的元素。
作为单元开发的一部分,您将要执行的测试是计算单元刚度矩阵的特征值和特征向量。对于 3D 结构元素,您应该得到正好六个零(或非常接近零)的特征值。这些的特征向量应该是三个方向上的刚体位移和围绕三个轴的刚体旋转(或这些运动的线性组合)。如果你得到超过六个,你的刚度矩阵没有必要的等级,因此是有缺陷的。听起来这可能是你的问题。(然而,一般来说,如果你得到少于六个零特征值,你的公式也是有缺陷的。)
那么什么会导致这个问题呢?
最可能的原因是形状函数在某些方面无效。第二种可能性是您没有准确地整合刚度矩阵中的各项。在我之前关于类似主题的帖子(3D Solid 8 Node FEM Matlab Code)中,我提到了使用太少集成点的问题。具体来说,问题在于这会导致零频率(零特征值)位移模式,由于特征向量的形状,通常被称为“沙漏”模式。
简单的解决方案是使用足够多的积分点来精确积分刚度矩阵项。还有一些巧妙的数值技巧可以用来消除这些沙漏模式,但仍然允许欠积分。此外,有趣的是,有时具有适当边界条件的缺陷单元网格会相互稳定,从而可以求解完整模型。但是,总的来说,这不是您想要依赖的东西。
元素刚度矩阵必须始终是奇异的,因为在推导它时,我们不施加任何约束或边界条件。因此,反转刚度矩阵以求解位移/位置矢量/自由度应该会产生不确定的结果。