与赫姆霍兹方程对应的二维有限元法

解拉普拉斯方程的二维有限元法~

9.4.1 基本原理
对于直流线源二维地电问题，电位满足的微分方程由（8.2.14）式表示

地球物理数据处理教程

上式是在取y轴平行于线源和地质体走向时得到的。为简单计，我们将电源挖掉，即计算区域内不包含电流源项，在地面边界取绝缘边界条件，并在其它边界取强加边界条件，例如在挖去电流边界处取理论计算值，在外边界取零，这时（8.2.14）式变为

地球物理数据处理教程

将计算区域进行有限单元剖分，在每个单元中设 σ 为常数，这样在每个单元内（9.4.1）式变成拉普拉斯方程

地球物理数据处理教程

这样，对于线源二维问题，求解u（x，z）的问题也就是求解满足上面边界条件的拉普拉斯方程的问题，这是电法正演中的最简单情况。
已经证明，这类问题可以等价于求下列泛函的极小值问题

地球物理数据处理教程

用有限单元法求解使泛函取极小值的u，其过程与一维有限元法类似。
9.4.2 区域剖分

图9.5 三角形单元

常用的办法是将D域剖分成许多小的三角形单元（当然也可以剖分成其他形状的单元，但三角形单元最简单），某一个三角形单元如图9.5所示。
当边界是曲线时，用三角元的一边近似，在遇到内部介质分界线时，不容许三角元跨越分界线；此外，不容许三角元顶点落在其它三角元的边上；还要避免出现太尖、太钝的三角元。在u变化大的地方，三角元密一些，反之稀一些。
各单元的顶点称为节点，将所有节点和三角元分别按一定顺序编号。为下面分析需要，对任一三角元e的三顶点（节点），按逆时针编号为i，j，m，其坐标为（xi，zi），（xj，zj），（xm，zm），其函数值为ui，uj，um。除第一类边界条件给定的边界节点上的函数值是已知的外，其余节点上的函数值都是待求的。这样一来，将连续函数u（x，y）的求解化成节点上离散函数值的求解。
9.4.3 线性插值
由于各三角形单元取得足够小，使得可以假定各单元内电位随坐标是线性变化的，即u（x，z）是线性函数，有
u=α1+α2x+α3z （9.4.3）
同时对三个顶点有

地球物理数据处理教程

按克莱姆法则解上述方程组可求得

地球物理数据处理教程

其中

地球物理数据处理教程

αi=xjzm-xmzj
αj=xmzi-xizm
αm=xizj-xjzi
bi=zj-zm
bj=zm-zi
bm=zi-zj
ci=xm-xj
cj=xi-xm
cm=xj-xi

地球物理数据处理教程

将α1、α2和α3代入（9.4.3）式中，整理可得

地球物理数据处理教程

或写为
u（x，z）=Ni（x，z）ui+Nj（x，z）uj+Nm（x，z）um（9.4.4）
式中：

地球物理数据处理教程

以上式中Δe为三角形单元的面积。
9.4.4 单元分析及总体合成
只研究（9.4.2）式对单元 e 的积分，记作 Je。将（9.4.4）式代入（9.4.2）式中并沿整个单元积分，为求得使（9.4.2）式取极小的函数 u，对任意节点微分（9.4.2）式可得

地球物理数据处理教程

由（9.4.5）式有关系

地球物理数据处理教程

而

地球物理数据处理教程

再由（9.4.5）式有关系

地球物理数据处理教程

将以上关系代入上式积分中，由于积分号内均为常数，可以提出积分号外，而

地球物理数据处理教程

最后可得

地球物理数据处理教程

同理可得

地球物理数据处理教程

或改写为

地球物理数据处理教程


地球物理数据处理教程

矩阵［k］e称为单元系数矩阵，其中元素

地球物理数据处理教程

由于krs=ksr，所以k阵为对称矩阵，［u］e是单元节点电位值组成的列向量。
如果节点r不属于单元e，则Je（u）中不含有ur，所以 =0，于是我们将（9.4.7）式扩展成

地球物理数据处理教程

（矩阵中的虚点均为零元素）
扩展后的矩阵是l0阶方阵，l0是节点总数；扩展后的列向量 ｛u｝ 由全体节点函数值组成。
将所有单元Je（u）合成，得到整个区域的J（u）

地球物理数据处理教程

它是所有节点函数值u1……ul0的函数，写成
J（u）=J（u1……ul0）
也可将J（u）看成变量u1……ul0的多元函数，所以泛函取极值相当于多元函数取极值

地球物理数据处理教程


图9.6 两个三角形单元的示意图

