地球物理中有限单元法.docVIP

地球物理算法技术（论文） 地球物理中的有限单元法 院 系:地球物理与信息技术院 姓 名:刘雅宁 学 号:2010120053 任课老师:张贵宾地球物理中的有限单元法 有限单元法的介绍 在地球物理理论计算中，存在着两类基本问题：正问题和反问题。给定场源的分布，求解场值的大小，这是正问题，或者称为正演问题。 地球物理正演的数值计算方法，种类很多，最常用的有：有限差分法和有限单元法。 有限单元法是50年代首先在弹性力学中发展起来的方法。主要优点是，适用于物性参数复杂分布的区域，但计算量大。随着计算机技术的发展，有限单元法在解决各个工程领域的许多数学物理问题中，得到了广泛的应用，称为一种高效、通用的计算方法。地球物理中的一些边值问题，也采用了有限单元法，解决了许多从前无法计算的地球物理问题。 有限单元法解决数学物理边值问题的基本思路和过程如下： 1、给出地球物理边值问题中的偏微分方程和边界条件（及初始条件）。这一点看起来似乎容易，但做起来并不容易，特别是边界条件的给定。只有对地球物理方法的原理和问题有深入的理解，才能给边值问题中的偏微分方程和边界条件以正确的描述。 2、将地球物理边值问题转变为有限元方程。实现这种转变的主要数学工具是变分法，用变分法得到的有限元法方程称为泛函极值问题。 3、用优先单元法解决泛函极值问题其步骤大致如下： 把研究区域剖分成有限个小单元，在每个单元上，把函数简化成线性函数、二次函数或高次函数，这称为单元上函数的插值。用简化后的函数计算每个单元上的泛函。各单元之间，通过单元间节点上的函数值相互联系起来。对各单元的泛函求和，获得整个区域上的泛函。这样，有限单元法将连续函数的泛函，离散成各单元节点上函数值得泛函。根据泛函取极值的条件，得到各节点的函数值应满足的线性代数方程组。解代数方程组，得到各节点的函数值。 有限单元法的主要优点是，适用于物性复杂分布的地球物理问题，而且，其解题过程也比较规范化。这些优点是有限单元法在地球物理中获得广泛的应用。但是，有限单元法是区域性方法，必须在全区域中进行剖分。剖分后的单元和节点数目多，最后得到的线性代数方程组很大。特别是三维问题和地球物理中常遇到的无解区域问题，需要中、大型计算机，才能完成有限单元法的计算。这是有限单元法的主要缺点。 基本步骤 用三角单元对区域进行剖分： 用三角单元对整个区域进行剖分，因为三角单元的边容易拟合地形线的形状。三角形的顶点称为节点，用节点上的离散场值来近似场值的连续分布。剖分后对单元和节点进行编号，次序号（i，j，m）按逆时针方向排列。第一类边界条件的节点上的场值是已知的，其余节点上的场值是待求的； 线性插值： 在三角单元内假定场值u是线性变化的，则u=ax+by+c，u还可表示为，其中是形函数， 单元分析： 将全区域的积分分解为单元e的积分之和 ， 继续推得单元e上的积分 其中是单元的u值列向量；是单元系数矩阵，，s，t=i，j，m，是对称矩阵。 总体合成: 将单元的场值列向量ue扩展成全体节点的场值向量u，u=(u1u2…ui…uj…um…und)，按照节点的总体序号，将单元系数矩阵中的各元放在的相应行与列的交叉位置上，其余位置的元为零，这样单元积分可写成 各单元积分相加时，只要将ke相加即可； 求变分： 通过以上四个步骤，已经将连续函数g的泛函，离散成各节点g值的多元函数：泛函的极值等于多元函数的极值，用多元函数求极值的方法，对上式求微分， 因为K是对称矩阵，有，所以，由于，所以由上式得。得到含有ND个元的ND个方程联立的线性代数方程组； 解线性代数方程组： 首先将第一类边界条件代入，通过定带宽储存的对称带型线性方程组，解方程组，得到各节点的u值，至此有限单元法的求解过程结束。 实现过程 为了验证有限单元法的有效性，我们设计一个规则形状的地下矿体，给出模型： 1、模型 密度均匀的水平半无限空间，一个均匀球体S，球体半径R=10m，剩余密度σ=1g/cm3，球心坐标（a，b）=（200，-100）。对于均匀球体来说，它与将其全部剩余质量集中在球心处的点的质量所引起的异常完全一样。若球心的埋藏深度为D，球的半径为R，剩余密度为σ，则它的剩余质量为M=(4πR3σ)/3,通过原点的任意水平剖面上则重力异常的解析解表达式为： 设测线长400m，高程变化（-200,300），地形设为一曲线：y=20*sin(0.02*x)+30，其中x为测线离原点的水平距离，y为高程，则S引起的重力异常为： 2、剖分 通过Matlab建模，得到地形曲线图，再用三角单元对划出的区域进行剖分（程序附后），剖分后对单元和节点进行编号，并将节点的xy坐标和单元节点号列表（表1和表3），分别放在XY（2，ND）和I3（3，N