基于二次场的二维大地电磁有限元法数值模拟
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第 35 卷第 8 期 2007 年 8 月
同 济 大 学 学 报 (自 然 科 学 版) JOURNAL OF TON GJ I UN IV ERSITY(NA TURAL SCIENCE)
Vol. 35 No . 8 Aug. 2007
基于二次场的二维大地电磁有限元法数值模拟
刘小军 ,王家林 ,于 鹏
(5)
9 Ez s 9y -
9 Eys 9z
=-
iωμ0 Hxs
(6)
式 (1) ～ (6) 中 : E 为电场值 ; H 为磁场值 ;下标 p 表
示一次场 ; s 表示二次场 ;ε为介电常数 ;μ0 为自由 空间的磁导率 ; Δσ为岩石电导率σ与不存在异常
体的电导率σ3 之差 ;ω为角频率.
将式 (1) 和 (2) 变形后代入式 (3) ,可得 TE 模式
(同济大学 海洋地质国家重点实验室 ,上海 200092)
摘要 : 引入基于计算二次场响应的方法 ,将其用于二维大地电磁数值模拟. 导出了二维大地电磁二次场的表达式 , 在使用有限元法解偏微分方程时 ,利用二次场偏微分方程与总场偏微分方程的相似性 ,直接采用总场法中的系数 矩阵作为二次场的合成矩阵. 将根据此原理所编制的二次场有限元法模拟软件 ,对几个典型模型进行了试算 ,并与 传统总场模拟法作了比较分析. 结果表明 ,直接计算二次场的结果精度更高 ,更接近真实解.
Exs.
从而提高了模拟精度.
图 1 层状模型解析法与总场法视电阻率曲线 Fig. 1 Comparison of apparent resistivity by using
analytical solution and total f ield FEM
3 模型试验与分析
为验证基于二次场有限元法 (简称二次场法) 的 计算效果 ,首先对总场有限元法 (简称总场法) 模拟 结果与解析解进行了比较 ,说明了二次场法相对总 场法计算精度高的主要原因. 对两个理论模型进行 了试算 ,并与总场法模拟计算结果作了比较 ,重点对 比分析两种方法受网格密度的影响. 3. 1 精度比较
图 1 为四层地电剖面模型 (图中所镶嵌的小图) 及测点 A 处视电阻率ρA 随周期变化曲线 ,两条曲 线分别为解析解与总场法的计算结果. 从中可以看 出 ,在短周期部分 ,两种方法求得的视电阻率值基本 相同 ,当周期 T > 500 s 时 ,两条曲线出现较明显的 偏差 ,低频段总场法的计算精度明显降低 ,最高相对 误差接近 5 %. 说明总场法模拟在低频段的计算效 果不太理想. 因此 ,本文提出将一次场从总场中分离 出来 ,采用精度较高的解析法求一次场 ,用有限元法 求二次场 ,相对提高了总场的精度. 图 2 中镶嵌的小图为层状介质中一个方形异常 体 ,异常体电阻率为 1 Ω·m ,曲线显示了模型中测点 A ′处采用二 次 场 法 和 总 场 法 计 算 的 视 电 阻 率 随 周 期变化的规律. 可以看出 ,在高频段 (低周期) ,两组 视电阻率曲线基本相似 ; 而在低频段 ( T > 400 s) , 两组曲线出现分离 ,总场法求出的视电阻率值相对 二次场的值有较明显的减小 ,相对误差最高达 3 %. 结合图 1 的结论 ,分析表明 ,二次场有限元法模拟的 结果更合乎实际. 主要原因是由于模拟计算中一次 场采用相对精确的解析法求得 ,相对降低了低频段 有限元法计算二次场时精度稍低带来的误差比重 ,
在 TE 极化模式时 ,远离异常体的外边界 Γ 上 不受异常体的影响 ,因此取边界处异常场 ExsΓ = 0. 由此可得 TE 极化模式时二次场的边值问题等价于 下列变分问题 :
∫∫ 1
F( Exs) = 2
1 Ω iωμ0 (
Δ
Exs) 2 - (σ -
∫∫ iωε) E2xs dΩ - ΩΔσExp ExsdΩ (9)
TE 模式
9 Exs 9z
=-
iωμ0 Hys
(1)
9 Exs 9y
= iωμ0 Hz s
(2)
9 Hzs 9y -
9 Hys 9z
=
(σ -
iωε) Exs +ΔσExp
(3)
TM 模式 9 Hxs 9z
=
(σ -
iωε)
Eys +ΔσEyp
(4)
9 Hxs 9y = -
(σ - iωε) Ez s - ΔσEzp
Leabharlann Baidu
下二次场的赫姆霍兹方程
9 1 9 Exs
9 1 9 Exs
9 y iωμ0 9 y + 9 z iωμ0 9 z -
(σ - iωε) Exs = ΔσExp
(7)
同理得 TM 模式下的方程
9
1 9 Hxs
9 y σ- iωε 9 y
9
1 9 Hxs
+ 9 z σ- iωε 9 z
- iωμ0 Hxs =
第 35 卷
电磁场模拟时有一些固有缺点[3 ] . 国内很多学者提 出一些改善计算精度和提高计算速度的方法 ,黄俊 革等人[4 ]提出了计算异常电位的有限元激发极化 数值模拟 ,采用异常电位法计算节点电位 ,保证了计 算精度 ;欧阳联华等人[5 ] 通过改变插值函数和剖分 方式来提高解的精确度 ,并且详细讨论了选择不同 插值函数和不同网格密度计算时对电磁场值的影 响 ;陈小斌等人[6 ] 提出直接迭代算法 ,直接在计算 网格的基本结构上 ,由泛函求极值过程中直接得出 完全类似于有限差分法中的显式迭代格式 ,这样避 免了刚度矩阵并提高了计算速度.
已有的有限元法模拟大地电磁场的改进方法 ,大 都是通过改变插值方式来提高计算速度和精度 ,而且 是基于计算总场的响应 ,受网格密度影响较大 ,且在 提高低频段的计算值精度等方面效果不太理想[2 - 3]. 针对上面问题 ,本文在前人研究的基础上 ,提出基于 二次场的大地电磁数值模拟方法 ,即将大地电磁场分 为两部分 :一部分为源 ,在一维均匀水平层状介质背 景中产生的电磁场 (简称为一次场) ,可由解析法导 出 ;另一部分为将地下异常体看作散射源 ,其产生的 二次场可通过麦克斯韦方程组导出其偏微分方程 ,采 用有限元法求解出二次场值. 通过理论模型试验发 现 ,二次场模拟方法比总场模拟效果更好 ,受网格密 度影响相对较低 ,而且较好地解决了低频段的场值计 算精度不高的问题 ,说明本思路是正确可行的.
L IU Xiaoj u n , W A N G J iali n , Y U Peng
( State Key Laboratory of Marine Geology , Tongji University , Shanghai 200092 , China)
Abstract : This paper present s an efficient met hod named secondary finite element met hod ( FEM) for two2dimensional magnetotelluric numerical modeling , which provides solution for secondary variations in t he field. To deduce t he rationale of t he new procedure ,t he total field is divided into primary and secondary field. Only t he secondary field equations which need to be solved for t he primary field can be obtained easily. The equations of t he secondary field can be solved effortlessly by finite element met hod when t he boundary conditions are simple. The effect of resolution is analysed to st udy t he response of synt hetic model by t he secondary field and a comparison is conducted between t he result s of secondary FEM and t he total field FEM. It is concluded t hat t he new arit hmetic of numerical modeling can be used to derive more precise result s in calculating t he secondary field and is almost not affected by t he density of mesh. It is expected to be used widely for t he inversion of magnetotelluric data. Key words :magnetotelluric ; finite element met hod ; total field ; secondary field ; numerical modeling
1 二次场方程的表达式
由电磁场理论 ,在非均匀地电结构中传播的大 地电磁场可分解为一次 (电或磁) 场与二次场之和. 其中一次场为无异常体时的场 ,二次场为异常体产 生的场. 二维地质构造情况下 ,取 x 轴平行于构造 走向 , y 轴垂直于走向方向 , z 轴垂直地面向下. 参 照文献[ 7 ]中二次电磁场方程的表达式 ,类推可得二 维大地电磁正演中 TE (电场分量平行于走向的模 式) 和 TM (磁场分量平行于走向的模式) 两种极化 模式下 ,二次场的麦克斯韦方程的表达式
Δσσ-iωiμωε0 Hxp
+
9 9z
Δσ σ - iωε
Eyp
(8)
式 (7) 和 (8) 中 ,所有一次场值都可由一维均匀水
平层状介质的解析法求得. 这样解两式就可得两种极
化模式下平行走向的二次场 Exs和 Hxs , 辅助二次场
Hys , Hzs和 Eys , Ezs可由式 (1) , (2) 和 (4) , (5) 求得 , 采
收稿日期 : 2006 - 06 - 06 基金项目 : 中科院知识创新工程资助项目 ( KZCX1 - SW - 18 - 01) 作者简介 : 刘小军 (1979 - ) ,男 ,湖北潜江人 ,博士生. E2mail :joneslxj @163. com
1 11 4
同 济 大 学 学 报 (自 然 科 学 版)
关键词 : 大地电磁 ; 有限元法 ; 总场 ; 二次场 ; 数值模拟 中图分类号 : P 631 文献标识码 : A 文章编号 : 0253 - 374X(2007) 08 - 1113 - 05
Secondary Field2Ba sed Two2Dimensional Magnetotelluric Numerical Modeling by Finite Element Method
ExsΓ = 0 δF ( Exs) = 0
第 8 期
刘小军 ,等 :基于二次场的二维大地电磁有限元法数值模拟
111 5
式中 ,Ω 为单元面积. 采用有限元法解式 (9) ,首先用矩形单元对区域
进行剖分 ,剖分时细划界面和物性变化较大的地方 , 其余则用相对较大的网格. 文中算例均采用矩形单 元 、双二次插值[11 - 12 ] . 单元积分及总体系数矩阵合 成参见文献 [ 11 ] ,可得二次场的线性方程组 Ks Exs = Ss. 式中 : Ks 为 M ×M 阶对称系数矩阵 , 其各项 元素与模型电导率及网格剖分有关 ; Exs为网格节点 上的二次场 ; Ss 为列向量 ,通过对单元内电导率和 一次电场积的面积积分后 ,再进行总体合成求得 ; M 为节点个数. 本文采用 L U (下三角矩阵 、上三角 矩阵) 分解法解方程组 ,可得到各节点的二次电场值
有限元法作为一种解偏微分方程的数值方法 , 已 广 泛 应 用 于 各 种 地 球 物 理 方 法 的 正 演 模 拟. Coggon[1 ]最早提出将有限元法应用电磁场与激发 极化场的模拟 ,并讨论了二维带地形的激发极化场
的模拟效果 ,随后 Wannamaker 等人[2 ] 将其应用到 大地电磁场的模拟中 ,并提出了带地形的模拟方法. 但有限元法是用非无限小的单元和有限网格来模拟 无限空间 ,本身就是一种近似. 因此 ,在对二维大地
用差分法求二次场的偏导数 ,两场之和即为总场.
2 基于二次场的有限元解法
从式 (7) 和 (8) 可以看出 ,两种极化模式下的二 次场偏微分方程与总场方程[7 - 10 ]相似 ,不同之处只 是在方程右端多了一个与一次场值和电导率差有关 的源项. 因而 ,在采用有限元法进行数值模拟时 ,总 场法的有限元系数矩阵 K 可直接用于求二次场的 方程中. 限于篇幅 ,只讨论 TE 极化模式情况.
同 济 大 学 学 报 (自 然 科 学 版) JOURNAL OF TON GJ I UN IV ERSITY(NA TURAL SCIENCE)
Vol. 35 No . 8 Aug. 2007
基于二次场的二维大地电磁有限元法数值模拟
刘小军 ,王家林 ,于 鹏
(5)
9 Ez s 9y -
9 Eys 9z
=-
iωμ0 Hxs
(6)
式 (1) ～ (6) 中 : E 为电场值 ; H 为磁场值 ;下标 p 表
示一次场 ; s 表示二次场 ;ε为介电常数 ;μ0 为自由 空间的磁导率 ; Δσ为岩石电导率σ与不存在异常
体的电导率σ3 之差 ;ω为角频率.
将式 (1) 和 (2) 变形后代入式 (3) ,可得 TE 模式
(同济大学 海洋地质国家重点实验室 ,上海 200092)
摘要 : 引入基于计算二次场响应的方法 ,将其用于二维大地电磁数值模拟. 导出了二维大地电磁二次场的表达式 , 在使用有限元法解偏微分方程时 ,利用二次场偏微分方程与总场偏微分方程的相似性 ,直接采用总场法中的系数 矩阵作为二次场的合成矩阵. 将根据此原理所编制的二次场有限元法模拟软件 ,对几个典型模型进行了试算 ,并与 传统总场模拟法作了比较分析. 结果表明 ,直接计算二次场的结果精度更高 ,更接近真实解.
Exs.
从而提高了模拟精度.
图 1 层状模型解析法与总场法视电阻率曲线 Fig. 1 Comparison of apparent resistivity by using
analytical solution and total f ield FEM
3 模型试验与分析
为验证基于二次场有限元法 (简称二次场法) 的 计算效果 ,首先对总场有限元法 (简称总场法) 模拟 结果与解析解进行了比较 ,说明了二次场法相对总 场法计算精度高的主要原因. 对两个理论模型进行 了试算 ,并与总场法模拟计算结果作了比较 ,重点对 比分析两种方法受网格密度的影响. 3. 1 精度比较
图 1 为四层地电剖面模型 (图中所镶嵌的小图) 及测点 A 处视电阻率ρA 随周期变化曲线 ,两条曲 线分别为解析解与总场法的计算结果. 从中可以看 出 ,在短周期部分 ,两种方法求得的视电阻率值基本 相同 ,当周期 T > 500 s 时 ,两条曲线出现较明显的 偏差 ,低频段总场法的计算精度明显降低 ,最高相对 误差接近 5 %. 说明总场法模拟在低频段的计算效 果不太理想. 因此 ,本文提出将一次场从总场中分离 出来 ,采用精度较高的解析法求一次场 ,用有限元法 求二次场 ,相对提高了总场的精度. 图 2 中镶嵌的小图为层状介质中一个方形异常 体 ,异常体电阻率为 1 Ω·m ,曲线显示了模型中测点 A ′处采用二 次 场 法 和 总 场 法 计 算 的 视 电 阻 率 随 周 期变化的规律. 可以看出 ,在高频段 (低周期) ,两组 视电阻率曲线基本相似 ; 而在低频段 ( T > 400 s) , 两组曲线出现分离 ,总场法求出的视电阻率值相对 二次场的值有较明显的减小 ,相对误差最高达 3 %. 结合图 1 的结论 ,分析表明 ,二次场有限元法模拟的 结果更合乎实际. 主要原因是由于模拟计算中一次 场采用相对精确的解析法求得 ,相对降低了低频段 有限元法计算二次场时精度稍低带来的误差比重 ,
在 TE 极化模式时 ,远离异常体的外边界 Γ 上 不受异常体的影响 ,因此取边界处异常场 ExsΓ = 0. 由此可得 TE 极化模式时二次场的边值问题等价于 下列变分问题 :
∫∫ 1
F( Exs) = 2
1 Ω iωμ0 (
Δ
Exs) 2 - (σ -
∫∫ iωε) E2xs dΩ - ΩΔσExp ExsdΩ (9)
TE 模式
9 Exs 9z
=-
iωμ0 Hys
(1)
9 Exs 9y
= iωμ0 Hz s
(2)
9 Hzs 9y -
9 Hys 9z
=
(σ -
iωε) Exs +ΔσExp
(3)
TM 模式 9 Hxs 9z
=
(σ -
iωε)
Eys +ΔσEyp
(4)
9 Hxs 9y = -
(σ - iωε) Ez s - ΔσEzp
Leabharlann Baidu
下二次场的赫姆霍兹方程
9 1 9 Exs
9 1 9 Exs
9 y iωμ0 9 y + 9 z iωμ0 9 z -
(σ - iωε) Exs = ΔσExp
(7)
同理得 TM 模式下的方程
9
1 9 Hxs
9 y σ- iωε 9 y
9
1 9 Hxs
+ 9 z σ- iωε 9 z
- iωμ0 Hxs =
第 35 卷
电磁场模拟时有一些固有缺点[3 ] . 国内很多学者提 出一些改善计算精度和提高计算速度的方法 ,黄俊 革等人[4 ]提出了计算异常电位的有限元激发极化 数值模拟 ,采用异常电位法计算节点电位 ,保证了计 算精度 ;欧阳联华等人[5 ] 通过改变插值函数和剖分 方式来提高解的精确度 ,并且详细讨论了选择不同 插值函数和不同网格密度计算时对电磁场值的影 响 ;陈小斌等人[6 ] 提出直接迭代算法 ,直接在计算 网格的基本结构上 ,由泛函求极值过程中直接得出 完全类似于有限差分法中的显式迭代格式 ,这样避 免了刚度矩阵并提高了计算速度.
已有的有限元法模拟大地电磁场的改进方法 ,大 都是通过改变插值方式来提高计算速度和精度 ,而且 是基于计算总场的响应 ,受网格密度影响较大 ,且在 提高低频段的计算值精度等方面效果不太理想[2 - 3]. 针对上面问题 ,本文在前人研究的基础上 ,提出基于 二次场的大地电磁数值模拟方法 ,即将大地电磁场分 为两部分 :一部分为源 ,在一维均匀水平层状介质背 景中产生的电磁场 (简称为一次场) ,可由解析法导 出 ;另一部分为将地下异常体看作散射源 ,其产生的 二次场可通过麦克斯韦方程组导出其偏微分方程 ,采 用有限元法求解出二次场值. 通过理论模型试验发 现 ,二次场模拟方法比总场模拟效果更好 ,受网格密 度影响相对较低 ,而且较好地解决了低频段的场值计 算精度不高的问题 ,说明本思路是正确可行的.
L IU Xiaoj u n , W A N G J iali n , Y U Peng
( State Key Laboratory of Marine Geology , Tongji University , Shanghai 200092 , China)
Abstract : This paper present s an efficient met hod named secondary finite element met hod ( FEM) for two2dimensional magnetotelluric numerical modeling , which provides solution for secondary variations in t he field. To deduce t he rationale of t he new procedure ,t he total field is divided into primary and secondary field. Only t he secondary field equations which need to be solved for t he primary field can be obtained easily. The equations of t he secondary field can be solved effortlessly by finite element met hod when t he boundary conditions are simple. The effect of resolution is analysed to st udy t he response of synt hetic model by t he secondary field and a comparison is conducted between t he result s of secondary FEM and t he total field FEM. It is concluded t hat t he new arit hmetic of numerical modeling can be used to derive more precise result s in calculating t he secondary field and is almost not affected by t he density of mesh. It is expected to be used widely for t he inversion of magnetotelluric data. Key words :magnetotelluric ; finite element met hod ; total field ; secondary field ; numerical modeling
1 二次场方程的表达式
由电磁场理论 ,在非均匀地电结构中传播的大 地电磁场可分解为一次 (电或磁) 场与二次场之和. 其中一次场为无异常体时的场 ,二次场为异常体产 生的场. 二维地质构造情况下 ,取 x 轴平行于构造 走向 , y 轴垂直于走向方向 , z 轴垂直地面向下. 参 照文献[ 7 ]中二次电磁场方程的表达式 ,类推可得二 维大地电磁正演中 TE (电场分量平行于走向的模 式) 和 TM (磁场分量平行于走向的模式) 两种极化 模式下 ,二次场的麦克斯韦方程的表达式
Δσσ-iωiμωε0 Hxp
+
9 9z
Δσ σ - iωε
Eyp
(8)
式 (7) 和 (8) 中 ,所有一次场值都可由一维均匀水
平层状介质的解析法求得. 这样解两式就可得两种极
化模式下平行走向的二次场 Exs和 Hxs , 辅助二次场
Hys , Hzs和 Eys , Ezs可由式 (1) , (2) 和 (4) , (5) 求得 , 采
收稿日期 : 2006 - 06 - 06 基金项目 : 中科院知识创新工程资助项目 ( KZCX1 - SW - 18 - 01) 作者简介 : 刘小军 (1979 - ) ,男 ,湖北潜江人 ,博士生. E2mail :joneslxj @163. com
1 11 4
同 济 大 学 学 报 (自 然 科 学 版)
关键词 : 大地电磁 ; 有限元法 ; 总场 ; 二次场 ; 数值模拟 中图分类号 : P 631 文献标识码 : A 文章编号 : 0253 - 374X(2007) 08 - 1113 - 05
Secondary Field2Ba sed Two2Dimensional Magnetotelluric Numerical Modeling by Finite Element Method
ExsΓ = 0 δF ( Exs) = 0
第 8 期
刘小军 ,等 :基于二次场的二维大地电磁有限元法数值模拟
111 5
式中 ,Ω 为单元面积. 采用有限元法解式 (9) ,首先用矩形单元对区域
进行剖分 ,剖分时细划界面和物性变化较大的地方 , 其余则用相对较大的网格. 文中算例均采用矩形单 元 、双二次插值[11 - 12 ] . 单元积分及总体系数矩阵合 成参见文献 [ 11 ] ,可得二次场的线性方程组 Ks Exs = Ss. 式中 : Ks 为 M ×M 阶对称系数矩阵 , 其各项 元素与模型电导率及网格剖分有关 ; Exs为网格节点 上的二次场 ; Ss 为列向量 ,通过对单元内电导率和 一次电场积的面积积分后 ,再进行总体合成求得 ; M 为节点个数. 本文采用 L U (下三角矩阵 、上三角 矩阵) 分解法解方程组 ,可得到各节点的二次电场值
有限元法作为一种解偏微分方程的数值方法 , 已 广 泛 应 用 于 各 种 地 球 物 理 方 法 的 正 演 模 拟. Coggon[1 ]最早提出将有限元法应用电磁场与激发 极化场的模拟 ,并讨论了二维带地形的激发极化场
的模拟效果 ,随后 Wannamaker 等人[2 ] 将其应用到 大地电磁场的模拟中 ,并提出了带地形的模拟方法. 但有限元法是用非无限小的单元和有限网格来模拟 无限空间 ,本身就是一种近似. 因此 ,在对二维大地
用差分法求二次场的偏导数 ,两场之和即为总场.
2 基于二次场的有限元解法
从式 (7) 和 (8) 可以看出 ,两种极化模式下的二 次场偏微分方程与总场方程[7 - 10 ]相似 ,不同之处只 是在方程右端多了一个与一次场值和电导率差有关 的源项. 因而 ,在采用有限元法进行数值模拟时 ,总 场法的有限元系数矩阵 K 可直接用于求二次场的 方程中. 限于篇幅 ,只讨论 TE 极化模式情况.