①单元节点编号：1、2、3；②单元节点编号：3、2、5
即分别对各单元求偏导数，然后合成。对每一个单元求偏导数都可获得形如（9.4.8）的式子。两个三角元的顶点不同，式（9.4.8）的矩阵中的非零元素所占的行、列也不同，但等式两边的列向量完全相同，所以将全部三角元的偏导数合成时，只要将矩阵中的对应元素加起来就可以了。例如：第①单元的节点编号为1、2、3，第②单元节点编号为3、2、5（逆时针顺序），如图9.6，这两个三角元的偏导数按如下方法合成：

地球物理数据处理教程


地球物理数据处理教程

将全部单元合成，有

地球物理数据处理教程

简记成［k］［u］=0。n阶矩阵［k］称为总体系数矩阵，它由许多对称矩阵合成，所以也是对称的。此外，数学上还可证明是正定的。从式（9.4.8）可见，非零元素只存在于三角元三顶点编号所对应行和列的九个交叉位置上，其他均为零元素，所以总体系数矩阵［k］中包含着大量零元素，称为稀疏矩阵。离主对角线最远的非零元素的位置取决于所有三角元顶点编号的最大差值R（绝对值）。
最后，解线性代数方程组，矩阵方程实际上就是m个线性代数方程，其中m是待求函数值的节点数，解这个方程组，即可求得这些节点上的函数值。

与赫姆霍兹方程对应的二维有限元法在电法勘探中有较广的使用范围，有重要的意义。对电阻率法，用点源二维有限元方法对不同的情况进行了试算和应用，取得了较好的效果。
9.6.1 理论对比
图9.19中示出了二层介质时偶极测深装置有限元法计算的视电阻率ρs曲线与理论曲线的对比，图中实线为理论曲线，黑点为计算结果，地电断面和装置均附在图中。由图可见，计算值与理论值符合很好，计算误差在1%以内。

图9.19 二层ρs偶极测深曲线

图9.20示出了对两种不同电阻率介质的垂直接触带上偶极测深视电阻率ρs曲线的计算结果，与理论曲线对比，计算误差在2%以内。图中实线为理论曲线，黑点为计算结果。

图9.20 垂直接触带ρs偶极测深曲线

9.6.2 模型试算结果
为检验前述算法，对大地水平，即在没有地形影响的情况下，设置了以下几种模型（图9.21、图9.22、图9.23、图9.24），每个模型的参数标注在模型下，采用对称四极测深和温纳装置进行了试算。其中对以上设计的前三种模型都采用对称四极斯伦贝尔装置，其最大电极距为25m，最小电极距为1m。后一种模型采用温纳装置，最大电极距为24m，最小电极距为1.5m。试算的结果如图所示。
模型1：设计了三层，第一层和第三层的电阻率都是100Ω· m，第二层的电阻率是10Ω·m，第二层的中心埋深h=4m。

图9.21 模型1与视电阻率等值线图


图9.22 模型2与视电阻率等值线图

模型2：在大地水平面下有一个形状为正方体的物体，其边长为6m，中心埋深h=5m，电阻率为10Ω·m，围岩电阻率为100Ω·m。
模型3：在大地水平面下有两个形状为正方体的物体，其边长为5m，中心埋深h=4.5m，电阻率为10Ω·m，围岩电阻率为100Ω·m。
模型4：在大地水平面下有三个截面形状为正方形的物体，其长为6m，高为6m，中心埋深为5m，其电阻率为10Ω·m，围岩电阻率为100Ω·m，用有限单元法计算的温纳装置下的视电阻率断面等值线图如图9.24所示。

图9.23 模型3与视电阻率等值线图


图9.24 模型4与温纳装置视电阻率断面等值线图

从计算的结果看是较好的，计算精度较高，视电阻率等值线图较好地反映了地下物体的电性分布。而且，这使结果的分析和解释变得直观和形象。
9.6.3 模拟实际电测深曲线的一个结果
四川省某地云母矿，产出于高阻伟晶岩脉内。其中，一个已知矿区的实测三极电剖面ρs曲线示于图9.25中，图下方为已知地质剖面，上方为实测ρs曲线，虚线为点源二维有限元的计算结果，图中每个点号之间的距离为10m，三极剖面的极距AO=80m。
图9.26中示出了有限元法计算时所取的地电断面，其上方亦为计算ρs曲线，该地电断面是根据实际地质情况、物性数据并考虑电测深曲线特点而选取的，其中取了5条高阻岩脉，还有一些高阻地表滚石层。

图9.25 某地云母矿上三极剖面ρs曲线

