有限元二点高斯积分（高斯积分与减缩积分）

在有限元分析中，为了计算各个单元的刚度矩阵和荷载向量，需要对单元上的积分进行数值近似。两种常见的积分方法是高斯积分和减缩积分。

高斯积分是一种数值积分方法，基于高斯-勒让德公式，用于近似计算函数在给定区间上的积分。在有限元分析中，高斯积分用于计算单元上的刚度矩阵和荷载向量。具体而言，对于每个单元，将其内部积分区域划分为一系列高斯积分点，然后根据高斯积分点的权重和位置，对积分表达式进行数值近似。通过对所有高斯积分点的贡献进行累加，可以得到单元上的刚度矩阵和荷载向量的近似值。

减缩积分（Reduced Integration）是一种在有限元分析中用于减少计算量和存储需求的技术，它采用了较少的高斯积分点来近似计算单元的刚度矩阵。减缩积分的核心思想是，通过在单元内部使用较少的高斯积分点来减少计算量，同时尽可能保持良好的数值稳定性。传统的高斯积分在每个单元上使用一组完整的高斯积分点，这可能会导致计算量过大，尤其是对于高阶元素或复杂几何形状的单元。减缩积分通过选择适当的积分点子集，例如在三维问题中选择积分点子集的体心点，来近似计算刚度矩阵。减缩积分的优点是减少了高斯积分点的数量，从而降低了计算量和存储需求。然而，它也可能导致计算结果的不精确性。由于减少了积分点的数量，减缩积分可能无法捕捉到一些细节信息，特别是在处理某些特殊情况（如奇异问题）时。总结来说，减缩积分是一种在有限元分析中用于减少计算量和存储需求的技术，它通过使用较少的高斯积分点来近似计算刚度矩阵。它可以在一定程度上降低计算开销，但也需要注意其对计算结果的影响。

以平面应力问题为例，使用mathematical计算了四边形单元的刚度矩阵，对比了精确解，高斯积分解，减缩积分解的差异性。

Clear["Global`*"] (*Parameters*) t=1; Em=1; v=0.2; (*mathematical计算默认都是精确解，但是我们提供的数据必须是分数；其次，对于大规模矩阵运算，精确解的求解计算非常慢*) Coord={{0,0},{3,2},{4,4},{2,5}}; (*Constitutive Law-Plain Stress*) Cc=(Em/(1 v))*{{(1/(1-v)),(v/(1-v)),0},{(v/(1-v)),(1/(1-v)),0},{0,0,1/2}}; (*Shape Functions*) N1=(1-eta)*(1-xi)/4;N2=(1-eta)*(1 xi)/4;N3=(1 eta)*(1 xi)/4;N4=(1 eta)*(1-xi)/4; Ni={N1,N2,N3,N4}; Head[eta]; (*没有赋值直接拿来用的变量属于symbol数据类型*) Head[N4]; (*symbol数据类型的组合是Times数据类型，实际上还是symbol类型，Times代表symbol类型之间运算后赋值给N1*) (*Coordinate Mapping*) CoordMap={x,y}; x=Sum[Ni[[i]]*Coord[[i,1]],{i,1,4}]; y=Sum[Ni[[i]]*Coord[[i,2]],{i,1,4}]; (*Jacobian Matrix*) NatCoord={xi,eta}; J=Table[D[CoordMap[[i]],NatCoord[[j]]],{i,1,2},{j,1,2}]; Jinv=Transpose[Inverse[J]]; H=Jinv//MatrixForm ; (*这里H仅仅是矩阵格式，不是原来的Jinv，这种格式就是用来显示的，一种代表形式的数据类型*) (*B Matrix*) B1={{1,0,0,0},{0,0,0,1},{0,1,1,0}}; B2={{Jinv[[1,1]],Jinv[[1,2]],0,0},{Jinv[[2,1]],Jinv[[2,2]],0,0},{0,0,Jinv[[1,1]],Jinv[[1,2]]},{0,0,Jinv[[2,1]],Jinv[[2,2]]}}; B2//MatrixForm; B3=Table[0,{i,1,4},{j,1,8}]; (*如果是四边形八节点的二次单元，这里改为{j,1,8}*) Do[B3[[j,2i-1]]=B3[[j 2,2i]]=D[Ni[[i]],NatCoord[[j]]],{i,1,4},{j,1,2}]; B3//MatrixForm; (*不要写成B=B1.B2.B3//MatrixForm*) B=B1.B2.B3; (*Stiffness Matrix*) IntFunc=(Transpose[B].Cc.B)*Det[J]; ke=t*NIntegrate[IntFunc,{xi,-1,1},{eta,-1,1}]; ke//MatrixForm; (*如果计算精确解，用Integrate，即，kn=t*Integrate[IntFunc,{xi,-1,1},{eta,-1,1}]，但是数据必须都是分数*) (*下面是用高斯积分,又算了一次刚度矩阵中的积分*) xiset={1/Sqrt[3],-1/Sqrt[3]}; etaset={1/Sqrt[3],-1/Sqrt[3]}; w={1,1}; (*权重系数列表*) kfull=t*Sum[w[[i]]*w[[j]]*(IntFunc/.{xi->xiset[[i]],eta->etaset[[j]]}),{i,1,2},{j,1,2}]; (*下面是用减缩积分,又算了一次刚度矩阵中的积分*) kreduced=t*2*2 IntFunc/.{xi->0,eta->0}; (*对比三种积分方法的差异*) ke//MatrixForm ; kfull//MatrixForm; kreduced//MatrixForm;

刚度矩阵精确结果

刚度矩阵高斯积分结果

刚度矩阵减缩积分结果

如上图所示，高斯积分的结果与解析解几乎一致，减缩积分的误差略大，可见减缩积分的确会降低精度。进一步改变单元形状发现，误差与单元不规则程度相关，例如，单元越接近矩形，精度越高。其中当单元形状为矩形时，高斯积分解等于精确解，这就是有限元软件中强调网格质量的原因（当然加密网格密度也是一种解决方法），或许很多人都因为单元扭曲严重而被终止或者警告分析。当然如果因为网格问题影响收敛时，不妨考虑以下方法：

1. 重新网格划分：如果单元扭曲严重，可能需要重新划分网格，以改善单元的形状和质量。使用更合适的网格生成算法或软件工具，可以生成更好的网格布局，减少扭曲现象的发生。
2. 使用更适合的单元类型：对于特定类型的问题，某些单元类型可能更适合，并能够减少单元扭曲。
3. 局部修复或调整：对于少数扭曲严重的单元，可以尝试通过局部调整来修复它们。这可能包括重新定义节点位置或移动单元的某些节点，以改善其形状。
4. 增强积分（Enhanced Integration）：增强积分是一种用于提高有限元分析精度的技术。通过在高应变或高梯度区域增加高斯积分点，可以提高对这些区域的数值近似。增强积分可以帮助解决扭曲导致的数值不稳定性问题。
5. 错误估计与自适应网格：使用错误估计技术来检测扭曲严重的单元，并基于估计的错误信息进行自适应网格细化或粗化。这可以帮助优化网格，减少扭曲并提高计算结果的准确性。

最近刚复习了一遍有限元分析基础，旨在希望加深理解原理，进一步更好的应用有限元软件，有不当之处，欢迎指正，后续会抽时间更新其他笔记，感兴趣的同学可以看看我的上一篇有限元编程之杆单元，希望大家一起学习交流。