单元积分点应力如何外插至节点上 | 数值实现篇
单元积分点应力如何外插至节点上 | 数值实现篇
在有限元分析中，积分点应力的外插至节点处是一个常见的需求。本文将详细介绍这一过程的数值实现方法，包括坐标系转换、应力插值公式以及具体的MATLAB实现代码。通过与Abaqus软件的结果对比，验证所实现算法的正确性。
应力外插
核心理念：坐标系的转换。
假设$\left( \xi ,\eta \right)$是母单元的自然坐标系，$\left( \hat{\xi},\hat{\eta} \right)$是由高斯积分点控制的坐标系（术语可能不专业），假设高斯积分方案为$2 \times 2$。坐标系转换关系：
$$(\xi,\eta)=(\hat{\xi},\hat{\eta})/\sqrt{3}\quad\mathrm{or}\quad(\hat{\xi},\hat{\eta})=(\xi,\eta)\sqrt{3}$$
单元内任一点的应力$\sigma_P$，由4个高斯积分点应力$\sigma_{Gi}$进行插值时，可表示为：
$$\sigma_{P}=\sum_{i=1}^{4}N_{i}\sigma_{Gi}=\begin{bmatrix}N_{1}(\hat{\xi},\hat{\eta})N_{2}(\hat{\xi},\hat{\eta})N_{3}(\hat{\xi},\hat{\eta})N_{4}(\hat{\xi},\hat{\eta})\end{bmatrix}\begin{bmatrix}\sigma_{G1}\\sigma_{G2}\\sigma_{G3}\\sigma_{G4}\end{bmatrix}$$
其中，$N(\hat{\xi},\hat{\eta})$是基于高斯积分点的形函数，第一个积分点的坐标在母单元坐标系下为（-1，-1），根据上述的坐标系转换的方式，在高斯积分点的坐标系下，第一个单元节点在高斯积分点坐标系下坐标为$\left( -\sqrt{3},-\sqrt{3} \right)$，将此坐标值代入第一个形函数，得$1 + \frac{\sqrt{3}}{2}$，相同的道理，可推导至四个节点在4个形函数下的$4 \times 4$外插矩阵：
$$\begin{bmatrix}\sigma_1\\sigma_2\\sigma_3\\sigma_4\end{bmatrix}=\begin{bmatrix}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}\end{bmatrix}\begin{bmatrix}\sigma_{G1}\\sigma_{G2}\\sigma_{G3}\\sigma_{G4}\end{bmatrix}$$
数值实现
借助以上理论，我们可以基于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号单元的节点应力$s_{11}$为例：
自研程序与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
end
end
完整版的代码，我将会放置在《有限元基础编程百科全书》有关平面单元的章节，有待更新~