1—实测值；2—有限元计算结果

图9.26 模拟时所采用的地电断面和有限元计算结果

ρ1=1000Ω·m，ρ2=26000Ω·m，ρ3=12000Ω·m
左面4个高阻岩脉下延深度160m

前面讲过，许多重要的电法问题，如点源二维电阻率法和激发极化法，二维大地电磁问题，二维的甚低频电磁场及线源问题等有关的电位或电磁场均可以用（8.3.26）式或（9.3.1）式二维赫姆霍兹方程来描述，因此与此相对应的有限元法在电法勘探数字模拟的实际应用上，具有重要的意义。

我们知道，由（9.3.1）式二维赫姆霍兹方程

地球物理数据处理教程

在第三类边界条件

地球物理数据处理教程

的解与满足泛函

地球物理数据处理教程

极小的函数相对应。

式中：u为目标函数；D为研究区域；Γ为D的边界。

9.5.1 系数矩阵的导出

当用有限元法求解泛函式（9.5.3）极小时，其步骤和上节中与拉普拉斯方程对应的二维有限元法相同。在D域的离散中仍然采用最简单的三角形剖分，划分的方法和原则均同上节所述，这里，我们研究某个三角形单元e，其顶点分别按逆时针编号，记为i，j，m，相应的坐标为（x_i，z_i）、（x_j，z_j）、（x_m，z_m），其函数值为u_i、u_j、u_m，设u=u（x，z）在单元内线性变化，同样可得到（9.4.4）式由三角形顶点函数值表示的单元e中的函数值为

u（x，z）=N_i（x，z）u_i+N_j（x，z）u_j+N_m（x，z）u_m （9.5.4）

式中

地球物理数据处理教程

地球物理数据处理教程

称为基函数，其中

α_i=x_jz_m-x_mz_j，b_i=z_j-z_m，c_i=x_j-z_j

α_j=z_mz_i-x_iz_m，b_j=z_m-z_i，c_j=x_i-z_m

α_m=x_iz_j-x_jz_i，b_m=z_i-z_j，c_m=x_j-z_i

Δe=

（b_ic_j-b_jc_i）为单元e的面积。

为了计算三角形中的泛函值，必须先计算一阶导数，利用（9.5.4）式可写出

地球物理数据处理教程

其中

地球物理数据处理教程

将（9.5.3）式的积分域D分成各小单元三角形之和，即

地球物理数据处理教程

这里J_e为一个小三角形单元e内的泛函

地球物理数据处理教程

设α、β、γ均为常数，可以提到积分号外，如果单元e在D域内部，则第二项的线积分为零，若三角形单元e有一个边在边界上，则要计算线积分部分。为了研究有一般性意义，不妨设单元e的一个边在D域的边界。

其目的在于求目标函数 u，使泛函 J（u）取极小，按照极小必要条件，应有下式成立：

地球物理数据处理教程

式中l₀为D域中所有节点总和。为了求出上式的具体形式，先考虑（9.5.5）式对u_l的一阶导数，为计算方便，不妨先设l=i。

地球物理数据处理教程

由前面引出的

，

的表达式，可求得

地球物理数据处理教程

同理，

地球物理数据处理教程

利用（9.5.4）式，可求得

地球物理数据处理教程

将上面四项代入（9.5.7）式，并考虑到

dxdz=Δe，在（9.5.7）式中，对中括号中的被积函数积分后可写为

地球物理数据处理教程

而（9.5.7）式被积函数的第二项为

地球物理数据处理教程

利用关系式

地球物理数据处理教程

其中

地球物理数据处理教程

可得

地球物理数据处理教程

为计算（9.5.7）式面积分中第三项，由表8-1知，f只有两种形式。当区域内无场源时f=0，否则f有Iδ（x）δ（z）的形式。这里我们考虑点源二维电阻率法的情况，f=

δ（x）δ（z），所以第三项可写为

地球物理数据处理教程

考虑两种情况：一是i点为供电点，即x₀=x_i，z₀=z_i，并由狄拉克函数的定义和N函数的特性^［8］N_i（x_i，z_i）=1，将这些关系代入可得

地球物理数据处理教程

这里Δ_i为i点周围小面积之和；二是i点不是供电点，则该项就为零。

图9.7 边界线段jm上u的变化示意图

对于边界单元，除了面积分部分以外，（9.5.7）式还有线积分部分：

地球物理数据处理教程

这里i可以代表任意节点，但在边界单元中往往以j、m节点代表边界节点。为了避免混淆，此处将上式中i写为l，即

地球物理数据处理教程

