1、1地球物理数值计算方法地球物理数值计算方法地球物理数值计算方法Numerical Methods in Geophysics王彦宾王彦宾 地球物理学系2008-2009学年第二学期2地球物理数值计算方法第六章 有限单元方法3有限单元法 地球内部介质，尤其是浅部，存在横向非均匀结构，包括分层、不规则形状的块体。由于解析方法不能给出这类复杂模型中的解，因此需要近似的数值求解方法。有限单元方法(Finite Element Method，FEM)是地球物理数值方法中另一种常用到的方法。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线性组合来逼
2、近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在地球物理领域，有限单元法应用于地震波场模拟、地球动力学模拟等。由于网格划分的灵活性，特别适用于非常复杂的模型，如自由地表、复杂边界模型。4有限单元法简介有限单元方法是将偏微分方程描述的连续问题进行离散求解的一种数值方法。其基本原理是：用简单的块体构造复杂的对象，
3、或将一个复杂的对象分为用以处理的小块体。例子：近似圆的面积 一个三角形的面积：圆的面积：N为三角形的个数，当N时，iiRSsin212)2sin(2121NNRSSNiiN2RSN5有限单元法简介有限单元方法包括以下基本步骤：1、根据变分原理或方程余量与权函数正交化原理，建立与微分 方程初边值问题等价的积分表达式，这是有限元法的出发 点。2、将研究对象剖分为简单的块体，即单元。3、描述每一个单元中的物理量，确定单元基函数，将各个单元 中的求解函数用单元基函数的线性组合表达式进行逼近。4、将所有单元通过节点组合到一起，按一定法则进行累加，形 成总体有限元方程。5、处理边界条件。6、采用适当的数值
4、计算方法求解根据边界条件修正的总体有限 元方程组，可求得各节点的函数值。7、计算感兴趣的单元内的函数值。6有限单元法简介为什么要用有限单元方法（有限单元方法的优势）：1、有限单元方法可以对任意形状的问题进行灵活剖分，特别适 用于非常复杂的模型，最早应用于结构力学。2、有限单元方法是工程中应用最广泛的计算及数值模拟方法。3、所需要的复杂节点生成可以利用图形界面软件（如CAD）实 现。4、有很多商用有限单元软件系统可以使用（如：ANSYS、ADINA、SMART）。7有限单元法简介有限单元方法的应用：1、力学、航空、土木工程、汽车工程2、结构分析（静力学、动力学、线形、非线性）3、热传导和流体力学
5、4、电磁学5、地质力学6、生物学8有限单元法简介有限单元方法在地球物理中的应用：1、地壳变形2、地球物理流体力学 地球动力学 地幔对流3、地电磁学4、波动传播5、强震地面运动、地震工程9有限单元法简介有限单元方法在地球物理中的应用：1、地壳变形2、地球物理流体力学 地球动力学 地幔对流3、地电磁学4、波动传播5、强震地面运动、地震工程10线性代数基础线性方程组：其中x1、x2、xn是待求变量，以上方程组写为矩阵形式：其中：bAx 11线性代数基础列向量：行向量：矩阵加减运算：其中：矩阵相乘运算：其中：其中：A：lxm矩阵，B：mxn矩阵，i=1,2,l，j=1,2,n一般地：但是12线性代数基
6、础特殊矩阵：矩阵的转置：对称矩阵：单位矩阵：并且有13线性代数基础特殊的行列式：方阵A的行列式是一个数，表示为 ，或 )det(AA14线性代数基础矩阵求逆：方阵A可逆的充分必要条件是它的行列式不为零，即 如果行列式为零，称A为奇异矩阵。矩阵求逆：对一个方阵A，如果存在一个矩阵B，使得则B是A的一个逆矩阵。当A可逆时，A的逆矩阵为：其中，A*为A的伴随矩阵，由A的代数余子式组成，0)det(AAIBAABAAA)det(11ijjiij)(MA115线性代数基础矩阵求逆：对逆矩阵有:线性方程组求解：对线性方程组 ，如果A可逆，则线性方程组求解的主要任务是求系数矩阵A的逆矩阵。常用求解技术，如：
7、高斯消元法(Gaussian elimination)（对给定的线性方程组施行初等变换，将其变成一个同解的阶梯形方程组，从而达到求解的目的。）迭代方法111)(ABABbAx bAx116线性代数基础正定矩阵：对于方阵A，如果对所有的非零的向量x，则A为正定矩阵。正定矩阵是非奇异矩阵。矩阵的微分、积分：设该矩阵对变量t的微分为：相应的积分为：0AxxT)()(tatijAdttdatdtdij)()(Adttadttij)()(A17有限单元法基础有限单元方法最早是为了解决弹性静力学问题而提出来的，因此一般按照方法的发展历史介绍其基本概念，如单元(elements)、刚度矩阵(stiffnes
8、s matrix)等等。需要的预备知识：胡克定律（Hookes Law）静力平衡原理功的概念应变能321FFFFtotFdsWDsF klijklijcW2118有限单元法基础一维问题：设有线性方程组设有向量y，。一般地，方程组两侧同乘y不改变它的解：考虑泊松方程：其中u是标量场，f是源项，一维情况下：泊松方程两侧乘以任意函数v(x)，nRyybyAx bAx)()(xfxu222x)()()()(xvxfxvxu19有限单元法基础一维问题：以上方程在整个求解区域D上求积分,为简单起见，定义区域D为0，1，则积分为区域的离散化(discretization)：为了求出近似解，需要对区域进行某种
9、形式的离散化，有限单元法中，函数值只定义在离散化后的离散点上（这一点和有限差分法类似）：DDxvxfxvxu)()()()(1010)()()()(dxxvxfdxxvxu20有限单元法基础一维问题：区域的离散化(discretization)：21有限单元法基础一维问题：有限单元法的中心思想是用选取的基函数的线性组合来近似表示函数值：其中N为区域上的离散点数，ci为系数（实数）。如果基函数 选取合理，ci等于离散点i处的真实函数值，意味着在离散点i处的基函数 ，其它点上的基函数值都为零。下面再看前面的积分方程：Niiicuu1i1i1010)()()()(dxxvxfdxxvxu22有限单元
10、法基础一维问题：对该方程的左侧作分部积分：设u在边界处的微分为零，得到：原来的方程变为：其中u为待求的未知函数值。101010)(udxvuvvdxu1010)(udxvvdxu1010fvdxvdxu23有限单元法基础一维问题：如果用近似值代替u，上式对近似值同样成立：其中：是我们利用选取的基函数展开后给出的近似函数值。V是任意选取的实数值函数，如果上式对任意函数成立，它对下式显然成立：因此对所有的基函数成立。jv1010fvdxvdxuNiiicu124有限单元法基础一维问题：将以上内容综合到一起：得到：jv1010fvdxvdxuNiiicu110101dxfdxckkniii25有限单
11、元法基础一维问题：系数ck为常数，所以对每一个基函数有：上式可以写为矩阵形式：得到解为：bi是基函数的系数，对于特殊选取的基函数，这些系数值正好是离散点i上的函数值。10101dxfdxckkiniikikigAbkiTikgbAkTikigAb126有限单元法基础一维问题：基函数：我们希望基函数具有如下性质：可以选择具有这种性质的任何函数，最简单的是如图所示的线性函数（蓝线：基函数值）i,jxxxxxjiifor 0for 1)(27有限单元法基础一维问题：基函数的梯度：为了构造刚度矩阵(stiffness matrix)，我们需要计算基函数（蓝线）的梯度（红线）：28有限单元法基础一维问题
12、：刚度矩阵(stiffness matrix)：知道以上性质的基函数后，我们可以计算矩阵形式方程中的Aij和gi：基函数是定义在区间0，1上的连续函数，10101dxfdxckkiniikikigAb10dxAkiik10dxfgkkelsewherexxxxxxxxxxxxxxxiiiiiiiiiii0forfor)(11111129有限单元法基础一维问题：刚度矩阵(stiffness matrix)：假设模型离散在等间隔节点上，设:基函数可以写为：1iiixxdxxxxelsewherexxxxxxxxxxxxxxxiiiiiiiiiii0forfor)(111111elsewheredx
13、xdxxxdxdxxxi00for10for1)(30有限单元法基础一维问题：刚度矩阵(stiffness matrix)：基函数的梯度为：1iiixxdxxxxelsewheredxxdxxdxdxxi00for/10for/1)(31有限单元法基础一维问题：刚度矩阵(stiffness matrix)：方程中Aij的计算：10dxAkiikdxdxdxdxdxdxdxdxAdxxxdxxxdxxx1111111111210111111dxdxdxdxdxdxdxdxAdxdxdxxxxdxx211 00222210222222222232有限单元法基础一维问题：刚度矩阵(stiffness
14、 matrix)：方程中Aij的计算：10dxAkiikdxdxdxdxdxdxdxdxAdxxxdxxxdxxx11 111111112102121121221AA33有限单元法基础一维问题：刚度矩阵(stiffness matrix)：方程中Aij的计算：10dxAkiik1112112111A34有限单元法基础一维问题：边界条件和源项：原始方程为：如果假设：最后变换为以下方程：加上边界条件后，设：其中u(0)和u(1)是区域0，1边界处的值。)()(xfxu10101dxfdxckkiniiNiiicu1NNiiiuucu)1()0(11235有限单元法基础一维问题：边界条件和源项：原始
15、方程为：加上边界条件后，设：最后变换为以下方程：可以写为矩阵形式：)()(xfxu10101101012)1()0(dxudxudxfdxcknkkkiniiNNiiiuucu)1()0(112kikigAbkiTikgbA36有限单元法基础一维问题：边界条件和源项：用图形形象表示（系统通过修改后的源项来考虑边界条件）：10101101012)1()0(dxudxudxfdxcknkkkiniikiTikgbA边界条件源项边界条件37有限单元法基础一维问题：数值算例（等间距离散点）：偏微分方程：区域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);边界条件：u(0)=u(1
16、)=0)()(xfxuMatlab 有限单元程序Matlab 有限差分程序38有限单元法基础一维问题：数值算例（等间距离散点）：偏微分方程：区域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);边界条件：u(0)=u(1)=0有限单元Matlab程序计算结果（蓝色）有限差分Matlab程序计算结果（红色）)()(xfxu39有限单元法基础一维问题：数值算例（等间距离散点）、边界条件不为零：偏微分方程：区域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);边界条件：u(0)=0.15 u(1)=0.05有限单元Matlab程序计算结果（蓝色）有限差分Mat
17、lab程序计算结果（红色）)()(xfxu40有限单元法基础一维问题：不等间距离散点时的刚度矩阵：10dxAkiik211021111021211211 111111111AhdxhdxhhdxdxAhhxxhxxiiiihhA11141有限单元法基础一维问题：数值算例（不等间距离散点）：偏微分方程：区域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);边界条件：u(0)=u0 u(1)=u1)()(xfxu刚度矩阵的计算程序42有限单元法基础一维问题：数值算例（不等间距离散点）：偏微分方程：区域:0,1;nx=100;dx=1/(nx-1);f(x)=(1/2);边界条件
18、：u(0)=0.15 u(1)=0.05)()(xfxu43有限单元法基础总结：1、有限单元方法的分析中，我们用一组正交基函数的线性组合来给出定义在区域D上的函数的近似值，线性组合的系数相当于某些离散节点上的函数值。2、对某些偏微分方程组，离散节点上的值可以通过求解线性方程组获得，求解线性方程组得过程包括矩阵（有时是稀疏矩阵）的求逆。3、边界条件在变换以后的方程中自然满足，这是有限单元方法的优势之一（和有限差分法相比）Niiicu1kTikigAb110101101012)1()0(dxudxudxfdxcknkkkinii44有限单元法-1D单元坐标变换：我们希望用一组基函数的线性组合来近似
19、给出区间a,b上的函数值u(x):其中i是位置xi处的离散节点（单元的边界）的编号。因为所有单元中的基函数具有相同的形式，因此为简单起见，我们将坐标系统变换到一个局部的坐标系统中：因此，单元定义在以下区间：niiicu1iiixxxx11 ,045有限单元法-1D单元1D线性单元：对于1D直线单元，没有形状的选择，但是单元的长度可以在区域内变化。我们希望在每一个单元内，函数u()用以下的线性函数来近似：离散节点定义在1,2=0,1处。我们要求：21)(ccu2122121111 uucccuuccuAuc 1101A46有限单元法-1D单元1D线性基函数：我们利用离散节点1,2处的函数值给出了
20、系数ci的值。接下来，我们利用离散节点上的值给出近似函数值：其中N1,2()是一维单元的线性基函数。)()()1()()(221121211NuNuuuuuuu47有限单元法-1D单元1D二次单元：我们希望在每一个单元内，函数u()用以下的二次函数来近似：离散节点定义在1,2,3=0,1/2,1处。我们需要：2321)(cccu32133212110.25c5.0cccuccucuAuc 242143001A48有限单元法-1D单元1D二次基函数：同样，我们可以用基函数的线性组合（系数为相应的三个离散节点上的值）近似给出函数值：注意这里我们每一个单元里用到了三个离散节点。31232221232
21、1)()2()44()231()(iiiNuuuucccu49有限单元法-1D单元1D三次基函数：用同样的过程，我们可以推导出三次基函数的形式：324323322321342321)(23)(2)(231)()(NNNNccccu50有限单元法-2D单元坐标变换：现在我们讨论二维（2D）问题的单元的几何特征和基函数，同样，我们希望在局部的坐标系统中讨论问题，以三角形为例：变换前变换后51有限单元法-2D单元坐标变换：任何一个三个角位置为 （反时针顺序）的三角形可以通过坐标变换变换为直角等腰三角形：如果=0，以上变化和1D坐标变换相同。现在我们希望用线性形式给出函数的近似：和1D的推导过程类似，
22、我们得到：3,2,1 ),(iyxPiii)()()()(1312113121yyyyyyxxxxxx321),(cccu31321211)1,0()0,1()0,0(ccuuccuucuu52有限单元法-2D单元系数：系数可以通过线性方程组求解：矩阵A为：31321211)1,0()0,1()0,0(ccuuccuucuuAuc 101011001A1101A包含1D部分53有限单元法-2D单元线性基函数：根据矩阵A和线性方程组，可以计算出系数，给出三角形单元的线性基函数：),(),(1),(321NNN54有限单元法-2D单元二次基函数：定义在三角形上的任意函数可以用二次函数近似为：在变换
23、后的坐标系里，形式为：和1D类似，这时需要另外的三 个点，共有六个点26524321),(yxyxyxyxu26524321),(ccccccu55有限单元法-2D单元二次基函数：为了计算系数，我们计算每一个节点上的值：矩阵A为：系数可以通过P点的值计算：631665432154214631342121162141414121214121cccuccccccucccucccucccucuAuc 400202444004004022400103004013000001A56有限单元法-2D单元二次基函数：基函数为：)1(4),(4),()1(4),()12(),()12(),()221)(1()
24、,(654321NNNNNN57有限单元法-2D单元二次基函数：基函数为：)1(4),(4),()1(4),()12(),()12(),()221)(1(),(654321NNNNNN58有限单元法-2D单元四边形单元：坐标变换：变换到局部坐标系统中。59有限单元法-2D单元四边形单元：线性基函数：得到矩阵A为：基函数为：4321),(ccccu1111100100110001A)1(),(),()1(),()1)(1(),(4321NNNN60有限单元法-2D单元四边形单元：二次基函数：得到矩阵A（8x8）和基函数，如：282726524321),(ccccccccu)1)(1(4),()2
25、21)(1)(1(),(51NN61有限单元法-2D单元总结：有限单元方法的基函数的求取规程：1、将坐标系统转换到局部的（单元的）坐标系统；2、对在整个单元内部取值的函数进行线性（二次、三次）函数的近似；3、利用插值条件（特殊选取的基函数保证相应点上的值为1，其它点上为0），得到系数值（离散节点上函数值的函数）4、利用求解线性方程组得到的系数推导出n个节点上的n个基函数。62有限单元法-例子声波方程：如何利用有限单元法求解与事件有关的问题？声波方程：其中v为声波速度。用和上述类似的思想，方程两侧同乘以任意选取的函数，然后在整个区域(如：0,1）上积分，做分部积分后，得：用基函数给出u的近似：f
26、uvufxuvtut2222222 10102102dxfdxuvdxujjjtNiiixtcuu1)()(63有限单元法-例子声波方程：这里的系数是时间t的函数，并且：代入积分方程，得到：写为矩阵形式：10102102dxfdxcvdxcjijiiijiitNiiitttxtcuu1222)()(M（质量矩阵）A（刚度矩阵）ggcAcMTTv2 64有限单元法-例子声波方程：对于特殊选取的基函数，系数值为离散节点上的函数值u。求解时间相关项：时间微分采用有限差分近似：得到tk+1时刻的解为：我们知道如何计算刚度矩阵A，如何计算质量矩阵M？gcAcMTTv2 gcAMkTkkkTvdtccc2
27、2112122112)()(kkkTTkdtvcccAgMc65有限单元法-例子声波方程：M为：基函数为：10102102dxfdxcvdxcjijiiijiit10dxMjiijelsewherehxhxxhhxxiiiii00for10for1)(1166有限单元法-例子声波方程：M为：基函数为：对角线元素Mii,i=2,n-1：10dxMjiijelsewherehxhxxhhxxiiiii00for10for1)(1133 1100221101iihhiiiiiihhdxhxdxhxdxMii67有限单元法-例子声波方程：68有限单元法-例子声波方程：规则网格（蓝色：FD，红色：FEM）：69有限单元法总结（方法的比较）：1、有限差分方法：2、有限单元方法：3、Fourier伪谱方法适用于时间有关的偏微分方程（地震波传播、地球物理流体动力学、电磁方程、穿地雷达等）特点：灵活、概念简单、容易编程、容易并行、规则网格、显式方法适用于静态的、时间有关的偏微分方程（地震波传播、地球物理流体动力学等所有问题）特点：隐式方法、需要计算矩阵求逆、不规则网格（灵活）、方法完善、复杂的算法和编程、工程问题正交基函数，可以看作有限差分的特例空间微分算子的精度是谱精度，适用于和时间有关的偏微分方程（地震波传播、穿地雷达等地球物理问题）特点：显式方法、规则网格、计算效率高、精度高
