计算地球物理课件-第3章-电磁场数值模拟.pptx

1、计算地球物理地球物理与信息工程学院 物探系周 辉2013年第三章 电磁场数值模拟 内容提要第一节 电磁响应的交错网格有限差分模拟第二节 大地电磁场有限元模拟=+=00ttDH JBEDB 在用时间域有限差分法（FDTDfinite-difference time-domain method）求解麦克斯韦方程时，只求解两个旋度方程，而两个散度方程自然得到满足。第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟E 电场强度H 磁场强度D 电感应强度B 磁感应强度J 电流密度q 自由电荷=0 DEBHJE地质雷达探测时，通常只记录与测线方向垂直的水平电场分地质雷达探测时，通常只记录与测线方向垂

2、直的水平电场分量，量，并假设并假设地下介质为二维的地下介质为二维的，故考虑，故考虑仅有一个电场分量的麦克仅有一个电场分量的麦克斯韦方程组。斯韦方程组。二维非均匀介质中的麦克斯韦方程组为二维非均匀介质中的麦克斯韦方程组为 1yxzyyEHHEItzx,1yxEHtz,1yzEHtx,第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟zx0,x z,x z二维介质 Yee 网格及电磁场空间配置.1 2 I 1 2 K z i k x Ey Hz Hx 第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟=+=CSSCSttdddtddtDH JBEHlDsJsElBs1yxzyyEHHE

3、Itzx 1 21 21 2(,)111(,)(,)22nnnxxxHi kHi kHi kzz 1 21 21 2(,)111(,)(,)22nnnzzzHi kHikHikxx 1 21(,)1(,)(,)nynnyyEi kEi kE i ktt 第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟利用利用 1 211(,)(,)(,)2nnnyyyEi kEi kEi k 得得 11 21 2121 21 213411(,)(,)(,)(,)2211(,)(,),22nnnnyyxxnnnzzyEi kAEi kAHi kHi kAHikHikA Ii k 第一节 电磁响应的交错网

4、格有限差分（时间域有限差分）模拟1yxEHtz,1yzEHtx 在在(1 2)nt时刻磁场分量的差分格式时刻磁场分量的差分格式 1 21 2511(,)(,)(,1)(,)22nnnnxxyyHi kHi kAEi kE i k 1 21 2611(,)(,)(1,)(,)22nnnnzzyyHikHikAE ikE i k 1.1 1.1 求解麦克斯韦方程组的时间域有限差分法求解麦克斯韦方程组的时间域有限差分法 第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟11 21 2121 21 213411(,)(,)(,)(,)2211(,)(,),22nnnnyyxxnnnzzyEi k

5、AE i kAHi kHi kAHikHikA Ii k1 21 2511(,)(,)(,1)(,)22nnnnxxyyHi kHi kAE i kE i k 1 21 2611(,)(,)(1,)(,)22nnnnzzyyHikHikAE ikE i k 求解步骤：求解步骤：（1 1）用用磁场的更新式磁场的更新式求求(1 2)nt时刻的磁场；时刻的磁场；（2 2）用用电场的更新式电场的更新式求求(1)nt时刻的电场。如此循环直至最大时刻。时刻的电场。如此循环直至最大时刻。1.1 1.1 求解麦克斯韦方程组的时间域有限差分法求解麦克斯韦方程组的时间域有限差分法 第一节 电磁响应的交错网格有限差

6、分（时间域有限差分）模拟设计设计 PMLPML的原理是：当内部介质和边界介质的波阻抗一样，的原理是：当内部介质和边界介质的波阻抗一样，即波阻抗匹配时，电磁波入射到边界时，电磁波只沿原方向传即波阻抗匹配时，电磁波入射到边界时，电磁波只沿原方向传播而不被反射回来。当然，必须要求电磁波在边界中快速衰减，播而不被反射回来。当然，必须要求电磁波在边界中快速衰减，当电磁波传播到最外侧的边界时，电磁波的能量衰减到当电磁波传播到最外侧的边界时，电磁波的能量衰减到 6060 分分贝以下，这时即使有边界反射，能量也贝以下，这时即使有边界反射，能量也极其极其微弱，对需要计算微弱，对需要计算波场的影响可以忽略。波场的

