有限元基础编程 | “选择性减缩积分”理论推导&数值实现

最新推荐文章于 2024-09-12 19:02:20 发布

易木木木响叮当

最新推荐文章于 2024-09-12 19:02:20 发布

阅读量874

点赞数 26

分类专栏：有限元编程文章标签：算法

本文链接：https://blog.csdn.net/qq_45715554/article/details/140918863

版权

有限元编程专栏收录该内容

11 篇文章 1 订阅

订阅专栏

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

本节主要分享选择性减缩积分技术，B-BAR的具体数值实现留作下一篇推文介绍，本次主要内容有：

何为体积自锁？
选择性减缩积分理论公式推导
数值案例
精度对比

何为体积自锁？

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

阐释体积自锁现象，参考《有限元法基础》P193

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

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

选择性减缩积分

公式推导

从这个名字来看，顾名思义就是要选择性的进行减缩积分计算，不要被名字所吓到，静下心一点一点学习这个概念，学习前辈们的智慧。

单元刚度矩阵 $K$ 可以分解为偏刚度矩阵 $K_{dev}$ （deviatoric stiffness matrix）和体积刚度矩阵 $K_{dil}$ （volumetric/dilatational
stiffness matrix），类似于弹性力学中将应力分解为偏应力和静水压力。

$K = K_{dev}+K_{dil}$

对于几乎不可压缩问题，体积模量趋于无穷大，体积刚度矩阵 $K_{dil}$ 远远大于偏刚度矩阵 $K_{dev}$ ，若要避免自锁（体积自锁或剪切自锁）现象，必须使得体积刚度矩阵阵 $K_{dil}$ 发生奇异时（有非零解），相关理论解释可翻阅《有限元法基础》P275。

在数值积分时，对于偏刚度矩阵 $K_{dev}$ 可采用完全积分，对于体积刚度矩阵 $K_{dil}$ 可采用减缩积分，保证其矩阵奇异。

对于体积自锁现象，用户也可使用全局减缩积分，但可能会出现虚假的沙漏模态，此时需要添加沙漏刚度避免沙漏现象，选择性减缩积分不必需要添加额外的沙漏刚度就可以保证解的收敛性。

接下来就详细展开，这两个刚度矩阵具体该如何求解。应力将其分解为体积应力（ $\sigma _{dil}$ ）和偏应力（ $\sigma _{dev}$ ）：
$\begin{cases} \sigma =\sigma _{dil}+\sigma _{dev}\\ \sigma _{dil}=D_{dil}\varepsilon _{dil}\\ \sigma _{dev}=D_{dev}\varepsilon _{dev}\\ \end{cases}$
其中， $D_{dil}$ 表示体积应变部分的弹性矩阵， $D_{dev}$ 表示偏应变部分的弹性矩阵。

为统一形式， $D_{dev}$ 和 $D_{dil}$ 可写为：

$\begin{cases} D_{dev}=2G\left( \boldsymbol{I}_{\boldsymbol{0}}-\small{\frac{1}{3}\boldsymbol{mm}^T} \right)\\ D_{dil\,\,}=\boldsymbol{m}K_{bulk}\boldsymbol{m}^T\\ \end{cases}$

其中， $G\,\,=\,\,\small{\frac{E}{2\left( 1+\nu \right)}}$ 为剪切模量， $K_{bulk}=\small{\frac{E}{3\left( 1-2\nu \right)}}$ 为体积模量，在三维问题中，

$\boldsymbol{m}=\left[ \begin{matrix} 1& 1& 1& 0& 0& 0\\ \end{matrix} \right] ^T$

$\boldsymbol{I}_{\boldsymbol{0}}=\left[ \begin{matrix} 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0.5& 0& 0\\ 0& 0& 0& 0& 0.5& 0\\ 0& 0& 0& 0& 0& 0.5\\ \end{matrix} \right]$

其显示表示式为：
$D_{de\nu}=2G\begin{bmatrix}\frac{2}{3}&-\frac{1}{3}&-\frac{1}{3}&0&0&0\\-\frac{1}{3}&\frac{2}{3}&-\frac{1}{3}&0&0&0\\-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}&0&0&0\\0&0&0&0.5&0&0\\0&0&0&0&0.5&0\\0&0&0&0&0&0.5\end{bmatrix}$

$D_{dil}=K_{bulk}\begin{bmatrix}1&1&1&0&0&0\\1&1&1&0&0&0\\1&1&1&0&0&0\\0&0&0&0&0&0\\0&0&0&0&0&0\\0&0&0&0&0&0\end{bmatrix}$

在二维问题中，

$\begin{cases} D_{dev}=2G\left( \boldsymbol{I}_{\boldsymbol{0}}-\small{\frac{1}{2}\boldsymbol{mm}^T} \right)\\ D_{dil\,\,}=\boldsymbol{m}K_{bulk}\boldsymbol{m}^T\\ \end{cases}$

$\boldsymbol{m}=\left[ \begin{matrix} 1& 1& 0\\ \end{matrix} \right] ^T$

$\boldsymbol{I}_{\boldsymbol{0}}=\left[ \begin{matrix} 1& 0& 0\\ 0& 1& 0\\ 0& 0& 0.5\\ \end{matrix} \right]$

在以上 $D_{dev}$ 的表达式中，会出现 $\frac{1}{2}$ 和 $\frac{1}{3}$ ，这是因为在做均匀应变处理，该状态下的正应变个数。

以上公式可能比较多，但都已给出显化形式，所以在程序中其实很好编，不需要过多的循环嵌套，把握好哪里需要完全积分哪里需要减缩积分即可：

% 偏刚度完全积分
for n=1:size(gaussWeights,1)
    % 高斯积分点位置
    xi_Gauss=gaussLocations_cols(n,1);
    eta_Gauss=gaussLocations_cols(n,2);
    % 形函数在自然坐标系下的偏导
    [~,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)';
    % 积分点循环计算单元刚度矩阵
    k = k+B.'*D*t*B*gaussWeights(n)*det(Jacob);
end
% 体积刚度减缩积分
xi_Gauss = 0.0;eta_Gauss = 0.0;gaussWeights = 4;
[~,naturalDerivatives]=shapeFunctionQuad(xi_Gauss,eta_Gauss,elemType);
[Jacob,XYderivatives]=Jacobian(elemNodeCoordinate,naturalDerivatives);
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)';
k = k+B.'*m*K*m'*t*B*gaussWeights*det(Jacob);

数值案例

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

$u_{r}=\frac{(1+\nu)a^{2}b^{2}}{E\left(b^{2}-a^{2}\right)}\left[\frac{p_{a}-p_{b}}{r}+(1-2\nu)\frac{p_{a}a^{2}-p_{b}b^{2}}{a^{2}b^{2}}r\right]$

1/4圆环有限元模型

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

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

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

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