有限元基础编程 ｜高阶单元计算环形区域惯性矩
本次给大家分享的是：如何使用有限元方法近似估计环形区域惯性矩？
主要围绕以下内容进行展开：
直角坐标-极坐标-等参元坐标的相互转换； 采用高斯积分法计算惯性矩积分公式； 粗网格（两个8节点单元）离散计算域； 细网格（八个8节点单元）离散计算域。
问题描述
采用高斯积分方法计算如下图所示的环形区域相对
整个区域数值积分计算
从上面的公式可以直观的看到，该积分是对二维区域进行积分，本小节就针对于高斯积分如何在二维区域进行积分详细展开讲讲，后续读者感兴趣，可基于相同的格式扩展至三维。
针对于环形区域，若使用直角坐标计算，积分上下限比较难取，可在极坐标系下进行较为方便的数值积分，也可以再转换到等参元坐标系下，进行高斯积分方法计算，坐标转换示意图如下图所示：
直角坐标系向极坐标系转换：
极坐标系向等参元坐标系转换：
大家看到这个公式可能会感到奇怪，以往的书中没有标明两者之间的转换关系，两者的坐标关系式怎么得来的呢？
细心的伙伴不妨再看一下，等参元系的范围是
极坐标系下的二重积分
在极坐标中，环形区域的微元面积
可进一步计算，
得到惯性矩的值为
等参元坐标系下的高斯积分
在等参元坐标系下，得到惯性矩公式：
进行高斯积分，可得
高斯点布置： ， 方向有两个高斯积分点， 方向有3个高斯积分点，参考一维情况下三点高斯积分点分布及权重：
| 阶数 | 高斯点 | 权重 |
|---|---|---|
| 1 | 0.0 | 2.0 |
| 2 | 0.57735027 | 1.0 |
| 3 | 0.77459667 | 0.55555556 |
| 0.0 | 0.88888889 |
带入后，方程可得：
有限元方法离散计算域
如下图所示，左图将环形区域离散为两个8节点单元，右图将环形区域离散为8个8节点单元。以粗网格模型为例，原环形区域的惯性矩可由离散后的区域的惯性矩相加和所得：
两个单元的惯性矩表示为：
对环形区域内的 坐标，用单元节点形函数插值得到：
微元面积 ，可表示为：
将上式代入后，可得到两个单元的惯性矩：
将上述积分进行离散操作：
其中，8节点平面单元形函数：
高斯积分精度
在进行积分点的选择时，需要保持良好的数值精度，数值分析的教材告诉我们：
假设被积函数是 阶的多项式，应用 点的高斯求积结果是精确的，在应用 个取样点时，实际上就是用 阶的多项式代替被积函数，数值积分结果的精确程度取决于多项式与被积函数曲线的拟合程度。
离散后的公式中，形函数是关于 和 的二次多项式，雅可比行列式由于“偏导操作”，对于 和 是一阶线性的，故该积分涉及的到的多项式阶次最高位5阶，需要在每个方向上布置3个积分点。
Matlab代码
clc
clear
format long e
%
global geom connec nel nne nnd RI RE
%
RI = 40; % Internal radius
RE = 70; % External radius
%
Eight_Q8 % Load input data
% Two_Q8
%
% Number of Gauss points
%
ngp = 3; % The polynomials involved are of degree 5
%
samp = gauss(ngp); % Gauss abscissae and weights
%
%
Ixx = 0.; % Initiliase the second moment of area to zero
%
for k=1:nel
coord = coord_q8(k,nne, geom, connec); % Retrieve the coordinates of
% the nodes of element k
X = coord(:,1); % X coordinates of element k
Y = coord(:,2); % Y coordinates of element k
for i=1:ngp
xi = samp(i,1);
WI = samp(i,2);
for j =1:ngp
eta = samp(j,1);
WJ = samp(j,2);
[der,fun] = fmquad(samp, i,j); % Form the vector of the shape functions
% and the matrix of their derivatives
JAC = der*coord; % Evaluate the Jacobian
DET =det(JAC); % Evaluate determinant of Jacobian matrix
Ixx =Ixx+ (dot(fun,Y))^2*WI*WJ*DET;
end
end
end
Ixx
求解结果如下：
| 精确解 | 粗网格 | 细网格 |
|---|---|---|
| 4211700 | 4205104 | 4211281 |
在使用有限元方法逼近的时候，我们可以看到细化后的网格，得出的数值精度会更高一点，若再进一步细化网格，精度会更高。
但若是模型较大时，一味的为了追求精度而增加网格数量，势必会引起计算成本的增加，计算成本与数值精度的取舍，是个值得思考的问题。
本文的主要参考内容及Matlab代码均来自Amar Khennane编写的《Introduction to Finite Element Analysis Using MATLAB and Abaqus》，想要进一步了解有限元编程的小伙伴可以入手，强烈推荐！
来源：易木木响叮当