7、影响可以忽略。PMLPML 吸收效果最好吸收效果最好。PMLPML 已经在声波、弹性波模拟中得到已经在声波、弹性波模拟中得到广泛应用。广泛应用。第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟以无吸收介质的电磁波传播为例说明以无吸收介质的电磁波传播为例说明 PMLPML 吸收边界。吸收边界。yxzEHHtzx （1 1）yxEHtz，yzEHtx （2 2）将将yE分解为分解为x x方向传播的部分和方向传播的部分和z z方向传播的部分，即方向传播的部分，即 yyxyzEEE 在边界中，电导率和磁导率分别记为在边界中，电导率和磁导率分别记为*,。第一节 电磁响应的交错网格有限差分（时间域

8、有限差分）模拟（1 1）式写为式写为 yxzxyxEHEtx yzxzyzEHEtz （2 2）式写为）式写为 *yxyzxzxEEHHtz *yxyzzxzEEHHtx 第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟*00zzxx*00zzxx*0*0,0,0zzxx *0,0,0zzxx*0,00,0 xxzz*0,00,0 xxzz 内部区域内部区域*0,0 x z*0,00,0 xxzz PML 第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟yxzEHHtzxyxEHtzyzEHtx yxzxyxEHEtx yzxzyzEHEtz*yxyzxzxEEHHtz*yxy

9、zzxzEEHHtx 根据阻抗匹配条件，可以得到如下关系根据阻抗匹配条件，可以得到如下关系 *xx，*yy 边界内的电导率是可变的。以左右边界为例，边界内的电导率是可变的。以左右边界为例，x的分布为的分布为 maxmax,0,11,1MxML xxxL xL xL xxNLxxNLxxNLxL x L L为为 PMLPML 的层数，的层数，N N为横向剖分的网格数，为横向剖分的网格数，M M为一常数，一般取为为一常数，一般取为3 3。第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟L L为为 PMLPML 的层数，的层数，N N为横向剖分的网格数，为横向剖分的网格数，M M为一常数，一

10、般取为为一常数，一般取为3 3。max为为 PMLPML 最外侧的电导率，它由边界的反射系数的大小决定。最外侧的电导率，它由边界的反射系数的大小决定。L L一般取一般取 8 8-1010 左右，当然，左右，当然，L L越大吸收效果越好，但计算效率就越大吸收效果越好，但计算效率就越低。越低。第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟第一节 电磁响应的交错网格有限差分（时间域有限差分）模拟波场的等时间剖面（波场的等时间剖面（snap shot）讨论二维大地电磁场的有限元模拟。大地电磁(MT)测深是一种以天然电磁场为场源，探测地下电性垂向变化的勘探方法，在构造单元划分、基底起伏形态研究

11、、沉积盆地电性分层、油气田探查、地热探查和地震预报中，得到了广泛的应用。大地电磁反演是当前地球物理反演研究的热门领域。正演是反演的基础，而大地电磁的正演，目前只对一维分层均匀介质有解析解。求解一维连续介质、二维和三维介质的大地电磁场的解析解是十分困难的，除少数简单情况外，一般都不能得到解析解。在实际工作中，常常遇到二维和三维介质，这时，必须依靠数值方法计算大地电磁场，最有效的数值方法是有限元法。第二节 大地电磁场有限元模拟 角频率为角频率为（时间因子为时间因子为i te）的定态电磁场的方程：）的定态电磁场的方程：（1 1）（2 2）当平面电磁波入射地面时，地下介质中的电磁波总以平面波当平面电磁

12、波入射地面时，地下介质中的电磁波总以平面波形式，几乎垂直地向下传播。将形式，几乎垂直地向下传播。将(1)和和(2)式按分量展开，并考虑到式按分量展开，并考虑到0z ，得两个独立的方程组，并以，得两个独立的方程组，并以 z 分量为准，分别命名为分量为准，分别命名为 E型和型和 H 型波。型波。第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟=+ttBEEHE 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 zE和和zH应满足的偏微分方程：应满足的偏微分方程：可统一表示成可统一表示成 对于对于E E型：型：对于对于 H H 型：型：第二节第二节 大地电磁场有限元模拟大地电磁场有限元模

