MT2D有限元单元分析

泛函公式

J(U)=∑V∫e12σ(∇U)2dV+∑V∫e12σk2U2dV−Iδ(A)UJ(U)=\sum_{V} \int_{e} \frac{1}{2} \sigma(\nabla U)^{2} d V+\sum_{V} \int_{e} \frac{1}{2} \sigma k^{2} U^{2} d V-I \delta(A) U J(U)=V∑∫e21σ(∇U)2dV+V∑∫e21σk2U2dV−Iδ(A)U
对于每一个单元，需要对上式的第一项和第二项的单元积分。最重要的是求∇U\nabla U∇U和UUU的积分计算。

三角单元的形函数

U的单元分析

假定在三角单元内的电位UUU是线性变化的，即
U=NiUi+NjUj+NmUmU=N_iU_i+N_jU_j+N_mU_mU=NiUi+NjUj+NmUm
其中，NnN_nNn为形函数，n=i,j,mn=i,j,mn=i,j,m
Nn=12Δ(an+bnx+cny)N_{n}=\frac{1}{2 \Delta}\left(a_{n}+b_{n} x+c_{n}y\right) Nn=2Δ1(an+bnx+cny)
由于形函数仅与空间坐标有关系，与其他因素无关，其中an,bna_{n},b_{n}an,bn公式如下：
{ai=xjym−yjxm,bi=yj−ym,ci=xj−xmaj=xmyi−ymxj,bj=ym−yi,cj=xm−xiam=xiyj−yixj,bm=yi−yj,cm=xi−xjΔ=12(bicj−bjci)\left\{\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}\right. ⎩⎪⎪⎨⎪⎪⎧ai=xjym−yjxm,bi=yj−ym,ci=xj−xmaj=xmyi−ymxj,bj=ym−yi,cj=xm−xiam=xiyj−yixj,bm=yi−yj,cm=xi−xjΔ=21(bicj−bjci)
利用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 onhp=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 <= 0message('出错了，妈的！！')endN1 = femN(p1,delta,a,b,c)N2 = femN(p2,delta,a,b,c)N3 = femN(p3,delta,a,b,c)
endfunction N = femN(p,delta,a,b,c)factor = 1.0/(2.0*delta);for i =1:3N(i) = (a(i)+b(i)*p(1)+c(i)*p(1))*factor;end
end

∇U\nabla U∇U的单元分析

因为，
∇U=∂U∂xı⃗+∂U∂yȷ⃗\nabla U=\frac{\partial U}{\partial x} \vec{\imath}+\frac{\partial U}{\partial y} \vec{\jmath} ∇U=∂x∂Uı+∂y∂Uȷ
所以有，
(∇U)2=(∂U∂x)2+(∂U∂y)2(\nabla U)^{2}=\left(\frac{\partial U}{\partial x}\right)^{2}+\left(\frac{\partial U}{\partial y}\right)^{2} (∇U)2=(∂x∂U)2+(∂y∂U)2
分开来看，
∂U∂x=∂∂x(∑n=i,j,mNnUn)=∑n=i,j,m∂Nn∂xUn=(∂N⃗∂x)TU⃗e=U⃗eT(∂N⃗∂x)\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) ∂x∂U=∂x∂(n=i,j,m∑NnUn)=n=i,j,m∑∂x∂NnUn=(∂x∂N)TUe=UeT(∂x∂N)
(∂U∂x)2=U⃗eT(∂N⃗∂x)(∂N⃗∂x)TU⃗e\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} (∂x∂U)2=UeT(∂x∂N)(∂x∂N)TUe
那么则有
(∂N⃗∂x)T=(∂Ni∂x,∂Nj∂x,∂Nm∂x)=12Δ(bi,bj,bm)\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) (∂x∂N)T=(∂x∂Ni,∂x∂Nj,∂x∂Nm)=2Δ1(bi,bj,bm)
同理，
(∂U∂y)2=U⃗eT(∂N⃗∂y)(∂N⃗∂y)TU⃗e\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} (∂y∂U)2=UeT(∂y∂N)(∂y∂N)TUe
其中，
(∂N⃗∂y)T=(∂Ni∂y,∂Nj∂y,∂Nm∂y)=12Δ(ci,cj,cm)\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) (∂y∂N)T=(∂y∂Ni,∂y∂Nj,∂y∂Nm)=2Δ1(ci,cj,cm)
通过上述表达式可求出：
∫e12σ(∇U)2dV=12σU⃗eT∫e(∂N⃗∂x)(∂N⃗∂x)T(∂N⃗∂y)(∂N⃗∂y)TdVU⃗e=12U⃗eTk1eU⃗e\int_{e} \frac{1}{2} \sigma(\nabla U)^{2} d V=\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} d V \vec{U}_{e}=\frac{1}{2} \vec{U}_{e}^{T} k_{1 e} \vec{U}_{e} ∫e21σ(∇U)2dV=21σUeT∫e(∂x∂N)(∂x∂N)T(∂y∂N)(∂y∂N)TdVUe=21UeTk1eUe
其中：
k1e=14Δ[bicibjcjbmcm][bibjbmcicjcm]=14Δbsct,s,t=i,j,mk_{1 e}=\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}, \quad s, t=i, j, m k1e=4Δ1⎣⎡bibjbmcicjcm⎦⎤[bicibjcjbmcm]=4Δ1bsct,s,t=i,j,m

