前面讲过，许多重要的电法问题，如点源二维电阻率法和激发极化法，二维大地电磁问题，二维的甚低频电磁场及线源问题等有关的电位或电磁场均可以用（8.3.26）式或（9.3.1）式二维赫姆霍兹方程来描述，因此与此相对应的有限元法在电法勘探数字模拟的实际应用上，具有重要的意义。
我们知道，由（9.3.1）式二维赫姆霍兹方程
地球物理数据处理教程
在第三类边界条件
地球物理数据处理教程
的解与满足泛函
地球物理数据处理教程
极小的函数相对应。
式中：u为目标函数；D为研究区域；Γ为D的边界。
9.5.1 系数矩阵的导出
当用有限元法求解泛函式（9.5.3）极小时，其步骤和上节中与拉普拉斯方程对应的二维有限元法相同。在D域的离散中仍然采用最简单的三角形剖分，划分的方法和原则均同上节所述，这里，我们研究某个三角形单元e，其顶点分别按逆时针编号，记为i，j，m，相应的坐标为（xi，zi）、（xj，zj）、（xm，zm），其函数值为ui、uj、um，设u=u（x，z）在单元内线性变化，同样可得到（9.4.4）式由三角形顶点函数值表示的单元e中的函数值为
u（x，z）=Ni（x，z）ui+Nj（x，z）uj+Nm（x，z）um （9.5.4）
式中
地球物理数据处理教程
地球物理数据处理教程
称为基函数，其中
αi=xjzm-xmzj，bi=zj-zm，ci=xj-zj
αj=zmzi-xizm，bj=zm-zi，cj=xi-zm
αm=xizj-xjzi，bm=zi-zj，cm=xj-zi
Δe=
为了计算三角形中的泛函值，必须先计算一阶导数，利用（9.5.4）式可写出
地球物理数据处理教程
其中
地球物理数据处理教程
将（9.5.3）式的积分域D分成各小单元三角形之和，即
地球物理数据处理教程
这里Je为一个小三角形单元e内的泛函
地球物理数据处理教程
设α、β、γ均为常数，可以提到积分号外，如果单元e在D域内部，则第二项的线积分为零，若三角形单元e有一个边在边界上，则要计算线积分部分。为了研究有一般性意义，不妨设单元e的一个边在D域的边界。
其目的在于求目标函数 u，使泛函 J（u）取极小，按照极小必要条件，应有下式成立：
地球物理数据处理教程
式中l0为D域中所有节点总和。为了求出上式的具体形式，先考虑（9.5.5）式对ul的一阶导数，为计算方便，不妨先设l=i。
地球物理数据处理教程
由前面引出的
地球物理数据处理教程
同理，
地球物理数据处理教程
利用（9.5.4）式，可求得
地球物理数据处理教程
将上面四项代入（9.5.7）式，并考虑到
地球物理数据处理教程
而（9.5.7）式被积函数的第二项为
地球物理数据处理教程
利用关系式
地球物理数据处理教程
其中
地球物理数据处理教程
可得
地球物理数据处理教程
为计算（9.5.7）式面积分中第三项，由表8-1知，f只有两种形式。当区域内无场源时f=0，否则f有Iδ（x）δ（z）的形式。这里我们考虑点源二维电阻率法的情况，f=
地球物理数据处理教程
考虑两种情况：一是i点为供电点，即x0=xi，z0=zi，并由狄拉克函数的定义和N函数的特性［8］Ni（xi，zi）=1，将这些关系代入可得
地球物理数据处理教程
这里Δi为i点周围小面积之和；二是i点不是供电点，则该项就为零。
图9.7 边界线段jm上u的变化示意图
对于边界单元，除了面积分部分以外，（9.5.7）式还有线积分部分：
地球物理数据处理教程
这里i可以代表任意节点，但在边界单元中往往以j、m节点代表边界节点。为了避免混淆，此处将上式中i写为l，即
地球物理数据处理教程
jm为边界线段。为了计算上式，同样设u在jm上为线性变化，如图9.7，则有
u=（1-t）uj+tum（0≤ t≤1）
记jm边长为
地球物理数据处理教程
则
地球物理数据处理教程
这样，对于j点线积分项可写为
地球物理数据处理教程
同样对于m点，有
地球物理数据处理教程
总结以上各式，对于面积分，讨论了l=i时
地球物理数据处理教程
式中单元系数矩阵的各元素为
地球物理数据处理教程
上式中只对边界单元才计算系数中最后的αγsi项。这里指的边界单元不包括地面边界，因为对地面边界，γ=0，无须计算边界项。
将D域中各单元的单元系数矩阵元素计算出来，再进行对应元素的叠加，可得到总体系数矩阵，利用（9.5.6）式建立方程，并将与电流I有关的项移到右端，得到
地球物理数据处理教程
式中I对应供电点所在节点。
计算出系数矩阵后，求解方程组（9.5.9），即可得到目标函数ul值（l=1，2，…l0）。
9.5.2 网格
解与赫姆霍兹方程相对应的有限元法，其步骤同一般前述的有限元法，但由于这个方法对电法实际问题有重要意义，所以要较详细地说明其中一些问题。
由前述可知，同样用三角形单元对研究区域进行剖分，节点与单元遍布整个D域，形成网格。这样将连续场函数的求解化为节点上离散函数值的求解。因此，网格的大小、网格的类型和网格的疏密，对函数值的计算会产生重要的影响。下面分别讨论。
9.5.2.1 网格的大小和疏密
网格大小就是所取泛函积分域D的大小，一般说来网格越大越好。对于微分方程边值问题的求解，只有给出正确的边界条件，才能求出域中比较精确的函数值。例如，对于使用第一类边界条件时，要求给定正确的边界函数值。但是，对于不均匀的地电断面，边界函数值和其他边界条件值，特别是地下部分，均无法正确求出。因此，往往采用地下边界值给零，或由均匀地电条件给值。这时就要求区域D的边界远离不均匀区，即要求网格大，否则就会影响计算的精度。但是，另一方面，如果网格内单元的大小不变的话，那么网格太大，势必要大量的增加节点数，从而需要更多的计算机内存和增加计算工作量，这是因为节点个数的增加，将使线性方程组的阶数增加，于是系数矩阵的形成，解方程的时间都要增加。程序中占机器内存最多的是系数矩阵的元素个数，而该数目是节点数的平方。由此可见，网格大小的选取兼顾计算精度和计算工作量几个方面由试验确定。
对于网格内部单元的大小，一般说来单元越小，计算精度越高。这是因为我们假定u函数在每个单元内呈线性变化。如果单元太大，实际函数便可能不满足这个条件，从而增加计算误差。另外，我们还假定单元内电性是均匀的，即电性参量σ为常数，这也要求单元较小，特别是要拟合复杂的地电断面和地形剖面，更需要划分得细致些，才能满足单元内均匀的条件。用点源二维有限元法计算了二层水平介质的视电阻率，选用了两种不同的单元大小，其计算精度如表9.1所列。表中相对误差指计算值与理论值的差与理论值之比，表中的计算结果也说明单元越小，计算效果越好。但同样，如果网格大小不变，单元变小则整个节点个数就增加，从而影响内存和计算量的增加。
表9.1 不同单元大小计算精度
为了克服网格大小和单元大小选择中精度和工作量之间的矛盾，我们采用非均匀的网格，即网格的中心部分单元小、节点密，边部单元大、节点稀，由中心到边缘单元逐渐放大，这样既保证了网格有足够的大小，又保证地电断面的复杂部位位于网格中心密区，以满足单元内电性均匀和u函数线性变化的条件。单元逐渐放大，是为了避免将三角形单元拉成窄长的形状。
9.5.2.2 网格的类型
通过研究发现网格类型对计算结果的精度有一定影响，现以二层介质地电断面的电阻率法计算为例，当用偶极测深装置时，不同网格的计算结果如表9.2所示，现讨论如下。
表9.2 网格类型Ⅲ水平二层计算的偶极剖面电阻率值与理论值
注：BM表示偶极剖面AB—MN装置B、M极之间的距离。
第一方案，按同一个方向划分三角形单元，网格如图9.8所示，计算结果发现，当供电点位于网格中心（如图9.8中第4号节点）时，左边（如图9.8 中第3 号节点）的电位值大于右边（如图9.8中第5号节点）的电位值，其具体数值列于表9.3中第一行，出现了电位值明显的不对称现象，用此电位值计算的偶极测深视电阻率曲线如图9.9中的黑圆点所示，与理论曲线相差较大（理论曲线为实线）。
图9.8 网格类型Ⅰ
表9.3 网格类型Ⅳ对应的水平二层计算的偶极剖面电阻率值与理论值
注：表上“左”“右”指与供电点相邻的一个节点。
图9.9 各种网格类型计算效果对比图
1—网格类型Ⅰ；2—网格类型Ⅱ；3—网格类型Ⅲ；4—网格类型Ⅳ
第二方案，以网格中心点为对称点，两边分别向不同方向划分，网格如图9.10所示。计算结果发现，当供电点位于网格中心时，左边和右边的电位值相等，完全对称；当供电点在网格中心的左边时，供电点以左的电位大于供电点以右的电位；当供电点位于网格中心的右边时，供电点以左的电位小于供电点以右的电位，而且供电点在网格两边对称位置上时，所出现的不对称现象刚好相反，其数值完全相同。如表9.3 中第二、三、四行所示，这与网格单元形状的对称性正好吻合，用这组电位值计算偶极测深视电阻率曲线，如图9.9中的“+”符号所示，与理论曲线相差仍较大。
第三方案，以网格中心为对称，两边用两个方向交替的波浪形划分三角单元，网格如图9.11所示。计算结果发现，当供电点在网格中心时，两边电位几乎相等，基本对称；当供电点位于网格的左边或右边时，不对称的现象较Ⅰ和Ⅱ有明显的减少，但仍不完全对称，并随着供电点向两侧移动而不对称现象加大，这种网格单元形状虽然各点均对称，但是各供电点两边所包涵的网格面积并不对称，试验所采用的小网格，更是能引起这种较小的不对称现象的原因，其具体例子见表9.3中第五、六、七、八、九、十行，用这组电位值计算偶极视电阻率曲线，如图9.9中的“▲”符号所示，且发现在不同的供电点上所得到的视电阻率是起伏变化的。以不同极距的偶极剖面的视电阻率值为例说明，水平二层的偶极剖面视电阻率值理论上应为某一常数，而从表9.2中可见，却为一组起伏的大小波浪相间的数值，与网格单元划分波浪形状有相似的规律。
图9.10 网格类型Ⅱ
图9.11 网格类型Ⅲ
图9.12 网格类型Ⅳ
第四方案，按照以网格中心为对称、各点均对称的两个方向交叉划分三角单元，简称交叉对称网格，如图9.12所示。计算结果表明，供电点在网格中心时严格对称，供电点在网格两边时也接近对称，其数值见表9.3中第十一、十二、十三行，用这组电位计算的偶极测深视电阻率曲线如图9.9中虚线所示，与理论曲线拟合最好。
对比上述四种单元形状，不难看出，有限单元的计算结果与单元的划分、形状的对称性均有一定的关系。而对比的结果发现，交叉对称网格是比较合理的网格，这种划分单元也较细致，对于模型复杂的地电断面是方便的，对模拟各种地形也是较容易的。不过遗憾的是，由于交叉形状使内部节点的个数相对的增加了接近一倍，势必会较多地增加计算工作量和内存。为此，对系数矩阵作了必要的处理，从而使计算工作量和内存基本上不增加。
在实际中所采用的网格如图9.13所示。
图9.13 点源二维电阻率法计算网格
从图9.13中可见，这种网格是综合分析了以上各方面试算的结果选择出来的比较合理的形式，网格的总面积、密区、稀区、总节点数、节点间隔距离均是可变的。
9.5.3 边界条件
为了确定具体的电位和电场分布，必须要提出定解边界条件。前面已提出有三类边界条件，在边界上给定函数值称为第一类边界条件，给定函数值的边界法向导数称为第二类边界条件，综合给定函数值和其法向导数的称为第三类边界条件。这些条件的给定依赖于具体的物理问题和函数特点等。边界条件正确与否对计算有重要的影响。这里仅以点源二维电阻率法为例，来说明边界条件的给定和处理方法。对于点源二维电阻率法，本节（8.2.16）式中目标函数为变换电位φ。由于第二类边界条件为自然边界条件，有限元解微分方程时自然满足，所以这里只讨论第一和第三类边界条件。
9.5.3.1 强加边界条件
第一类边界条件又称为强加边界条件，我们通常采用两种方法给定，其一是将除地面以外的地下边界上的φ（x，k，z）取为零。这样在程序处理中最方便，但要求将网格取得足够大，否则便会由此带来计算误差。因为理论上只有在远离源无穷远处φ函数才等于零。其二是将地下边界上φ值取为均匀介质的φ值。均匀介质任一点P（x，z）的φ值由（8.2.21）式计算，即
φ（x，k，z）|Γ=AK0（k·r）
式中r=
当一部分边界节点电位已知时，这时要对方程组（9.5.9）进行加工，设（9.5.9）式写为以下形式
Aφ=G
设网格的总节点为l0，内部待求电位的节点个数为l1，在边界上已知电位的节点个数为l0-l1。为书写简单计，设已知电位的节点编号在最后（实际计算往往不必如此），则上述方程组可写为
地球物理数据处理教程
式中
地球物理数据处理教程
在实际计算中将已知电位的节点在A阵中划行划列，即将该节点对应的对角线元素设为1，将对应的行和列的其他元素全设为零。
9.5.3.2 混合边界条件
第三类边界条件又称为混合边界条件。柯岗的试验证明，在半无限空间中半水平圆柱体的异常电位（在均匀一次场中），用水平圆柱体截面直径的10倍宽度的网格进行计算，用自然边界条件计算时结果高于真值（理论值），而用地下边界给零电位条件计算时结果低于真值，如图9.14。而这两种计算结果的平均值与真值非常靠近，说明应该联合使用第一和第二类边界条件，这就形成了第三类边界条件。
图9.14 半无限空间水平圆柱体的异常电位（均匀一次场）
1—理论曲线；2—自然边界条件计算结果；3—地下边界给零的计算结果；4—2、3两种结果的平均
图中0.3为半圆截面的边界，3.0为网格边界
对均匀半空间点源二维电阻率问题，混合边界条件由（8.2.23）式给定，即
地球物理数据处理教程
但是对于非均匀的地电断面，便导不出上式条件。为此，采用更一般形式的混合边界条件，即
地球物理数据处理教程
式中λ为与地电断面不均匀性等有关的一个系数。这时（9.3.6）式成为
地球物理数据处理教程
对均匀介质，考虑到电源点附近误差在计算中的传播和反傅氏变换的实际等问题，还是不应该取λ=1。当然，无论是均匀介质还是不均匀介质，λ值均难于从理论上确定。因此，将λ视为特定系数，由取得最佳计算结果的试验求出。大量的研究发现λ取值在0～1范围内，均值约为1/2。
图9.15示出了垂直接触带中梯ρs曲线的计算结果。由图可见λ=1/2时的计算结果与理论ρs值误差在1%以内；而λ=1时的计算误差可达4%，在界面附近ρs曲线形态有畸变。
为了对比强加边界条件和混合边界条件对计算精度的影响，选用了节点个数为69×13的网格对均匀介质试验了这两种边界条件，采用地下边界赋零的条件时计算结果与理论值相对误差为2% ～5%，而取混合边界条件时误差为1%左右。
图9.15 垂直接触带中梯ρs曲线
AB=38，MN=1，ρ1=100Ω·m，ρ2=50Ω·m
1—理论曲线；2—λ=0.5时的计算结果；3—λ=1时的计算结果
在式（9.5.12）边界条件中，γ由（9.2.23）式给出，其中包括了第二类零阶和一阶修正贝塞尔函数。
9.5.4 系数矩阵的简化
图9.16 网格中的任一个小长方形示意图
选用交叉对称网格来离散泛函积分域，由于交叉点使计算工作量和内存均要增加，通过下述的公式变形，将线性方程组未知数中包括交叉点的电位φ的那些项消去，使未知数的个数降为不包括交叉点的节点个数，这就使计算工作量和内存基本上不增加。
设网格中任一个小长方形如图9.16所示，各节点在总的网格中的序号不妨设为1、2、11、12，交叉节点设为α，小长方形的长为DX，宽为DZ，长方形分成4个小三角形，三角形2、α、12为Δ1；2、α、1为Δ2；1、α、11为Δ3；11、α、12为Δ4。它们的面积均相同，即为
Δ1=Δ2=Δ3=Δ4=DX·DZ/4=Δ
它们的导电率分别为σ1、σ2、σ3、σ4。
按照前面公式（9.5.8）的形式，对于点源二维电阻率法，采用混合边界条件λ
在Δ3中
地球物理数据处理教程
同样在Δ2中
地球物理数据处理教程
从图9.16中可见，系数矩阵中k1，1只有在Δ2和Δ3两个三角形中才对它有贡献，因此k1，1应为二者之和
地球物理数据处理教程
类似的可写出其他各对角线元素如下
地球物理数据处理教程
对kαα，应由4个三角形对它作贡献：
在Δ1中
地球物理数据处理教程
在Δ2中
地球物理数据处理教程
在Δ3中
地球物理数据处理教程
在Δ4中
地球物理数据处理教程
求和即得
地球物理数据处理教程
对于非对角线元素分别如下
在Δ2中
地球物理数据处理教程
类似的可写出
地球物理数据处理教程
下面列出与α点有关的系数
在Δ3中
地球物理数据处理教程
在Δ2中
地球物理数据处理教程
在Δ1中
地球物理数据处理教程
在Δ4中
地球物理数据处理教程
将所有贡献叠加起来，即得各系数的表达式如下：
地球物理数据处理教程
下面将线性方程组中，以这5个节点的φ为未知数的部分列出来。
设对应于节点1、2、11、12、α的φ分别为φ1、φ2、φ11、φ12、φα，线性方程组中有关部分为：
地球物理数据处理教程
通过消元一次，将φα在（9.5.14）式中下面4个方程中消去，于是得到新的一组方程如下
地球物理数据处理教程
去掉第一个方程式，即得到仅包括φ1、φ2、φ11、φ12四个未知数的4个方程式，将交叉节点α的φα从未知数中去掉了，于是系数矩阵和最后解方程均用不着考虑交叉节点α，仅仅是在形成系数矩阵的各元素时，按新的系数表达式（9.5.15）计算即可。
式（9.5.15）新的一组系数表达式应为：
地球物理数据处理教程
地球物理数据处理教程
相应的混合边界条件的修改项如下（取λ=
对应于左边界：
k′1，1=k′1，1+α·DZ·σ2/3
k′2，2=k′2，2+α·DZ·σ2/3 （9.5.17）
k′1，2=k′1，2+α·DZ·σ2/6
对应于右边界：
k′11，11=k′11，11+α·DZ·σ4/3
k′12，12=k′12，12+α·DZ·σ4/3 （9.5.18）
k′11，12=k′11，12+α·DZ·σ4/6
对应于底边界：
k′1，1=k′1，1+α·DX·σ3/3
k′11，11=k′11，11+α·DX·σ3/3 （9.5.19）
k′1，11=k′1，11+α·DX·σ3/6
通过上述处理以后，采用交叉对称网格，计算工作量和内存基本不增加，所以计算时间是较少的，精度亦较高。
9.5.5 系数矩阵的特点、 存放方式及解方程方法的对比
由前所述，可知系数矩阵是一个阶数等于节点总数的方阵，其元素个数为节点数的平方，而且矩阵中大量的元素为0，是一个大型对称正定矩阵，且非零元素分布在对角线附近的一个条带内，如图9.17所示，该矩阵对应于图9.13的网格划分形式。网格节点个数为190个，该矩阵的行为190行，列为190列，其中非零元素集中在对角线附近的24个元素以内，这类矩阵数学上称带状矩阵，将对角线元素到该行最远距离的一个非零元素之间的总元素个数称为半带宽（包括对角线元素和最远的那个元素在内），以S表示，如图9.17中S应为12。
网格节点编号如图9.13所示，是先沿z方向排列，然后再沿x方向排列，于是带宽S应为z方向节点个数加2。根据系数矩阵的这些特点，已试验了以下三种存放方式以及解方程的方法；
（1）采用全部元素存放，用二维数组形式，用一般的高斯消元法求解方程，求解的精度较高，误差小于10-9，但要求的内存太大，计算花时间较多。
（2）只存放下三角中带宽以内的元素，用一维数组按行顺序存放，采用高斯消元法求解，这种存放虽然省了内存，却增加了寻找各元素在一维数组中位置的时间，因而也不实用。
图9.17 系数矩阵元素分布示意图
0表示零元素；Δ表示非零元素
图9.18 系数矩阵带宽存放示意图
图中S表示带宽
（3）只存放上三角或下三角带宽内元素，用二维数组按图9.18形式存放，二维数组的列为半带宽的个数，行为节点个数，图9.18中虚线部分的元素充零。
这种存放既节省内存单元，取数也比较方便，相应解方程的方法，通过高斯消元法、LU、LDLT分解法、LLT分解法试算结果发现，它们的精度均较高，误差小于 10-9；速度以LLT分解法为最快。
通过试算研究确定，在实用程序中选用第三种存放方式，即二维带宽存放，如图9.18所示，解方程的方法采用LLT分解法，而且对于不同的供电点，采用分解一次系数矩阵（即将供电点近似地视为位于网格中心），不同的供电点仅仅改变右端项，分别进行回代即可。
本回答被网友采纳