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

5月前浏览3019

本次给大家分享的是：有限元数值编程中面对几乎不可压问题时，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方法不需要使用减缩积分也能达到比较好的精度。

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

表3.5是针对于本次问题的精度对比，不代表所有问题中选择性减缩积分都比带沙漏控制的减缩积分方案精度逊色，在应用实际问题时，还需具体问题具体分析，本着学习的初衷可将这些方案都学习一下，实际用的过程中，还需灵活运用。