对于第二项，
∫e12σk2U2dV=12σk2U⃗eT∫eN⃗N⃗TdVU⃗e=12U⃗eTk2eU⃗e\int_{e} \frac{1}{2} \sigma k^{2} U^{2} d V=\frac{1}{2} \sigma k^{2} \vec{U}_{e}^{T} \int_{e} \vec{N} \vec{N}^{T} d V \vec{U}_{e}=\frac{1}{2} \vec{U}_{e}^{T} k_{2 e} \vec{U}_{e} ∫e21σk2U2dV=21σk2UeT∫eNNTdVUe=21UeTk2eUe
又由，
∫eNiaNjbNmcdV=2Δa!b!c!(a+b+c+2)!\int_{\mathrm{e}} N_{i}^{a} N_{j}^{b} N_{m}^{c} d V=2 \Delta \frac{a ! b ! c !}{(a+b+c+2) !} ∫eNiaNjbNmcdV=2Δ(a+b+c+2)!a!b!c!
其中，
k2e=112[211121112]k_{2 e}=\frac{1}{12}\left[\begin{array}{lll}{2} & {1} & {1} \\ {1} & {2} & {1} \\ {1} & {1} & {2}\end{array}\right] k2e=121⎣⎡211121112⎦⎤

function [GG NN] = GGNN(delta,b,c)factor = 0.25/delta;for i=1:3for j=1:3GG(i,j)=factor*(b(i)*b(j)+c(i)*c(j));if i==jNN(i,j)= delta/6.0;elseNN(i,j)=delta/12.0;endendend
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

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

大地电磁正演程序MT2D主程序分析
0 相关背景信息见大地电磁二维正演程序--详细介绍(https://blog.csdn.net/spvfly/article/details/92795423). 程序的框架如下图 1 主程序MT2 ...
视电阻率正演matlab,大地电磁测深一维正演——地电学实验报告分析.doc
大地电磁测深一维正演--地电学实验报告分析实验报告课程名称: 地电学课题名称: 大地电磁层状模型数值模拟实验专业: 地球物理学姓名: xx 班级: 06xxxx 完成日期: 2016 ...
matlab 地震正演程序,seismic-forward 地球物理勘探中，基础的地震波正演模拟程序，包括五个（高阶中心差分 matlab 266万源代码下载- www.pudn.com...
文件名称: seismic-forward下载收藏√ [ 5 4 3 2 1 ] 开发工具: Visual C++ 文件大小: 2985 KB 上传时间: 2014-08-06 下载次数 ...
二维大地电磁有限元数值模拟矩形+线性插值
% smt2djxfwd01.m % 二维大地电磁有限元数值模拟矩形+线性插值 clear close all clc format short g % % 读频率数据 %% [nfreq,freq ...
球重力异常matlab程序,球体重力异常正演程序介绍.docx
<应用地球物理学>课程作业基于MATLAB的球体重力正演程序实验报告 1一程序简介本程序基于MATLAB软件的GUI模块编写,旨在实现球体重力正演结果的可视化分析.MATLAB是一个高级 ...
大地电磁二维正演程序--详细介绍
0 相关背景 1.开发平台vs+intel MKL + Eigen 2.正演方法为有限单元法 3.程序分为两部分: 网格剖分---使用Triangle程序进行正演计算---使用任政勇教授开源的MT2 ...
matlab地震波褶积正演程序直达波有频散
clc; clear; h=300;%界面埋深300m V1=2800;%上层速度 dx=5; %道间距 V2=3500;%下层速度 trace_num=20; fs=1000; N=12000; T ...
TEM一维正演matlab,大地电磁学chp3一维正演.ppt
大地电磁学Geo-electromagnetismMagnetotellurics(MT) 地球物理专业用成都理工大学 Mao Lifeng 2010年9月1日第三章均匀水平层状介质大地电磁测深 ...
用matlab作地震波vsp图,《多层介质vsp正演方法研究》开题报告.doc
<多层介质vsp正演方法研究>开题报告.doc 本科毕业设计(论文)开题报告题目多层介质VSP正演方法研究学生姓名院(系)油气资源学院专业班级指导教师完成时间2009年3月15日要求1.开 ...

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