有限元计算过程中积分点应力如何外插至节点处？【公式推导篇】

本次分享的是：有限元计算过程中，单元积分点应力如何外推至节点？

有关积分点与节点的概念可点击跳转阅读历史推文：有限元基本概念-【节点和积分点】，现科普一下Q4单元、Q8单元、Q9单元的形函数和高斯积分方案。

Q4单元

Q8/9单元

应力外插

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

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

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

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

对于Q8、Q9单元，依然可采用    高斯积分方案（减缩积分）。

相应形函数外插矩阵：

公式推导

为了便于公式的推导，可借助Mathematica符号计算软件，编写如下代码：

(*应用2*2高斯积分方案*)
(*定义形函数（Q4）*)
l1[\[Xi]_] := (1 - \[Xi])/2
l2[\[Xi]_] := (1 + \[Xi])/2
N1[\[Xi]_, \[Eta]_] := l1[\[Xi]]*l1[\[Eta]]
N2[\[Xi]_, \[Eta]_] := l2[\[Xi]]*l1[\[Eta]]
N3[\[Xi]_, \[Eta]_] := l2[\[Xi]]*l2[\[Eta]]
N4[\[Xi]_, \[Eta]_] := l1[\[Xi]]*l2[\[Eta]]

(*将单元坐标系转换至高斯积分点坐标系*)
nodeQ4 = Sqrt[3] {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}};
nodeQ8 = Sqrt[3] {
            {-1, -1}, {1, -1}, {1, 1}, {-1, 1},(*corner nodes*)
                 {0, -1}, {1, 0}, {0, 1}, {-1, 0}     (*mid-
    side nodes*)
                 };
nodeQ9 = Sqrt[3] {
            {-1, -1}, {1, -1}, {1, 1}, {-1, 1},(*corner nodes*)
            {0, -1}, {1, 0}, {0, 1}, {-1, 0},(*mid-side nodes*)
            {0, 0}                               (*center node*)
            };
(*计算高斯积分点参考坐标系下的形函数外插矩阵*)
extrapolationMatrix = 
  Table[{N1[\[Xi], \[Eta]], N2[\[Xi], \[Eta]], N3[\[Xi], \[Eta]], 
      N4[\[Xi], \[Eta]]} /. {\[Xi] -> p[[1]], \[Eta] -> p[[2]]}, {p, 
     nodeQ4}] // Expand;
(*矩阵形式展现*)
MatrixForm[extrapolationMatrix]

 

 

 

以上是有关应力外插的理论知识，如何将之添加到我们的有限元代码中，我会在后续的推文中一步一步数值实现，感谢你的阅读！