续集(一维弹簧单元的直接刚度方法)Python编程和ABAQUS结果对比

大家好！

上一个帖子链接：(一维弹簧单元的直接刚度方法)Python编程和ABAQUS结果对比

上一个帖子我们对比了基于直接刚度法来求解得到一维弹簧单元的刚度矩阵和ABAQUS提取出来的不太一样。

今天我来详细讲一下其中的原理。

例子：还是上一个帖子的例子，它是一个非齐次边界条件的单自由度弹簧系统。

我们知道这个公式:

在整个系统来看，此时

所以我们可以写出

***注意一下，这里的刚度矩阵 [k] 的行列式 |K| =0, 是没有逆矩阵的。

现在我们的目的是想求出U2,U3,U4 这三个位置位移，我们改写一下这个线性方程组

然后移项化简

这时，我们可以删掉U=0的行，以及对应的 [K] 中的列

整理一下

再把求得的位移反带入公式中

这个是解线性方程组的直接解法，利用了矩阵的变换，结果是精确解。在过程中我们发现，原来不可逆的【K】矩阵经过删除行列之后变成了可逆的矩阵。

然而在ABAQUS中，不是这样处理的。

在这一步的时候，我们的解法已经介绍。然而，ABAQUS 运用了补偿法这一巧妙的解法。在边界的节点上补偿一个刚度为kb的弹簧，其中Kb为大刚度系数，具体在公式中体现如下

不用怀疑，理论来讲，方程组中的未知数U2,U3,U4,F1x,F5x的结果没变。这个时候【K】的行列式|K|≠0，于是【K】有逆矩阵，我们可以直接通过解矩阵方法求解位置向量{U}，

在这里就要注意了，假设我们设Kb = 10^36 N/mm ，我们可以忽略F1x和F5x，所以求得的解都是近似解，解的精确程度取决于Kb取值的大小，Kb越大，结果越精确。

此时再把{U}反带入

求得{F}。

我们编程的结果如下所示

ABAQUS结果和编程结果对比

所以ABAQUS提取的整体刚度矩阵实际是经过补偿后的刚度矩阵，严格来说是不正确的，但是并不影响力和位移的结果。

正确的刚度矩阵如下

如果疑问，欢迎交流和指正。

补偿法的部分代码如下