基于B_BAR技术的选择性积分法计算C3D8单元的刚度矩阵(Python版)

本次分享的内容是基于B_BAR技术的选择性积分法计算C3D8单元的刚度矩阵。

前言

在以往的自编有限元程序中，按照教科书中的C3D8单元刚度矩阵计算公式，对其进行编程求解，最终的位移场结果和Abaqus的几乎一致，但是唯独刚度矩阵出现很大差异，反复检查，确认程序的正确性，但就是刚度矩阵不一样，当时一度纳闷，后来在Abaqus的官网中找出了根本原因。官方解释如下图：

大致意思就是，使用了the selectively reduced-integration technique，有什么好处呢？

1. prevent mesh locking
2. provides accurate solutions in incompressible or nearly incompressible cases

中间操作了哪个步骤呢：the strain-displacement relation (β-matrix) is modified.

OK，知道了以上原因后，就顺藤摸瓜，查阅文献，弄明白什么是：the selectively reduced-integration technique，最终在Hughes老爷子的书里面找到相关介绍：

网上也有一篇文档介绍BBAR技术：

木木就以上文献，给大家做了个总结，还是那句话，计算原理不需要知道太清楚，会在程序中实现行。

BBAR理论

对于C3D8单元，原有的应变关系矩阵B应该是这样的：

对于某个节点的B矩阵：

现在将    矩阵做一个变换：

其中，

OK，介绍到这里就好了，后续精彩理论内容感兴趣的可以详查Hughes老爷子的书，对于单元刚度矩阵计算的程序，只需要以上公式即可。

在积分点循环时，使用    进行全积分计算，    进行减缩积分计算，也就是说，对    进行8个高斯积分点循环，对于    进行1个高斯积分点计算。

 

python程序

def calElemStiff2(XCoords, nu):
    """
    基于B_BAR技术的选择性积分法计算C3D8单元的刚度矩阵
    """
    GaussPoints = np.array([[-1, -1, 1],
                           [1, -1, 1],
                           [1, 1, 1],
                           [-1, 1, 1],
                           [-1, -1, -1],
                           [1, -1, -1],
                           [1, 1, -1],
                           [-1, 1, -1]])
    GaussPoints = 1.0 / np.sqrt(3) * GuassPoints
    Ke = np.zeros([24, 24], dtype=float)
    CC = C0(nu)
    for i in range(8):
        # the deviatoric part of B 进行全积分计算
        xi, eta, zeta = GaussPoints[i]
        B_Dev, detJ = BMechDev(XCoords, xi, eta, zeta)
        Ke += np.dot(np.dot(B_Dev.T, CC), B_Dev) * detJ
    # the dilatational part of B 进行减缩积分计算
    xi, eta, zeta = 0.0, 0.0, 0.0
    w1 = 8.0
    B_Dil, detJ = BMechDil(XCoords, xi, eta, zeta)
    Ke += w1 * np.dot(np.dot(B_Dil.T, CC), B_Dil) * detJ
    return Ke

最终与Abaqus刚度矩阵内的元素最大相差max. difference: 5.828670879282072e-16，几乎保持一致！