三维单元的选择性减缩积分技术【理论+数值实现】

19小时前浏览5

本期内容将为大家分享一个解决有限元领域中“自锁”问题的有效方法—选择性减缩积分。

在之前的推文中有较为详细的介绍二维单元的选择性减缩积分技术，本次将从另外一个视角分享三维单元的数值实现，主要内容如下：

单元刚度矩阵分解
应力分解
弹性矩阵
matlab代码片
数值案例

单元刚度矩阵分解

建议读者在学习本小节之前，先跳转阅读二维单元的选择性减缩积分，本小节对三维部分再做详细介绍。

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

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

在数值积分时，对于偏刚度矩阵可采用完全积分，对于体积刚度矩阵可采用减缩积分，保证其矩阵奇异，选择性减缩积分不必需要添加额外的沙漏刚度可以保证解的收敛性。

应力分解

在二维问题中我们将应力将其分解为体积应力（）和偏应力（）：

其中，表示体积应变部分的弹性矩阵，表示偏应变部分的弹性矩阵。

在三维问题中我们也可以将应力分解为拉梅(Lame)形式 ：

其中，是一个四阶张量，拉梅常数，。

弹性矩阵

这时的可表示为：

可表示为：

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

matlab代码片

% 拉梅常数
lamda = prop.v * prop.E/((1-2*prop.v)*(1+prop.v));
nu = prop.E / (2 * (1 + prop.v));

% 体积应变部分的弹性矩阵 Ddil
Ddil = lamda * [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];
% 偏应变部分的弹性矩阵 Ddev    
Ddev = nu * [2 0 0 0 0 0;
             0 2 0 0 0 0;
             0 0 2 0 0 0;
             0 0 0 1 0 0;
             0 0 0 0 1 0;
             0 0 0 0 0 1];

% Gauss 积分点数量为 2x2x2
ngp = 2; % 使用 2x2x2 积分点
samp = gauss3D(ngp); % 获取高斯积分点的坐标和权重

% 偏刚度矩阵 Kdev 采用完全积分
for ig = 1:ngp
    xi_Gauss = samp(ig, 1); % 第 ig 个积分点的 xi 坐标
    wi = samp(ig, 2);       % 第 ig 个积分点的权重
    for jg = 1:ngp
        eta_Gauss = samp(jg, 1); % 第 jg 个积分点的 eta 坐标
        wj = samp(jg, 2);        % 第 jg 个积分点的权重
        for kg = 1:ngp
            gama_Gauss = samp(kg, 1); % 第 kg 个积分点的 gamma 坐标
            wk = samp(kg, 2);         % 第 kg 个积分点的权重
            
            % 计算 B 矩阵和 Jacobian 矩阵
            [B, Jacob] = Bmatrix3D(elenode, elemNodeCoordinate, DOF, xi_Gauss, eta_Gauss, gama_Gauss);

            % 使用高斯积分公式更新偏刚度矩阵
            % kdev = ∫ B^T * Ddev * B * det(Jacob) dΩ
            k = k + B.' * Ddev * B * det(Jacob) * wi * wj * wk;
        end
    end
end

% 体积刚度矩阵 Kdil 采用减缩积分，仅在单元中心进行计算
xi_Gauss = 0.0; eta_Gauss = 0.0; gama_Gauss = 0.0; W = 8;
[B, Jacob] = Bmatrix3D(elenode, elemNodeCoordinate, DOF, xi_Gauss, eta_Gauss, gama_Gauss);
% kdil = ∫ B^T * Ddil * B * det(Jacob) dΩ
k = k + B.' * Ddil * B * det(Jacob) * W;