有限元基础编程 ｜高阶单元计算环形区域惯性矩

公众号、B站：[易木木响叮当]
关注可了解更多的有限元数值仿真技巧。问题或建议，请公众号留言;
如果你觉得木木同学对你有帮助，欢迎赞赏。

本次给大家分享的是：如何使用有限元方法近似估计环形区域惯性矩？

主要围绕以下内容进行展开：

问题描述

采用高斯积分方法计算如下图所示的环形区域相对 轴的惯性矩，计算公式：

整个区域数值积分计算

从上面的公式可以直观的看到，该积分是对二维区域进行积分，本小节就针对于高斯积分如何在二维区域进行积分详细展开讲讲，后续读者感兴趣，可基于相同的格式扩展至三维。

针对于环形区域，若使用直角坐标计算，积分上下限比较难取，可在极坐标系下进行较为方便的数值积分，也可以再转换到等参元坐标系下，进行高斯积分方法计算，坐标转换示意图如下图所示：

直角坐标系向极坐标系转换：

极坐标系向等参元坐标系转换：

大家看到这个公式可能会感到奇怪，以往的书中没有标明两者之间的转换关系，两者的坐标关系式怎么得来的呢？

细心的伙伴不妨再看一下，等参元系的范围是 ，可将其代入上述表达式，得到环形区域的范围。

极坐标系下的二重积分

在极坐标中，环形区域的微元面积 可以表示为 ，得到惯性矩公式：

可进一步计算，

得到惯性矩的值为 。

等参元坐标系下的高斯积分

在等参元坐标系下，得到惯性矩公式：

进行高斯积分，可得

高斯点布置： ， 方向有两个高斯积分点， 方向有3个高斯积分点，参考一维情况下三点高斯积分点分布及权重：

带入后，方程可得：

有限元方法离散计算域

如下图所示，左图将环形区域离散为两个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

求解结果如下：

在使用有限元方法逼近的时候，我们可以看到细化后的网格，得出的数值精度会更高一点，若再进一步细化网格，精度会更高。

但若是模型较大时，一味的为了追求精度而增加网格数量，势必会引起计算成本的增加，计算成本与数值精度的取舍，是个值得思考的问题。

本文的主要参考内容及Matlab代码均来自Amar Khennane编写的《Introduction to Finite Element Analysis Using MATLAB and Abaqus》，想要进一步了解有限元编程的小伙伴可以入手，强烈推荐！

完整代码可在公众号后台回复：8节点单元，即可自动获取。