使用有限元近似估计区域总降水量 | Hammer积分

今日分享的主要内容：如何使用三角形单元Hammer积分有限元近似估计区域总降水量

以往我们在讲解有限元的时候，总是在计算模型的位移结果，现在给大家来点新鲜的！对于区域内某点的降雨量，我们通过构造试函数，基于有限元的离散思想，进一步求得整体的降雨量。

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

问题描述

如下图所示，在一个三角形区域内布置6个雨量计，雨量计坐标及测得的降雨量如下表所示

估算该地区的总降雨量及降雨面积，采取两种方案计算：

1. 4个线性三角形单元；

2. 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.5    25.   ; ... %  Node 2   
        110.    35.   ; ... %  Node 3   
        87.5    70.   ; ... %  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 = [1   2   6; ... % Element 1
          2   3   4; ... % Element 2
          2   4   6; ... % Element 3
          6   4   5];    % Element 4
%
AT = 0. ; % Initialise total area to zero
QT = 0. ; % Initialise total rainfall to zero
%
for i=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.5    25.   ; ... %  Node 2   
        110.    35.   ; ... %  Node 3   
        87.5    70.   ; ... %  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 = [1   2   3    4    5    6]; ... % Element 1
%
AT = 0. ; % Initialise total area to zero
QT = 0. ; % Initialise total rainfall to zero
%
for i=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));
        for j=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    ，结果一致。