有限元基础编程 | B-BAR技术理论公式推导&数值实现

本次给大家分享的是：有限元数值编程中面对几乎不可压问题时，B-BAR技术如何施展？

低阶单元完全积分在分析几乎不可压缩问题（泊松比接近0.5）时，易出现体积自锁，精度是非常差的！为适应此类情况，前辈们拓展了积分方案（减缩积分）以及各种数值技术来适应不可压问题，其中选择性减缩积分和B-BAR技术比较具有代表性。

本节主要分享B-BAR技术，主要内容有：

• 何为体积自锁？
• B-BAR理论公式推导
• 数值案例
• 精度对比

何为体积自锁？

如下图所示的平面应变问题（    ），体积应变应为    ，对于不可压缩或者几乎不可压缩问题（如橡胶材料）来说，其泊松比无限接近于0.5，体积模量接近于无穷大，单元的体积在受力过程中几乎保持不变。

 

上图中的1、2、5号节点的自由度被固定，若要使单元1的体积保持不变，节点6不能沿着x方向移动，同时为了保持单元2的体积不变，节点6同样不能沿着y方向移动，因此就造成了节点6被固定，同样的道理，延伸至其他节点也是被固定，模型在受力过程中，网格各个节点均不能发生移动，产生了"闭锁"现象，有限元世界里常称之为"体积自锁"（volumetriclocking）或"网格自锁"（mesh locking）。

以上是从物理现象角度出发，解释何为体积自锁，其实背后有其数学物理方程来解释，本文档不打算记录该物理过程，感兴趣的读者可翻阅张雄老师的《有限元法基础》。

B-BAR

选择性减缩积分将应力分解为体积应力部分和偏应力部分，适用于各向同性的弹性材料，难以扩展到弹塑性材料、轴对称、各向异性、大变形等问题上。为增加普适性，Hughes将几何矩阵    分解为体应变部分    和偏应变部分    ，然后对体应变部分    做出某种改进（记为    ），称为B-bar方法，即：

公式推导

单元的几何矩阵    可写为：

其中，    表示单元节点个数，三维问题中，子块矩阵    可写为：

其中，    ，    表示与节点    关联的形函数，    表示自由度号，体应变部分    可表示为：

偏应变    ：

为适应不可压缩材料，体应变部分    需要加一个bar：

   的取值可以为单元形心处的值，即：

本节的程序中采用了形心处的值，读者也可采用单元体积内的平均值：

最终的几何矩阵    可转化为    ：

最终的单元刚度矩阵可写为：

写到这里，程序已经可以顺着公式往下编写了，如果需要对矩阵    的各个元素进行具体细化，可以进一步写为：

其中，

对于二维问题，

体应变    可表示为：

此处的    表示该状态下包含两个正应变，后续的推导过程与三维状态一致。

计算    的matlab代码如下：

if ID == 1 % 平面应力状态
   D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2];
elseif ID == 2 % 平面应变状态
   D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2];
end

for n=1:size(gaussWeights,1)
    xi_Gauss=gaussLocations_cols(n,1);
    eta_Gauss=gaussLocations_cols(n,2);
    % 原始B矩阵
    [B,Jacob] = Bmatrix(elenode,elemNodeCoordinate,elemType,xi_Gauss,eta_Gauss);
    % 偏矩阵B
    Bdil = Bdilmatrix(elenode,elemNodeCoordinate,elemType,xi_Gauss,eta_Gauss);
    % 偏矩阵B-bar，取单元形心处的值(0,0)
    BdilBar = Bdilmatrix(elenode,elemNodeCoordinate,elemType,0,0);
    % 体积矩阵
    Bdev = B-Bdil;
    % 最终的B-bar
    Bbar = Bdev+BdilBar;
    k = k+Bbar.'*D*t*Bbar*gaussWeights(n)*det(Jacob);
end

function [B,Jacob] = Bmatrix(elenode,elemNodeCoordinate,elemType,xi_Gauss,eta_Gauss)
    B=zeros(3,2*elenode);
    % 形函数在自然坐标系下的偏导
    [~,naturalDerivatives]=shapeFunctionQuad(xi_Gauss,eta_Gauss,elemType);
    % 计算雅可比矩阵及形函数在物理坐标系下的偏导
    [Jacob,XYderivatives]=Jacobian(elemNodeCoordinate,naturalDerivatives);
    % 计算几何矩阵B
    B(1,1:2:end) = XYderivatives(:,1)';
    B(2,2:2:end) = XYderivatives(:,2)';
    B(3,1:2:end) = XYderivatives(:,2)';
    B(3,2:2:end) = XYderivatives(:,1)';
end

数值案例

数值案例仍然采用沙漏控制章节的案例，现考虑如下图所示的    圆环，下边界约束    方向自由度，左边界约束    方向自由度。内壁（    ）承受大小为1的压强    ，外壁（    ）承受大小为0的压强    ，弹性模量    取13，泊松比    取0.499，单元厚度为1，将之作为平面应变问题看待，径向位移的解析解可表示为：

 

分别使用完全积分、带沙漏控制的全局减缩积分、选择性减缩积分方案、B-bar方法进行有限元分析计算，计算结果汇总于如下图表。

本次案例为了模拟几乎不可压缩材料，将泊松比设置为0.499，从图3.36(a)和表 3.5可以看出，采用完全积分方案计算时，内壁节点的位移精度已经完全失真，可初步得出结论：线性单元完全积分时不能用于模拟几乎不可压材料的力学行为。图3.36(b)-3.36(d)和表3.5可以看出，带沙漏控制的减缩积分方案在应用几乎不可压问题时精度良好，选择性减缩积分和B-bar方法在不带沙漏控制的前提下均可以较好地解决几乎不可压问题，并且B-bar方法不需要使用减缩积分也能达到比较好的精度。

如果对于程序的精度还存在疑问，可进行分片试验，方法在前面章节均已给出，本小节的选择性减缩积分方案可以通过分片试验，读者可自行验证。

 

声明：本篇理论知识来源有张雄老师的《有限元法基础》、Hughes老爷子的《Thefinite element method linear static dynamic finite elementanalysis》，在后台回复Hughes，即可自动获取老爷子的著作。

本次分享的理论解释和matlab代码均上传至《有限元基础编程百科全书》最新版PDF中，后台回复：星球，即可加入获得最新版PDF，对代码的任何疑问欢迎在星球中进行提问。

《有限元基础编程百科全书最新进度》

本章节的缩略图：