MT2D大地电磁有限元正演程序--单元分析

MT2D有限元单元分析

泛函公式

$$J(U)=\sum _V{\int }_e\frac{1}{2}\sigma (\nabla U{)}^{2}dV+\sum _V{\int }_e\frac{1}{2}\sigma k^2{U}^{2}dV-I\delta (A)U$$
对于每一个单元，需要对上式的第一项和第二项的单元积分。最重要的是求$\nabla U$和$U$的积分计算。

三角单元的形函数

U的单元分析

假定在三角单元内的电位$U$是线性变化的，即
$$U=N_i{U}_i+N_j{U}_j+N_m{U}_m$$
其中，$N_n$为形函数，$n=i,j,m$
$$N_n=\frac{1}{2\Delta }\left(a_n+b_{n}x+c_{n}y\right)$$
由于形函数仅与空间坐标有关系，与其他因素无关，其中$a_n,b_n$公式如下：
$$\{\begin{array}{l}{a}_i=x_j{y}_m-y_j{x}_m,b_i=y_j-y_m,c_i=x_j-x_m\\ a_j=x_m{y}_i-y_m{x}_j,b_j=y_m-y_i,c_j=x_m-x_i\\ a_m=x_i{y}_j-y_i{x}_j,b_m=y_i-y_j,c_m=x_i-x_j\\ \Delta =\frac{1}{2}\left(b_i{c}_j-b_j{c}_i\right)\end{array}$$
利用matlab简单实现计算单元的算法

function fem
    % 假设有一个三角单元，其三个点的坐标分别如下：
    p1 = [.5 .5];
    p2 = [0 0];
    p3 = [1 0];
    hs=scatter([p1(1);p2(1);p3(1)],[p1(2);p2(2);p3(2)]);
    hold on
    hp=patch([p1(1);p2(1);p3(1)],[p1(2);p2(2);p3(2)],[.7 .7 .7]);

    a(1) = p2(1)*p3(2)-p2(2)*p3(1);
    a(2) = p3(1)*p1(2)-p3(2)*p1(1);
    a(3) = p1(1)*p2(2)-p1(2)*p2(1);
    b(1) = p2(2)-p3(2);
    b(2) = p3(2)-p1(2);
    b(3) = p1(2)-p2(2);
    c(1) = p3(1)-p2(1);
    c(2) = p1(1)-p2(1);
    c(3) = p2(1)-p1(1);
    delta = 0.5*(b(1)*c(2)-b(2)*c(1));
    if delta <= 0
        message('出错了，妈的！！')
    end

    N1 = femN(p1,delta,a,b,c)
    N2 = femN(p2,delta,a,b,c)
    N3 = femN(p3,delta,a,b,c)
end

function N = femN(p,delta,a,b,c)
    factor = 1.0/(2.0*delta);
    for i =1:3
        N(i) = (a(i)+b(i)*p(1)+c(i)*p(1))*factor;
    end
end

$\nabla U$的单元分析

因为，
$$\nabla U=\frac{\partial U}{\partial x}\vec{ı}+\frac{\partial U}{\partial y}\vec{ȷ}$$
所以有，
$$(\nabla U{)}^2={\left(\frac{\partial U}{\partial x}\right)}^2+{\left(\frac{\partial U}{\partial y}\right)}^2$$
分开来看，
$$\frac{\partial U}{\partial x}=\frac{\partial }{\partial x}\left(\sum _{n=i,j,m}{N}_n{U}_n\right)=\sum _{n=i,j,m}\frac{\partial N_n}{\partial x}{U}_n={\left(\frac{\partial \vec{N}}{\partial x}\right)}^T{\vec{U}}_e={\vec{U}}_e^T\left(\frac{\partial \vec{N}}{\partial x}\right)$$
$${\left(\frac{\partial U}{\partial x}\right)}^2={\vec{U}}_e^T\left(\frac{\partial \vec{N}}{\partial x}\right){\left(\frac{\partial \vec{N}}{\partial x}\right)}^T{\vec{U}}_e$$
那么则有
$${\left(\frac{\partial \vec{N}}{\partial x}\right)}^T=\left(\frac{\partial N_i}{\partial x},\frac{\partial N_j}{\partial x},\frac{\partial N_m}{\partial x}\right)=\frac{1}{2\Delta }\left(b_i,b_j,b_m\right)$$
同理，
$${\left(\frac{\partial \mathrm{U}}{\partial \mathrm{y}}\right)}^2={\vec{U}}_e^T\left(\frac{\partial \vec{N}}{\partial y}\right){\left(\frac{\partial \vec{N}}{\partial \mathrm{y}}\right)}^{\mathrm{T}}{\vec{U}}_e$$
其中，
$${\left(\frac{\partial \vec{N}}{\partial y}\right)}^T=\left(\frac{\partial N_i}{\partial y},\frac{\partial N_j}{\partial y},\frac{\partial N_m}{\partial y}\right)=\frac{1}{2\Delta }\left(c_i,c_j,c_m\right)$$
通过上述表达式可求出：
$${\int }_e\frac{1}{2}\sigma (\nabla U{)}^{2}dV=\frac{1}{2}\sigma {\vec{U}}_e^T{\int }_e\left(\frac{\partial \vec{N}}{\partial x}\right){\left(\frac{\partial \vec{N}}{\partial x}\right)}^T\left(\frac{\partial \vec{N}}{\partial y}\right){\left(\frac{\partial \vec{N}}{\partial y}\right)}^{T}dV{\vec{U}}_e=\frac{1}{2}{\vec{U}}_e^T{k}_{1e}{\vec{U}}_e$$
其中：
$$k_{1e}=\frac{1}{4\Delta }\left[\begin{array}{cc}{b}_i& c_i\\ b_j& c_j\\ b_m& c_m\end{array}\right]\left[\begin{array}{ccc}{b}_i& b_j& b_m\\ c_i& c_j& c_m\end{array}\right]=\frac{1}{4\Delta }{b}_s{c}_t,s,t=i,j,m$$

对于第二项，
$${\int }_e\frac{1}{2}\sigma k^2{U}^{2}dV=\frac{1}{2}\sigma k^2{\vec{U}}_e^T{\int }_e\vec{N}{\vec{N}}^{T}dV{\vec{U}}_e=\frac{1}{2}{\vec{U}}_e^T{k}_{2e}{\vec{U}}_e$$
又由，
$${\int }_{\mathrm{e}}{N}_i^a{N}_j^b{N}_m^{c}dV=2\Delta \frac{a!b!c!}{(a+b+c+2)!}$$
其中，
$$k_{2e}=\frac{1}{12}\left[\begin{array}{lll}2& 1& 1\\ 1& 2& 1\\ 1& 1& 2\end{array}\right]$$

function [GG NN] = GGNN(delta,b,c)
    factor = 0.25/delta;
    for i=1:3
        for j=1:3
            GG(i,j)=factor*(b(i)*b(j)+c(i)*c(j));
            if i==j
                NN(i,j)= delta/6.0;
            else
                NN(i,j)=delta/12.0;
            end
        end
    end
end

以上为三角单元的单元分析的原理和matlab的代码。

MT2D中的单元分析代码如下，可以通过上述的matlab代码和公式去理解

#ifndef  _FEM_H
#define  _FEM_H

// C++ includes
#include <cassert>
#include "tri.h"

class FEM
{
 public: 
   FEM(Tri& tri);   % 构造函数
   ~FEM(); 			% 解析函数
   void init(); 
   void N(Point p, double n[3]);  % 求形函数 
   void Grad_N(Point p, Point grad_n[3]); 
   void Grad_N_dot_N(Matrix3D& GG,  Matrix3D& NN); %求解形函数的二阶导数相乘，和N*N  
 
 public:
   Tri*                          tri;
   double                        s;
   double                        a[3], b[3], c[3];   
};


inline
FEM::FEM(Tri& tri_):
  tri  (&tri_)
{
  this->init();
}

inline
FEM::~FEM()
{
}

inline
void FEM::init()
{
  Point& v0 = *tri->get_node(0);
  Point& v1 = *tri->get_node(1);
  Point& v2 = *tri->get_node(2);
  a[0]= v1(0)*v2(1) - v1(1)*v2(0);
  a[1]= v2(0)*v0(1) - v2(1)*v0(0);
  a[2]= v0(0)*v1(1) - v0(1)*v1(0);
  b[0]= v1(1)-v2(1);
  b[1]= v2(1)-v0(1);
  b[2]= v0(1)-v1(1);
  c[0]= v2(0)-v1(0);
  c[1]= v0(0)-v2(0);
  c[2]= v1(0)-v0(0);
  
  double size = 0.5*(b[0]*c[1]-b[1]*c[0]);
  this-> s    = size;
  assert(size>0.0);
  
  return;
}

inline
void FEM::N(Point p, double N[3])
{
  assert(this->s>0.);
  const double factor= 1.0/(2.0*this->s);
  for(int i=0; i<3; i++) N[i]= (a[i]+b[i]*p(0)+c[i]*p(1))*factor;
  return;
}



inline
void FEM::Grad_N(Point p, Point grad_n[3])
{
  assert(s>0.);
  const double factor= 1.0/(2.0*this->s);

  for(int i=0; i<3; i++) {
    grad_n[i](0) = b[i]*factor;
    grad_n[i](1) = c[i]*factor;    
  }
  return;
}


inline
void FEM::Grad_N_dot_N(Matrix3D& GG,  Matrix3D& NN)
{
  GG.setZero(); NN.setZero();

  for(int i=0; i<3; i++)
   for(int j=0; j<3; j++) 
     GG(i,j)= (0.25/s)*(b[i]*b[j]+c[i]*c[j]);
   
  for(int i=0; i<3; i++)
   for(int j=0; j<3; j++) {
     if(i==j) NN(i,j)= s/6.0;
     else NN(i,j)= s/12.0;
   }   
  return;
}


#endif // _FEM_H