jm为边界线段。为了计算上式，同样设u在jm上为线性变化，如图9.7，则有

u=（1-t）u_j+tu_m（0≤ t≤1）

记jm边长为

地球物理数据处理教程

则

地球物理数据处理教程

这样，对于j点线积分项可写为

地球物理数据处理教程

同样对于m点，有

地球物理数据处理教程

总结以上各式，对于面积分，讨论了l=i时

的各项计算公式，对l=j，l=m时类似的方法可得相应的计算公式。加上对边界单元的讨论，可以归纳如下，对于任何一个小三角形元，

（l=i，j，m）的计算，有以下矩阵形式

地球物理数据处理教程

式中单元系数矩阵的各元素为

地球物理数据处理教程

上式中只对边界单元才计算系数中最后的αγs_i项。这里指的边界单元不包括地面边界，因为对地面边界，γ=0，无须计算边界项。

将D域中各单元的单元系数矩阵元素计算出来，再进行对应元素的叠加，可得到总体系数矩阵，利用（9.5.6）式建立方程，并将与电流I有关的项移到右端，得到

地球物理数据处理教程

式中I对应供电点所在节点。

计算出系数矩阵后，求解方程组（9.5.9），即可得到目标函数u_l值（l=1，2，…l₀）。

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）|_Γ=AK₀（k·r）

式中r=

，x、z为供电点到边界P点的坐标差。这种取法也要求网格较大，因为忽略了实际介质的非均匀性对网格边界节点的影响。当用强加边界条件时，在解方程组求得电位φ以前，要对方程组进行必要的处理。

当一部分边界节点电位已知时，这时要对方程组（9.5.9）进行加工，设（9.5.9）式写为以下形式

Aφ=G

设网格的总节点为l₀，内部待求电位的节点个数为l₁，在边界上已知电位的节点个数为l₀-l₁。为书写简单计，设已知电位的节点编号在最后（实际计算往往不必如此），则上述方程组可写为

地球物理数据处理教程

式中

为节点i上的已知电位值。实际未知电位为φ₁…φ_l1的节点个数共l₁个。这样可将方程组的常数项右移叠加在右端项中，即可得下面方程组

地球物理数据处理教程

在实际计算中将已知电位的节点在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，ρ₁=100Ω·m，ρ₂=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为Δ₁；2、α、1为Δ₂；1、α、11为Δ₃；11、α、12为Δ₄。它们的面积均相同，即为

Δ₁=Δ₂=Δ₃=Δ₄=DX·DZ/4=Δ

它们的导电率分别为σ₁、σ₂、σ₃、σ₄。

按照前面公式（9.5.8）的形式，对于点源二维电阻率法，采用混合边界条件λ

+γφ=0，这时求取各三角形单元对系统矩阵的贡献如下：

在Δ₃中

地球物理数据处理教程

同样在Δ₂中

地球物理数据处理教程

从图9.16中可见，系数矩阵中k_1，1只有在Δ₂和Δ₃两个三角形中才对它有贡献，因此k_1，1应为二者之和

地球物理数据处理教程

类似的可写出其他各对角线元素如下

地球物理数据处理教程

对k_αα，应由4个三角形对它作贡献：

在Δ₁中

地球物理数据处理教程

在Δ₂中

地球物理数据处理教程

在Δ₃中

地球物理数据处理教程

在Δ₄中

地球物理数据处理教程

求和即得

地球物理数据处理教程

对于非对角线元素分别如下

在Δ₂中

地球物理数据处理教程

类似的可写出

地球物理数据处理教程

下面列出与α点有关的系数

在Δ₃中

地球物理数据处理教程

在Δ₂中

地球物理数据处理教程

在Δ₁中

地球物理数据处理教程

在Δ₄中

地球物理数据处理教程

将所有贡献叠加起来，即得各系数的表达式如下：

地球物理数据处理教程

下面将线性方程组中，以这5个节点的φ为未知数的部分列出来。

设对应于节点1、2、11、12、α的φ分别为φ₁、φ₂、φ₁₁、φ₁₂、φ_α，线性方程组中有关部分为：

地球物理数据处理教程

通过消元一次，将φ_α在（9.5.14）式中下面4个方程中消去，于是得到新的一组方程如下

地球物理数据处理教程

去掉第一个方程式，即得到仅包括φ₁、φ₂、φ₁₁、φ₁₂四个未知数的4个方程式，将交叉节点α的φ_α从未知数中去掉了，于是系数矩阵和最后解方程均用不着考虑交叉节点α，仅仅是在形成系数矩阵的各元素时，按新的系数表达式（9.5.15）计算即可。