13、拟 为了求解方程为了求解方程，还必须给出边界条件。，还必须给出边界条件。外边界条件外边界条件 E E 型波型波取图所示的研究区域。上边界取图所示的研究区域。上边界 ABAB 离地面足够远，使异离地面足够远，使异常场在常场在 ABAB 上为零，以该处的上为零，以该处的u u为为 1 1 单位：单位：第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 zxy 外边界条件外边界条件 E E 型波型波 下边界下边界 CDCD 以下为均质岩石，局部不均匀体的异常场在以下为均质岩石，局部不均匀体的异常场在 CDCD 上上为零，电磁波在为零，电磁波在 CDCD 以下的传播方程为以下的传播方程为 （-y

14、方向传播）方向传播），0u是常数，是常数，2ki。对低频电磁波，对低频电磁波，ki。对对u u求导求导ukuy。因为。因为 CDCD 处处yn，所以，所以 CDCD 处的边界处的边界条件：条件：第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 zxy 外边界条件外边界条件 E E 型波型波 取左右边界取左右边界 AD,BCAD,BC 离局部不均匀体足够远，电磁场在离局部不均匀体足够远，电磁场在 ADAD，BCBC 上左右对称，其上的边界条件是上左右对称，其上的边界条件是 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 zxy 外边界条件外边界条件 H H 型波型波取所示的研究区域

15、。取所示的研究区域。上边界上边界 ABAB 直接取在地面上，直接取在地面上，并以该处的并以该处的u u为为 1 1 单位：单位：这是因为对大地电磁波来说，有这是因为对大地电磁波来说，有1，有有 所以，空气中的所以，空气中的zH近似为常量。这样，可取地面的近似为常量。这样，可取地面的u u=1=1。其它条其它条件件同同E E波。波。第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 zxy 内边界条件内边界条件 E E 型波型波由于由于zE沿着介质分界面，根据电场切向分量连续性可沿着介质分界面，根据电场切向分量连续性可知，在两种介质的分界面知，在两种介质的分界面1上，有上，有 12uu 由图

16、，磁场的切向分量由图，磁场的切向分量 根据介质分界面上磁场切向分量的连续性可知，介质两侧的根据介质分界面上磁场切向分量的连续性可知，介质两侧的1zEin是连续的，即是连续的，即 zEn n 与与 x x 的夹角为的夹角为90，方向余弦为，方向余弦为sin，n 与与 y y 的夹角为的夹角为，方向余弦，方向余弦为为cos，HtHt 括号内为方向导数。括号内为方向导数。y1 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 内边界条件内边界条件 H H 型波型波由于由于 HzHz 沿着介质分界面，根据磁场切向分量连续性可沿着介质分界面，根据磁场切向分量连续性可知，在两种介质的分界面知，在两种

17、介质的分界面1上，有上，有 12uu 由图，电场切向分量由图，电场切向分量 根据介质分界面上电场切向分量的连续性可知，介质两侧的根据介质分界面上电场切向分量的连续性可知，介质两侧的1zHin是连续的，即是连续的，即 1 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 zHzxy综合以上讨论，边值问题归纳为综合以上讨论，边值问题归纳为 2 1 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 1211112211220,nSnuuuu vdv uuvvudv u duvvudv u dv uvvduvvudv udv uduvvuvdd AAAASAdV 为为简简单单起起见见，假假

18、定定在在区区域域中中只只有有一一个个不不均均匀匀体体。用用满满足足边边界界条条件件 v v=0 0（在在 A AB B 上上）的的任任意意 v v 乘乘微微分分方方程程两两侧侧并并积积分分 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 1111121211211211uu vduuvdvduvvudnnuuuvdvdvdnnnuvvuduvduvvudn 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 外法线方向相反110ABBCCDDACDuuu vdvduvvudnuuuuvdvdvdvdnnnnuvvuduvduvvudnvuuv d ,0,CDnA ddkuvd A

