使用有限元近似估计区域总降水量 | Hammer积分
今日分享的主要内容：如何使用三角形单元Hammer积分有限元近似估计区域总降水量
以往我们在讲解有限元的时候，总是在计算模型的位移结果，现在给大家来点新鲜的！对于区域内某点的降雨量，我们通过构造试函数，基于有限元的离散思想，进一步求得整体的降雨量。
本文的主要参考内容及Matlab代码均来自Amar Khennane编写的《Introduction to Finite Element Analysis Using MATLAB and Abaqus》，想要进一步了解有限元编程的小伙伴可以入手，强烈推荐！
问题描述
如下图所示，在一个三角形区域内布置6个雨量计，雨量计坐标及测得的降雨量如下表所示
| 雨量计 | x(km) | y(km) | 降水量（mm） |
|---|---|---|---|
| 1 | 15 | 15 | 20 |
| 2 | 62.5 | 25 | 15 |
| 3 | 110 | 35 | 10 |
| 4 | 87.5 | 70 | 20 |
| 5 | 65 | 105 | 30 |
| 6 | 40 | 60 | 25 |
估算该地区的总降雨量及降雨面积，采取两种方案计算：
4个线性三角形单元；
1个二次三角形单元。
问题求解
我们首先构造一个试探函数 ，用于描述 区域内某点的降雨量，总降雨量可表示为：
三角形节点分布如下，左为线性三节点，右图为二次六节点：
4个线性三角形单元
对于某个三角形单元内部某点的降雨量
该单元降雨量可表示为：
在计算三角形单元上的数值积分，比较方便的方法是引进面积坐标。
线性三角形单元面积坐标
如下图所示，三角形单元内，任意一点面积坐标
任何方向移动点
以上就是给大家引入的面积坐标地概念，要耐心往下看哦，方便数值积分的地方马上就来啦~
三角形单元面积坐标的数值积分（面积积分）
对于线性三角形单元的面积坐标阶次均为1。其中三角形面积公式，可由节点坐标计算：
故单元降雨量又可写为：
Matlab代码
nel = 4% Total number of elements
nnd = 6% Total number of nodes
nne = 3% Number of nodes per element
%
% Coordinates of the rain guages (nodes)in km
%
geom = [15.15. ; ... % Node 1
62.525. ; ... % Node 2
110.35. ; ... % Node 3
87.570. ; ... % Node 4
65.105. ; ... % Node 5
40.60.] ; % Node 6
%
% Precipiations recorded by the rain guages
%
q = [20.; 15.; 10.; 20.; 30.; 25.] ;
%
% Connectivity
%
connec = [126; ... % Element 1
234; ... % Element 2
246; ... % Element 3
645]; % Element 4
%
AT = 0. ; % Initialise total area to zero
QT = 0. ; % Initialise total rainfall to zero
%
fori=1:nel
%
% for each element retrieve the x and y coordinates of its nodes
%
xi = geom(connec(i,1),1); yi = geom(connec(i,1),2);
xj = geom(connec(i,2),1); yj = geom(connec(i,2),2);
xk = geom(connec(i,3),1); yk = geom(connec(i,3),2);
%
% Retrieve the precipitations recorded at its nodes
%
qi = q(connec(i,1)); qj =q(connec(i,2)); qk =q(connec(i,3));
%
% calculate its area
%
A = (0.5)*det([1 xi yi;...
1 xj yj;...
1 xk yk]);
AT = AT + A;
%
% Estimate quatity of rain over its area
%
Q = (qi+qj+qk)*A/3;
QT = QT + Q;
end
%
AT
QT
1个二次三角形单元
我们可以构造六节点单元内部某点的试函数：
形函数如下式：
该区域的面积为：
Hammer数值积分
我从书里面挑了两个具有代表性的Hammer积分点位置，如下图所示，与高斯点分布不同哦~，在本题中，用到了7个积分点，具体的位置，我还没找到，直接按照书上写的数值来吧，会用就行......
积分点坐标即相应的权重系数可在书中查找，此处不再赘述。
总的降雨量：
其中雅可比矩阵：
Matlab代码
nel = 1% Total number of elements
nnd = 6% Total number of nodes
nne = 6% Number of nodes per element
npt=7;
samp=hammer(npt);
%
% Coordinates of the rain guages (nodes)in km
%
geom = [15.15. ; ... % Node 1
62.525. ; ... % Node 2
110.35. ; ... % Node 3
87.570. ; ... % Node 4
65.105. ; ... % Node 5
40.60.] ; % Node 6
%
% Precipiations recorded by the rain guages
%
q = [20.; 15.; 10.; 20.; 30.; 25.] ;
%
% Connectivity
%
connec = [123456]; ... % Element 1
%
AT = 0. ; % Initialise total area to zero
QT = 0. ; % Initialise total rainfall to zero
%
fori=1:nel
%
% for each element retrieve the vector qe containing the precipitations
% at its nodes as well as the matrix coord containing
% the x and y coordinates of its nodes
%
for k=1: nne
qe(k) = q(connec(i,k));
forj=1:2
coord(k,j)=geom(connec(i,k),j);
end
end
for ig = 1:npt
WI = samp(ig,3);
[der,fun] = fmT6_quad(samp, ig);
JAC = der*coord;
DET = det(JAC);
%
% calculate its area
%
AT = AT+ WI*DET;
%
% Estimate quatity of rain over its area
%
QT = QT + WI*dot(fun,qe)*DET;
end
end
%
AT
QT
其中，hammer.m函数用于计算Hammer积分点及权重值，函数fmT6_quad.m用于计算形函数偏导。
两种方案计算得到的区域面积为3775 ，总降雨量为75500 ，结果一致。