式（9.5.15）新的一组系数表达式应为：

地球物理数据处理教程

地球物理数据处理教程

相应的混合边界条件的修改项如下（取λ=

）。

对应于左边界：

k′_1，1=k′_1，1+α·DZ·σ₂/3

k′_2，2=k′_2，2+α·DZ·σ₂/3 （9.5.17）

k′_1，2=k′_1，2+α·DZ·σ₂/6

对应于右边界：

k′_11，11=k′_11，11+α·DZ·σ₄/3

k′_12，12=k′_12，12+α·DZ·σ₄/3 （9.5.18）

k′_11，12=k′_11，12+α·DZ·σ₄/6

对应于底边界：

k′_1，1=k′_1，1+α·DX·σ₃/3

k′_11，11=k′_11，11+α·DX·σ₃/3 （9.5.19）

k′_1，11=k′_1，11+α·DX·σ₃/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、LDL^T分解法、LL^T分解法试算结果发现，它们的精度均较高，误差小于 10^-9；速度以LL^T分解法为最快。

通过试算研究确定，在实用程序中选用第三种存放方式，即二维带宽存放，如图9.18所示，解方程的方法采用LL^T分解法，而且对于不同的供电点，采用分解一次系数矩阵（即将供电点近似地视为位于网格中心），不同的供电点仅仅改变右端项，分别进行回代即可。

---

冀州市18745612413： matlab编程,有限元方法求解Helmholtz方程的边值问题用,矩阵法 - ？
温艺双姜： clc clear tic; % n 是行与列划分的格子数,对整个[0,1]*[0,1]有n^2个划分,can为方程中的参数k %这里我们用1,2,3按逆时针来表示一个三角形的各个顶点 % s是一个n^2*10的关联矩阵,s(i,1)表示第i个三角形,s(i,2),s(i,3),s(i,4)分别表示第i个三角形...

冀州市18745612413： 剪力墙受力特点及构造要求 - ？
温艺双姜： 在水平荷载作用下,这类剪力墙截面上的正应力分布略偏离了直线分布的规律,变成了相当于在整体墙弯曲时的直线分布应力之上叠加了墙肢局部弯曲应力,当墙肢中的局部弯矩不超过墙体整体弯矩的...

冀州市18745612413： 极限平衡法包括哪些?极限平衡法与有限元法的差异 - ？
温艺双姜： 目前常用的二维极限平衡分析方法有:瑞典法(亦称作Fellenious法)、简化Janbu法、Bishop简化法、严格Janbu法、Lowe-Karafiath(罗厄)法、美国陆军工程师团法、Morgenstern-Price法、Spencer法、垂直条分Sarma法、斜条分Sarma法...

冀州市18745612413： Galerkin方法与有限元法的关系 - ？
温艺双姜： 提出一维弹性动力问题的局部Petrov-Galerkin方法,这是一种真正的无网格方法.这种方法采用移动最小二乘函数采近似解变量,并且采用移动最小二乘近似函数的权函数作为加权残值法的权函数.文中对形成的离散动力学方程用Newmark方法求...

冀州市18745612413： 元分析的步骤 - ？
温艺双姜： 有限元分析(FEA,Finite Element Analysis)的基本概念是用较简单的问题代替复杂问题后再求解.它将求解域看成是由许多称为有限元的小的互连子域组成,对每一单元假定一个合适的(较简单的)近似解,然后推导求解这个域总的满足条件...

冀州市18745612413： 有限积分法和有限差分法 - ？
温艺双姜： 1.1 概念 有限差分方法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用.该方法将求解域划分为差分网格,用有限个网格节点代替连续的求解域.有限差分法以Taylor级数展开等方法,把控制方程中的导数用网格节点上的函数值...

冀州市18745612413： 化学.吉布 斯赫姆霍兹方程？
温艺双姜： 因为△H和△S随温度的变化很小,为了简化,计算时带入标准状态下的值.在物理化学课程中可以学习二者随温度而改变的方程.满意请采纳,谢谢^_^

冀州市18745612413： 国外注塑模具发展现状 - ？
温艺双姜： 随着全球经济的发展,新的技术革命不断取得新的进展和突破,技术的飞跃发展已经成为动世界经济增长的重要因素.市场经济的不断发展,促使工业产品越来越向多品种、小批量、高质量、低成本的方向发展,为了保持和加强产品在市场上的...

冀州市18745612413： “怎样理解数字图像处理的两个观点???解释一下.... 两个观点分别为“离散方法观点”、“连续方法观 - ？
温艺双姜： 离散 看成是模拟的. 连续看成是数字的