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

易公子 2023年1月26日浏览：1434 收藏：2

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

今日分享的主要内容：如何使用三角形单元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个二次三角形单元。

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

问题求解

我们首先构造一个试探函数，用于描述区域内某点的降雨量，总降雨量可表示为：

Q_{T}=\int_{A} Q(x, y) \mathrm{d} A

三角形节点分布如下，左为线性三节点，右图为二次六节点：

使用有限元近似估计区域总降水量 | Hammer积分的图2 — 三角形单元节点分布

4个线性三角形单元

对于某个三角形单元内部某点的降雨量 $\bar{Q}(x, y)$ ，采用节点近似作为试函数：

\bar{Q}(x, y)=N_{i}(x, y) Q_{i}+N_{j}(x, y) Q_{j}+N_{k}(x, y) Q_{k}

该单元降雨量可表示为：

Q_{e}=\int_{A_{1}} \bar{Q}(x, y) \mathrm{d} A=Q_{1} \int_{1} N_{1}(x, y) \mathrm{d} A+Q_{2} \int_{A_{2}} N_{2}(x, y) \mathrm{d} A+Q_{3} \int_{A_{2}} N_{3}(x, y) \mathrm{d} A

在计算三角形单元上的数值积分，比较方便的方法是引进面积坐标。

线性三角形单元面积坐标

如下图所示，三角形单元内，任意一点面积坐标 $L_1、L_2、L_3$ 可表示为：

\begin{aligned} L_{1} & =\frac{\text { Area }_{O 23}}{\text { Area }_{123}} \\ L_{2} & =\frac{\text { Area }_{O 13}}{\text { Area }_{123}} \\ L_{3} & =\frac{\text { Area }_{O 12}}{\text { Area }_{123}} \end{aligned}

任何方向移动点 $O$ ，都将导致面积坐标 $L_1、L_2、L_3$ 随着 $x、y$ 的线性变化，符合形函数特征，且 $L_{1}+L_{2}+L_{3}=1$ ，即

\begin{array}{l}L_1=N_1(\xi,\eta)=1-\xi-\eta\\ L_2=N_2(\xi,\eta)=\xi\\ L_3=N_3(\xi,\eta)=\eta\end{array}

以上就是给大家引入的面积坐标地概念，要耐心往下看哦，方便数值积分的地方马上就来啦~

三角形单元面积坐标的数值积分（面积积分）

\int\limits_A L_i^\alpha L_j^\beta L_k^\gamma\mathrm dA=\dfrac{\alpha!\beta!\gamma!}{(\alpha+\beta+\gamma+2)!}2A

对于线性三角形单元的面积坐标阶次均为1。其中三角形面积公式，可由节点坐标计算：

A=\dfrac{1}{2}\text{det}\begin{bmatrix}1&x_1&y_1\\ 1&x_2&y_2\\ 1&x_3&y_3\end{bmatrix}

故单元降雨量又可写为：

\bar{Q}_{e}=\left(Q_{1}+Q_{2}+Q_{3}\right) \frac{A_{e}}{3}

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个二次三角形单元

我们可以构造六节点单元内部某点的试函数：

\overline{Q}(x,y)=\sum\limits_{i=1}^{6}N_i(\xi,\eta)Q_i

形函数如下式：

\left\{\begin{array}{l} N_1(\xi, \eta) \\ N_2(\xi, \eta) \\ N_3(\xi, \eta) \\ N_4(\xi, \eta) \\ N_5(\xi, \eta) \\ N_6(\xi, \eta) \end{array}\right\}=\left\{\begin{array}{c} -(1-\xi-\eta)(1-2(1-\xi \eta)) \\ 4 \xi(1-\xi-\eta) \\ -\xi(1-2 \xi) \\ 4 \xi \eta \\ -\eta(1-2 \eta) \\ 4 \eta(1-\xi-\eta) \end{array}\right\}

该区域的面积为：

A=\int_A \mathrm{~d} A=\int_0^{+1} \int_0^{1-\xi} \operatorname{det}(J(\xi, \eta)) \mathrm{d} \xi \mathrm{d} \eta=\sum_k^{n p t} W_k \operatorname{det}\left(J\left(\xi_k, \eta_k\right)\right)

Hammer数值积分

我从书里面挑了两个具有代表性的Hammer积分点位置，如下图所示，与高斯点分布不同哦~，在本题中，用到了7个积分点，具体的位置，我还没找到，直接按照书上写的数值来吧，会用就行......

使用有限元近似估计区域总降水量 | Hammer积分的图4 — Hammer积分点位置

积分点坐标即相应的权重系数可在书中查找，此处不再赘述。

总的降雨量：

\begin{aligned} Q_e & =\int_1 \bar{Q}(x, y) \mathrm{d} A \\ & =\int_0^{+1} \int_0^{1-\xi}\left(\sum_{i=1}^6 N_i(\xi, \eta) Q_i\right) \operatorname{det}[J(\xi, \eta)] \mathrm{d} \xi \mathrm{d} \eta \\ & =\sum_k^{n p t} W_k\left(\sum_{i=1}^6 N_i\left(\xi_k, \eta_k\right) Q_i\right) \operatorname{det}\left[J\left(\xi_k, \eta_k\right)\right] \end{aligned}

其中雅可比矩阵：

[J]=\left[\begin{array}{ll} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{array}\right]=\left[\begin{array}{cc} \sum_{i=1}^6 \frac{\partial N_i}{\partial \varepsilon_i} x_i & \sum_{i=1}^6 \frac{\partial N_i}{\partial \xi} y_i \\ \sum_{i=1}^6 \frac{\partial N_i}{\partial \eta} x_i & \sum_{i=1}^6 \frac{\partial N_i}{\partial \eta} y_i \end{array}\right]

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