二维电导率各向异性介质电磁感应效应的有限元算法.docx

二维电导率各向异性介质电磁感应效应的有限元算法李予国摘要二维电性各向异性大地电磁场通过有限元方法计算获得模型由包含各向异性块的层状介质构造组成每个块或者层可能是由 3×3 的电导率张量限定正演问题可以归结为两个耦合的关于垂直-平行于场分量的 和 它们由有限元法数值地解出线性有限元的解系是用预条件共轭梯度法求取的之后，地xEH表垂直-平行的场分量 和 ，利用样条插值通过对 和 进行数值微分得到y xEH二维有限元算法通过对比二维有限差分的解证明是正确的三个类型的模型的大地电磁响应用来证明各向异性的影响：水平，垂直，倾斜各向异性第四个模型用来模拟在剪切和俯冲带各向异性的影响这些模型的响应模拟了电阻率曲线在长周期和具有明显的主对角线元素的张量主抗存在室的分离，正如以前观测到的1、简介 近些年来，越来越多的关注投入到了对电性各向异性的电磁感应方面上的研究，尤其是为了试图完全理解对更长周期的大地电磁的观测数据加拿大地区的大地电磁测量显示了电性的各向异性在下地壳和上地幔（Kellett 等，1992；Mareschal 等 1995） 大的来自围绕德国深钻点（KTB）各向异性大地电磁曲线解释为上地壳高度的电各向异性（Eisel 和 Haak 1999） 。

Rasmussen (1988)使用一个深地壳层内的各向异性模型来解释瑞典南部的大地电磁横断面数据层状构造的各向异性效应首先被 O’Brien 和 Morrison (1967)研究Reddy 和 Rankin (1975)首先开始研究二维各向异性模型，他们只考虑水平各向异性的工作最近，Osella和 Martinelli (1993)计算了一些圆滑不规则边界有一定的主轴旋转角的大地电磁响应Schmucker (1994)提出一个计算不均匀薄的覆盖体在层状半空间的电磁感应效应，其中可以含有一个或者两个各向异性的电导层用有限元方法，Pek 和 Verner (1997) 还有 Weidelt (1996) 模拟了广义二维、三维各向异性构造，分别具有任意的主轴旋转角本文中，二维模型被重新研究，但是使用有限元（FE）的方法首先，我们详细描述有限元程序的数值近似之后，我们展示有限元程序模拟不同的简单测试模型的响应我们的结果和 Pek& Verner (1997)的有限差分程序进行了对比最后，我们计算三种类型的各向异性介质的大地电磁响应：水平，垂直，以及倾斜各向异性我们最后以模拟一个地质剪切带或者板块俯冲带的地质背景下的模型结束。

本文这样的构造是为了展示一个有限元策略去实现二维各向异性介质的电磁感应尽管有限差分对于这类问题是有效的，有限元解仍然是需要的，因为有限元和有限差分都有自己的特殊的优势还可以彼此互相核对另外，有限元方法可以模拟非矩形线性构造的真实地球构造2、边界值问题考虑到图 1 中的二维模型阴影异常区嵌入到一个简单的层状构造中，这个构造包括 n 层而且底部的层延伸到无穷远为了简单，异常区域没有显示内边界作为一个定向独立的电导率 ，它至少和周边的第 j 层的电导率张量 不同然后我们的程序允许更一般的模型，异常区域可以细分成不同的不均匀块j而且可以接触不同的层异常区域延伸到 y 正向到无穷远，例如模型可以出现到一个新的正常构造在这种状态下算法排除源的感应模型在走向方向 x 方向是不变的感应电磁场在 x 方向也是不变的，尽管如此场的矢量也有三个分量对于极化的磁源平行或者垂直走向，由于各向异性因此，通常区分 TE 和 TM 模式的各向同行构造就变得不合理了由此就有结果，当主磁场矢量垂直轴向，电场将有主场方向的分量，导致电荷产生于界面（除了倾斜各向异性的例子） 假设一个时谐变量 ，在准静态近似的电磁场的情况下iwte(1), oEiHE其中 是自由空间的磁导率，且oxyxzzxyz为电性电导率张量。

它是对称的，当旋转到主轴方向（x',y',z'）上去，可以表示为'''0'xyz图 1.本文中考虑的二维各向异性模型在特别的二维例子中，方程（1）消退为常数电导率的均匀块 , (2), 3; yzoxxyozEiwHi (4), 5, (6). yzxxyxzxyxyyzzzzHE 7很明显如果垂直-平行分量 和 得到了，剩下的分量 和 ， 和 可由 和 的空间xEHyEzyHzxE导数得到方程（2）-（7）可以得到两个关于分量 和 的二阶偏微分方程：x200yzzyyxzyx10, (8)()()() , 9-, (-)/,()/xxxxxxxCABiwyzHHiDADBD 其 中 xzyzC1BA从方程（8）和（9 ）很明显看出各向异性和独立的垂直-平行分量 和 通过一阶偏微分方程耦合xEH在一起因此，在这种情况下不存在独立的 TE、TM 模式对于广义各向异性情况因此，这些方程必须同时利用 和 求解xEH之后叙述了不同的空间形式的各向异性的模拟结果。

