繼上次的推文：有限元計算過程中積分點應力如何外插至節點處？【公式推導篇】，本次分享單元積分點應力外插至節點處的數值實現過程。

數值實現

借助以上理論，我們可以基于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的節點位移場：

Abaqus位移場結果

自研程序位移場結果

在位移場一致的前提下，我們再來對標應力結果。以常見的mises應力為例：

Abaqus位移應力場結果

自研程序應力場結果

結果是一致的，說明了程序的正確性。

如果我們還想看一下細節方面的，以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
...

完整版的代碼，我將會放置在《有限元基礎編程百科全書》有關平面單元的章節，有待更新~

覺得本篇推文對你有幫助的話，可以動動的小手一鍵三連（點贊?在看?分享）哦~