单元积分点应力如何外插至节点上 | 数值实现篇

本次分享单元积分点应力外插至节点处的数值实现过程。

应力外插

核心理念：坐标系的转换。

假设    是母单元的自然坐标系，    是由高斯积分点控制的坐标系（术语可能不专业），假设高斯积分方案为    。坐标系转换关系：

单元内任一点的应力    ，由4个高斯积分点应力    进行插值时，可表示为

其中，    是基于高斯积分点的形函数，第一个积分点的坐标在母单元坐标系下为（-1，-1），根据上述的坐标系转换的方式，在高斯积分点的坐标系下，第一个单元节点在高斯积分点坐标系下坐标为    ，将此坐标值代入第一个形函数，得    ，相同的道理，可推导至四个节点在4个形函数下的    外插矩阵：

数值实现

借助以上理论，我们可以基于matlab平台编制以下代码段：

% 将积分点应力外插至单元节点上，这里只列举了Q4的情况
for i = 1:3
  StressElem(e,:,i) = [1+0.5*sqrt(3) -0.5 1-0.5*sqrt(3) -0.5;
  -0.5 1+0.5*sqrt(3) -0.5 1-0.5*sqrt(3);
  1-0.5*sqrt(3) -0.5 1+0.5*sqrt(3) -0.5;
  -0.5 1-0.5*sqrt(3) -0.5 1+0.5*sqrt(3)]*...
  [stress(e,1,i);stress(e,2,i);stress(e,3,i);stress(e,4,i)];
end

对标Abaqus

模型材料参数为普通的线弹性材料，单元类型选择CPS4，网格划分及边界条件设置如下：

在结果对标过程中，可以先对比自研程序与Abaqus的节点位移场：

 

 

在位移场一致的前提下，我们再来对标应力结果。以常见的mises应力为例：

 

 

结果是一致的，说明了程序的正确性。

如果我们还想看一下细节方面的，以1号单元的节点应力s11为例：

自研程序与Abaqus的结果也是一致的，在提取Abaqus单元节点应力时，应该将应力平滑选项取消勾选，即：

单元积分点应力外插matlab函数

function [StressElem,StressNode] = QuadNodeStress(node, element, prop, U, averageType,elemType,guassType)
% 通过节点位移计算节点应力，正应力：Sxx、Syy、Sxy、VonMises
% 增加节点应力均匀化标识：averageType，==1时，采用绕节点直接平均，==2时采用绕节点面积加权平均
    E = prop(1);
    NU = prop(2);
    ID = prop(4);

    [numberNodes, ~] = size(node);
    [numberElements, ~] = size(element);
    StressElem = zeros(numberElements, 3); % 只计算出正应力Sxx、Syy、Sxy即可
    StressNode = zeros(numberNodes, 4);
    WeightSum = zeros(numberNodes, 1);  % 用于加权平均的权重总和
    % 根据平面应力/应变状态ID选择应力-应变矩阵
    if ID == 1
        D = (E/(1-NU^2)) * [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

    % quadrature according to quadType
    [gaussWeights,gaussLocations_cols]=gauss(guassType);
    stress = zeros(numberElements,size(gaussLocations_cols,1),3);
    StressElem = zeros(numberElements,4,3);
    elementDof = zeros(1,2*4);
    % 遍历所有单元计算单元应力
    for e = 1:numberElements
        indice = element(e,:);
        elementDof(1:2:end)=2*indice-1;
        elementDof(2:2:end)=2*indice;
        elementNode = element(e, :);
        elemNodeCoordinate = node(elementNode, :);
        elenode = length(elemNodeCoordinate);
        B=zeros(3,2*elenode);
        for q = 1:size(gaussWeights,1)
            xi_Gauss=gaussLocations_cols(q,1);
            eta_Gauss=gaussLocations_cols(q,2);
            % shape functions and derivatives
            [shapeFunction,naturalDerivatives]=shapeFunctionQuad(xi_Gauss,eta_Gauss,elemType);
            % Jacobian matrix, inverse of Jacobian,
            % derivatives w.r.t. x,y
            [Jacob,XYderivatives] = Jacobian(elemNodeCoordinate,naturalDerivatives);
            A = det(Jacob)*4;
            % B matrix
            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)';
            % element deformation
            strain = B*U(elementDof);
            stress(e,q,1:3) = D*strain;
        end

        % 计算单元应力
        % 将积分点应力外插至单元节点上，这里只列举了Q4的情况
        for i = 1:3
            StressElem(e,:,i) = [1+0.5*sqrt(3) -0.5 1-0.5*sqrt(3) -0.5;
            -0.5 1+0.5*sqrt(3) -0.5 1-0.5*sqrt(3);
            1-0.5*sqrt(3) -0.5 1+0.5*sqrt(3) -0.5;
            -0.5 1-0.5*sqrt(3) -0.5 1+0.5*sqrt(3)]*...
            [stress(e,1,i);stress(e,2,i);stress(e,3,i);stress(e,4,i)];
        end
...

完整版的代码，我将会放置在《有限元基础编程百科全书》有关平面单元的章节，有待更新~