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

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，其坐标为（x_i，z_i），（x_j，z_j），（x_m，z_m），其函数值为u_i，u_j，u_m。除第一类边界条件给定的边界节点上的函数值是已知的外，其余节点上的函数值都是待求的。这样一来，将连续函数u（x，y）的求解化成节点上离散函数值的求解。

9.4.3 线性插值

由于各三角形单元取得足够小，使得可以假定各单元内电位随坐标是线性变化的，即u（x，z）是线性函数，有

u=α₁+α₂x+α₃z （9.4.3）

同时对三个顶点有

地球物理数据处理教程

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

地球物理数据处理教程

其中

地球物理数据处理教程

α_i=x_jz_m-x_mz_j

α_j=x_mz_i-x_iz_m

α_m=x_iz_j-x_jz_i

b_i=z_j-z_m

b_j=z_m-z_i

b_m=z_i-z_j

c_i=x_m-x_j

c_j=x_i-x_m

c_m=x_j-x_i

地球物理数据处理教程

将α₁、α₂和α₃代入（9.4.3）式中，整理可得

地球物理数据处理教程

或写为

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

式中：

地球物理数据处理教程

以上式中Δe为三角形单元的面积。

9.4.4 单元分析及总体合成

只研究（9.4.2）式对单元 e 的积分，记作 J_e。将（9.4.4）式代入（9.4.2）式中并沿整个单元积分，为求得使（9.4.2）式取极小的函数 u，对任意节点微分（9.4.2）式可得

地球物理数据处理教程

由（9.4.5）式有关系

地球物理数据处理教程

而

地球物理数据处理教程

再由（9.4.5）式有关系

地球物理数据处理教程

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

地球物理数据处理教程

最后可得

地球物理数据处理教程

同理可得

地球物理数据处理教程

或改写为

地球物理数据处理教程

地球物理数据处理教程

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

地球物理数据处理教程

由于k_rs=k_sr，所以k阵为对称矩阵，［u］_e是单元节点电位值组成的列向量。

如果节点r不属于单元e，则J_e（u）中不含有u_r，所以

地球物理数据处理教程

（矩阵中的虚点均为零元素）

扩展后的矩阵是l₀阶方阵，l₀是节点总数；扩展后的列向量 ｛u｝ 由全体节点函数值组成。

将所有单元J_e（u）合成，得到整个区域的J（u）

地球物理数据处理教程

它是所有节点函数值u₁……u_l0的函数，写成

J（u）=J（u₁……u_l0）

也可将J（u）看成变量u₁……u_l0的多元函数，所以泛函取极值相当于多元函数取极值

地球物理数据处理教程

图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是待求函数值的节点数，解这个方程组，即可求得这些节点上的函数值。