19、第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 与与定定解解问问题题对对应应的的由由虚虚功功原原理理得得到到的的变变分分问问题题：01,CDvuuv dkuvduAB 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 泛泛函函变变分分的的求求法法：0J uJ uu 例例 1 1：求求泛泛函函 120J uu dx的的变变分分。12001100022J uuudxuuu dxu udx 例例 2 2：求求泛泛函函 1220J uuudx的的变变分分。12200J uuuuudx 1022,J uu uuu dxuu 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 第二节第二

20、节 大地电磁场有限元模拟大地电磁场有限元模拟 121111221122110nnuuuuudu uu uuu du u du uuu du u du u du uuu du udu udu uuu duun 1122ududu uuu dn 用用 u 的变分的变分 u 乘微分方程两侧乘微分方程两侧并积分并积分 vvv AAAnA dd A第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 1100CDCDuuuududu uuudnku udu uuudu uuudku ud 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 222200111102222CDCDCDuuu u dk

21、u udu uu uu u dku udxxyyuuudku dxy 故故由由位位能能原原理理所所得得泛泛函函可可定定义义为为 22211122201,CDF uuudku dF uuAB 为简单起见，假定在区域中只有一个不均匀体为简单起见，假定在区域中只有一个不均匀体。构造泛函构造泛函 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 其变分为其变分为 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 的第二和四项为零。的第二和四项为零。第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 是区域的外边界，是区域的外边界，I是内边界，是内边界，12,n n分别是分别是12,的外

22、法向，的外法向，方向正好相反。根据内边界条件，方向正好相反。根据内边界条件，有有 即内边界条件不出现在即内边界条件不出现在 I u的变分中的变分中，自动满足，自动满足。2 1 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 因此因此 uI uudn 将下列边界条件将下列边界条件 代入得代入得 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 所以，边值问题所以，边值问题 与下列变分问题等价与下列变分问题等价 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 双线性插值在求解二维大地电磁变分双线性插值在求解二维大地电磁变分问问题中的应用。用矩形题中的应用。用矩形单元对区域进行

23、剖分，并给每个单元的电阻率单元对区域进行剖分，并给每个单元的电阻率赋值。矩形单元赋值。矩形单元的宽度为的宽度为a a，高度为，高度为b b。单元的四个节点的编号单元的四个节点的编号按规律编排按规律编排。在单。在单元内进行双线性插值。将变分问题中的区域积分分解为各单元积元内进行双线性插值。将变分问题中的区域积分分解为各单元积分之和：分之和：右侧最后一项积分只对右侧最后一项积分只对 CDCD 边界上的单元进行。边界上的单元进行。第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 1234567891011122221111222Teeeeeuuuddxdyxy u K u 1eijKK 41i

24、iiuN u，（iN为形函数为形函数，1114iiiN）0,1 1,0 1 2 3 4 第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 4jjiiijejjiieNNNNKdxxyyNNNNabd dxxyy。单元积分单元积分 2421224211242362124eabK。第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 当单元的当单元的_12边落在无穷远边界上时，线积分边落在无穷远边界上时，线积分 其中其中3eijKK由下式给出由下式给出 320010000000061002ekbK。第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 1eK，2eK和和3eK都是都是 4

25、 4 阶矩阵，将阶矩阵，将其其扩展成全体节点的矩阵扩展成全体节点的矩阵1eK，2eK，3eK。然后将各单元的扩展矩阵相加，。然后将各单元的扩展矩阵相加，F（u）变成变成 其中其中123eeeCDKKKK是总体系数矩阵。是总体系数矩阵。对对 F（u）式求式求变分，得变分，得 由由u的任意性，得的任意性，得 Ku0。解线性代数方程组前，将解线性代数方程组前，将 ABAB 线上的边界值代入。线上的边界值代入。u 代表各节代表各节点的点的zH或或zE。至此，有限单元法求解。至此，有限单元法求解 u 的过程全部结束。的过程全部结束。用用zH或或zE可进一步求得视电阻率。可进一步求得视电阻率。第二节第二节 大地电磁场有限元模拟大地电磁场有限元模拟 第三节第三节 大地电磁场矢量有限元模拟大地电磁场矢量有限元模拟 eyl 1234123441,eeeiiix yEEN 123411,2211,22eeyyececeeeeyyeeececxxeeeexxllyyyyllllxxxxllNiNiNjNj exl xy