如图 1 一个均匀的异常区域具有各向同性的正常结构的围岩， 一个或者两个平行于参考主轴（x，y，z）坐标系的主轴被考虑了进去2.1 水平各向异性当 ，主轴 Z'是垂直的，另外两个主轴 x'和 y'在水平面（x,y）具有相对于 x 轴的垂直角0xzy方程（8）和（ 9）就消退为220 01(), (10)1() ()0, xyxxy xyxxzyHEEiwzHiwEyz感应方程变得比方程（8）和（9）更加简单，但是垂直-平行分量 和 仍然耦合在一起通过关于xHz 的一阶偏导所以方程（10）和（11）必须联立在一起去为了得到 和 2.2 倾斜各向异性当 ，电导率张量主轴 x'是水平的在走向方向，另外两个主轴 y'和 z'在垂直面（y,z）具0xyz有相对于 y 轴的倾斜角 现在方程（8）和（9 ）就解耦为两个独立的模式20, (12)() 3xxEiwH方程（12）对应于 E 极化模式可以通过二维各向同行程序解得- 只需要电导率替换为 作为标量电导率x然而，TM 模式的表达式子仍然比较复杂2.3 垂直各向异性如果 ，电导率张量的三个主轴真好和坐标系（x，y，z）是一致的。

如果再有0xyz，这种结果关于 z 轴轴对称，水平方向的电导率 完全不同于垂直方向的电导率 xh hv这样（8）和（9）仍然是耦合在一起的，我们有 200, (14)11()(), 5xhxxxvEiwHiwHyz尽管 B 极化模式感应方程比方程（9）更加简单，它不能通过二维各向同性程序解得，因为电导率在水平方向和垂直方向是不同的如果两个方向的电导率相同，方程（15）消退为各向同性的例子之后的边界条件应用：模型的外边界条件，设定为狄里克雷边界条件，由对应的一维层状模型的左右边界条件在模型的顶部和底部由左右边界的值进行线性插值得到（Pek 和 Verner 1997） 在内边界上，切向电场和磁场的分量， 和 必须连续根据图 1， 和 为tEtHtEtHcosin,.tyz重构方程式子之间的关系至关重要设 n 为不均匀区域的外侧的单位矢量，且有p, (16)xyxzAEeB和 来自方程（8）和（9） ， 和 为 y 轴和 z 轴的单位向量相应的这样，实质性的偏导数后，yez表达式为 pn (17)xtHE和 01 (8)xthiwn为发现的切向分量。

方程（17）是自动满足的，对自后的偏导是很重要的结合这些表达式基本的积分-微分方程（23）和（24）可以得到浓缩的形式，这些形式可以直接正演数值模拟对待更多细节看李（2000） 3 有限元方法由方程（8）和（9 ）组成的数值近视问题，基于有限元方法垂直-水平分量 和 跟电导构造相xEH切，所以处处连续于是，在有限元中假设所有分量穿过单元都是满足连续这种近似建立在模型 域，整天包围二维不均匀体，区域延伸的足够远使得异常场在边界处很小为了避免绝缘的空气层在地面上产生奇异性，我们假设空气层具有很小的电导率，但是非零，实际上小于 在我们的计算中数120Sm值计算结果显示电导率到 301Sm加权余量方法用来获得从偏微分方程（8）和（9）得到的积分方程方程(8)乘以一个关于电场的任意小的的变量 ，在模型域 Ω 里的积分：xE201( )0 (19)xxxHCABdiwyz上述方程中，第一项包含二阶偏导数，可以用格林函数简化2 uvdvvn其中 指示模型域 的边界之后方程（19）可以写成等价形式0 01 10 (2) xxxxxxHEECEdAEdBddiwyziwn 类似地，方程（9）乘以一个磁场分量任意小的变量 ,之后用高斯公式改变一下xuvdnuvd加上方程（16）这样得到积分方程0() ()0 (21)xxxxxHHiwHpdpnd 这里我们仍然用这个公式 (2) xxHn模型域 可被分割为矩形或者是三角形单元。

矩形单元已被李（2000）年描述过以下的部分，展示三角形单元公式（20）和（21）分解到每个单元带有指示标 e=1,2,…ne0 01 1() (23)ene nene nxxxx txi i i iHEdCEdAEBdHEdwyzw  01 111() 4 neexxx xtxe eHwp 其中 表示为三角形单元 e 的边界在形成方程（23）的时候，我们利用方程（18）时替换被积函数为e.用一种简单的方法方程（24）由方程（17）获得 txHE现在边界条件必须被利用在内部边界切向电场和磁场（ ）是被要求的因为每个边界在集,xEH成的过程中被遍历两次且方向相反，则在单元边界的线积分的总和为零在外边界，因为设置的狄利克雷边界条件，电场和磁场的变量 和 等于 0.因此鲜鸡粉也等于 0.这样以来，方程（23）和（24）xExH最终消退为 01()0 (25)ne nenexxxxi i iEdCdAEBdwyz01 11() 6 nexx xe eHwpH 根据有限元的线性近似方法，我们假设在每个三角形单元内， 和 都为关于 y,z 的线性函数：xE3 31 1(,), (,) (27) xixiEyzNyzN其中 和 为第 i 个顶点的电场和磁场对应的坐标为（ ） ，i=1,2,3 为三角形 （图 2） ，iiH,iyz13为线性基函数。

它们由一下定义：i1(),12,3 (8)2iiiiNaybzc其中 121(), 29为 三 角 形 的 面 积 , （ ）33123221213,, (30)1,, 2azbycyz方程（25）和（26）在一个单元内的面积分用方程（27）- （32）估算所有单元的积分都可以组装成两个线性方程系详细内容见附录 A联立这些方程系，我们可以最终写下方程以矩阵的形式：0 (3)KU其中12,xEkH和 是 阶方阵（ 是结点的数目） ， 和 不对称的方阵，满足 因此矩阵 K12dnd12k12Tk是一个 2 的对称方阵它是稀疏的复元素U 是 2 列向量，它是包含电场和磁场的未知数方程dn（33）代入外边界条件，在内部的点上 和 可以用共轭梯度方法解出来xE图 2 三角形单元切面共轭梯度。