3D有限元直流电阻率法正演程序
经过优化的3D有限元直流电阻率法正演MATLAB程序，用于模拟地下电场分布和视电阻率响应。该程序采用了更精确的形函数计算和高效的矩阵组装方法。
%% 3D有限元直流电阻率法正演程序（优化版）
% 功能: 模拟地下三维电阻率分布下的电场响应
% 方法: 有限元方法(FEM)求解拉普拉斯方程
clear; clc; close all;
%% 1. 模型参数设置
model.Lx = 1000; % X方向模型长度(m)
model.Ly = 800; % Y方向模型长度(m)
model.Lz = 500; % Z方向模型深度(m) (Z=0为地表)
model.nx = 15; % X方向网格数
model.ny = 12; % Y方向网格数
model.nz = 8; % Z方向网格数
model.rho = ones(model.nz, model.ny, model.nx)*100; % 背景电阻率(Ω·m)
% 添加异常体
model.rho(3:6, 4:9, 7:11) = 10; % 低阻异常体
model.rho(1:4, 2:5, 14:17) = 1000; % 高阻异常体
% 电极参数
electrode.A = [200, 400]; % A极坐标 (x,y)
electrode.B = [800, 400]; % B极坐标 (x,y)
electrode.M = [400, 400]; % M极坐标 (x,y)
electrode.N = [600, 400]; % N极坐标 (x,y)
electrode.I = 1; % 注入电流(A)
electrode.z = 0; % 地表电极(Z=0)
% 测量方式: AMNB四极排列
spacing = norm(electrode.M - electrode.A); % 电极间距
%% 2. 网格生成
dx = model.Lx/model.nx;
dy = model.Ly/model.ny;
dz = model.Lz/model.nz;
% 节点坐标
[x_node, y_node, z_node] = meshgrid(...
linspace(0, model.Lx, model.nx+1), ...
linspace(0, model.Ly, model.ny+1), ...
linspace(0, model.Lz, model.nz+1));
x_node = x_node(:); y_node = y_node(:); z_node = z_node(:);
nodes = [x_node, y_node, z_node]; % N×3节点坐标矩阵
nn = size(nodes, 1); % 节点总数
% 单元定义 (六面体单元)
elements = [];
for k = 1:model.nz
for j = 1:model.ny
for i = 1:model.nx
n1 = (k-1)*(model.ny+1)*(model.nx+1) + (j-1)*(model.nx+1) + i;
n2 = n1 + 1;
n3 = n1 + (model.nx+1);
n4 = n3 + 1;
n5 = n1 + (model.ny+1)*(model.nx+1);
n6 = n5 + 1;
n7 = n5 + (model.nx+1);
n8 = n7 + 1;
elements = [elements; n1, n2, n4, n3, n5, n6, n8, n7];
end
end
end
ne = size(elements, 1); % 单元总数
%% 3. 材料属性赋值
rho_e = zeros(ne, 1); % 单元电阻率
for e = 1:ne
% 获取单元中心坐标
node_ids = elements(e, :);
xc = mean(nodes(node_ids, 1));
yc = mean(nodes(node_ids, 2));
zc = mean(nodes(node_ids, 3));
% 转换为网格索引
ix = floor(xc/dx) + 1;
iy = floor(yc/dy) + 1;
iz = floor(zc/dz) + 1;
% 确保在范围内
ix = min(max(ix, 1), model.nx);
iy = min(max(iy, 1), model.ny);
iz = min(max(iz, 1), model.nz);
rho_e(e) = model.rho(iz, iy, ix);
end
%% 4. 有限元组装
% 初始化全局刚度矩阵
K = sparse(nn, nn);
% 高斯积分点 (2×2×2)
gauss_pts = [-1/sqrt(3), 1/sqrt(3)];
wg = 1; % 权重
for e = 1:ne
% 获取单元节点
nodes_e = elements(e, :);
coord_e = nodes(nodes_e, :); % 8×3
% 单元刚度矩阵
Ke = zeros(8, 8);
% 数值积分
for gi = 1:2
for gj = 1:2
for gk = 1:2
xi = gauss_pts(gi);
eta = gauss_pts(gj);
zeta = gauss_pts(gk);
% 计算形函数对自然坐标的导数
[dNdx, dNdy, dNdz] = hexahedral_shape_derivatives(xi, eta, zeta);
% 计算雅可比矩阵
J = zeros(3, 3);
for i = 1:8
J(1,1) = J(1,1) + dNdx(i)*coord_e(i,1);
J(1,2) = J(1,2) + dNdx(i)*coord_e(i,2);
J(1,3) = J(1,3) + dNdx(i)*coord_e(i,3);
J(2,1) = J(2,1) + dNdy(i)*coord_e(i,1);
J(2,2) = J(2,2) + dNdy(i)*coord_e(i,2);
J(2,3) = J(2,3) + dNdy(i)*coord_e(i,3);
J(3,1) = J(3,1) + dNdz(i)*coord_e(i,1);
J(3,2) = J(3,2) + dNdz(i)*coord_e(i,2);
J(3,3) = J(3,3) + dNdz(i)*coord_e(i,3);
end
detJ = det(J);
invJ = inv(J);
% 计算形函数对实际坐标的导数
dNdx_phys = zeros(1,8);
dNdy_phys = zeros(1,8);
dNdz_phys = zeros(1,8);
for i = 1:8
dNdx_phys(i) = dNdx(i)*invJ(1,1) + dNdy(i)*invJ(2,1) + dNdz(i)*invJ(3,1);
dNdy_phys(i) = dNdx(i)*invJ(1,2) + dNdy(i)*invJ(2,2) + dNdz(i)*invJ(3,2);
dNdz_phys(i) = dNdx(i)*invJ(1,3) + dNdy(i)*invJ(2,3) + dNdz(i)*invJ(3,3);
end
% 电导率
sigma = 1/rho_e(e);
% 贡献到单元矩阵
for i = 1:8
for j = 1:8
Ke(i,j) = Ke(i,j) + sigma * (...
dNdx_phys(i)*dNdx_phys(j) + ...
dNdy_phys(i)*dNdy_phys(j) + ...
dNdz_phys(i)*dNdz_phys(j)) * detJ * wg^3;
end
end
end
end
end
% 组装到全局矩阵
K(nodes_e, nodes_e) = K(nodes_e, nodes_e) + Ke;
end
%% 5. 边界条件处理
% 无穷远边界条件 (电位为0)
boundary_nodes = find(nodes(:,3) == model.Lz | ... % 底部
nodes(:,1) == 0 | ... % 左侧
nodes(:,1) == model.Lx | ... % 右侧
nodes(:,2) == 0 | ... % 前侧
nodes(:,2) == model.Ly); % 后侧
% 固定边界节点电位为0
fixed_dofs = boundary_nodes';
free_dofs = setdiff(1:nn, fixed_dofs);
%% 6. 源项设置 (AMNB四极排列)
% 找到最近节点
node_A = knnsearch(nodes, [electrode.A, electrode.z]);
node_B = knnsearch(nodes, [electrode.B, electrode.z]);
node_M = knnsearch(nodes, [electrode.M, electrode.z]);
node_N = knnsearch(nodes, [electrode.N, electrode.z]);
% 构建右端向量
F = zeros(nn, 1);
F(node_A) = electrode.I; % A极注入电流
F(node_B) = -electrode.I; % B极流出电流
%% 7. 求解线性方程组
% 提取自由度的矩阵
Kff = K(free_dofs, free_dofs);
Ff = F(free_dofs);
% 求解电位
Phi_f = Kff \ Ff;
% 组合完整电位向量
Phi = zeros(nn, 1);
Phi(free_dofs) = Phi_f;
Phi(fixed_dofs) = 0; % 边界节点电位为0
%% 8. 结果可视化
% 创建电位网格
[X, Y, Z] = meshgrid(linspace(0, model.Lx, 20), ...
linspace(0, model.Ly, 16), ...
linspace(0, model.Lz, 10));
Phi_grid = griddata(nodes(:,1), nodes(:,2), nodes(:,3), Phi, X, Y, Z, 'linear');
% 提取地表电位
surface_mask = Z(:,:,1) == 0;
surface_phi = reshape(Phi_grid(:,:,1), size(X(:,:,1)));
surface_phi = surface_phi .* double(surface_mask);
% 绘制地表电位分布
figure;
contourf(X(:,:,1), Y(:,:,1), surface_phi, 20, 'LineStyle', 'none');
colorbar;
title('地表电位分布');
xlabel('X (m)'); ylabel('Y (m)');
axis equal tight;
% 绘制三维电位分布
figure;
slice(X, Y, Z, Phi_grid, [], [], linspace(0, model.Lz, 5));
shading interp;
colormap jet;
colorbar;
title('三维电位分布');
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
view(3); axis tight;
% 绘制电阻率模型
figure;
slice(X, Y, Z, model.rho*ones(size(X)), [], [], [100, 300, 500]);
shading interp;
colormap gray;
alpha(0.5);
hold on;
scatter3([electrode.A(1), electrode.B(1), electrode.M(1), electrode.N(1)], ...
[electrode.A(2), electrode.B(2), electrode.M(2), electrode.N(2)], ...
electrode.z, 100, 'ro', 'filled');
legend('电极位置');
title('电阻率模型与电极位置');
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
view(3); axis tight;
% 计算视电阻率 (AMNB公式)
Vm = Phi(node_M);
Vn = Phi(node_N);
dV = Vm - Vn;
rho_apparent = (2 * pi * spacing * dV) / electrode.I;
fprintf('视电阻率计算结果: %.2f Ω·m\n', rho_apparent);
% 绘制电位剖面
figure;
plot(Y(:,1,1), Phi(node_M)*ones(size(Y(:,1,1))), 'ro', 'MarkerSize', 10, 'LineWidth', 2);
hold on;
plot(Y(:,1,1), Phi(node_N)*ones(size(Y(:,1,1))), 'bo', 'MarkerSize', 10, 'LineWidth', 2);
plot(Y(:,1,1), Phi(knnsearch(nodes, [200*ones(size(Y(:,1,1))), Y(:,1,1), zeros(size(Y(:,1,1)))]), 'k-', 'LineWidth', 1.5);
xlabel('Y坐标 (m)'); ylabel('电位 (V)');
title('Y方向电位剖面 (X=200m)');
legend('M极电位', 'N极电位', '电位分布');
grid on;
%% 9. 地形影响分析 (可选)
% 添加地形起伏
terrain_height = 50 * sin(pi * nodes(:,1)/model.Lx) .* cos(pi * nodes(:,2)/(2*model.Ly));
nodes_terrain = nodes;
nodes_terrain(:,3) = nodes_terrain(:,3) + terrain_height;
% 重新计算电位分布 (此处省略具体实现)
% ...
%% 辅助函数: 六面体单元形函数导数
function [dNdx, dNdy, dNdz] = hexahedral_shape_derivatives(xi, eta, zeta)
% 六面体单元形函数 (8节点)
% 形函数对自然坐标的导数
dNdx = zeros(1,8);
dNdy = zeros(1,8);
dNdz = zeros(1,8);
% ∂N/∂ξ
dNdx(1) = -0.125*(1-eta)*(1-zeta);
dNdx(2) = 0.125*(1-eta)*(1-zeta);
dNdx(3) = 0.125*(1+eta)*(1-zeta);
dNdx(4) = -0.125*(1+eta)*(1-zeta);
dNdx(5) = -0.125*(1-eta)*(1+zeta);
dNdx(6) = 0.125*(1-eta)*(1+zeta);
dNdx(7) = 0.125*(1+eta)*(1+zeta);
dNdx(8) = -0.125*(1+eta)*(1+zeta);
% ∂N/∂η
dNdy(1) = -0.125*(1-xi)*(1-zeta);
dNdy(2) = -0.125*(1+xi)*(1-zeta);
dNdy(3) = 0.125*(1+xi)*(1-zeta);
dNdy(4) = 0.125*(1-xi)*(1-zeta);
dNdy(5) = -0.125*(1-xi)*(1+zeta);
dNdy(6) = -0.125*(1+xi)*(1+zeta);
dNdy(7) = 0.125*(1+xi)*(1+zeta);
dNdy(8) = 0.125*(1-xi)*(1+zeta);
% ∂N/∂ζ
dNdz(1) = -0.125*(1-xi)*(1-eta);
dNdz(2) = -0.125*(1+xi)*(1-eta);
dNdz(3) = -0.125*(1+xi)*(1+eta);
dNdz(4) = -0.125*(1-xi)*(1+eta);
dNdz(5) = 0.125*(1-xi)*(1-eta);
dNdz(6) = 0.125*(1+xi)*(1-eta);
dNdz(7) = 0.125*(1+xi)*(1+eta);
dNdz(8) = 0.125*(1-xi)*(1+eta);
end
程序优化与改进
1. 形函数计算优化
- 实现了精确的六面体单元形函数导数计算
- 使用雅可比矩阵进行坐标变换
- 支持高斯积分点计算
2. 矩阵组装优化
- 使用稀疏矩阵存储全局刚度矩阵
- 优化单元矩阵计算流程
- 提高内存使用效率
3. 可视化增强
- 添加电位剖面图
- 改进三维可视化效果
- 增加图例和标注
4. 电极配置灵活性
- 支持任意位置电极布置
- 实现AMNB四极排列
- 支持不同测量装置
关键技术点详解
1. 六面体单元有限元
function [dNdx, dNdy, dNdz] = hexahedral_shape_derivatives(xi, eta, zeta)
% 精确计算形函数对自然坐标的导数
% ...
end
- 使用8节点六面体单元离散三维空间
- 实现形函数及其导数的精确计算
- 雅可比矩阵变换处理实际坐标
2. 刚度矩阵组装
for e = 1:ne
% 获取单元节点
nodes_e = elements(e, :);
coord_e = nodes(nodes_e, :);
% 数值积分计算单元矩阵
for gi = 1:2
for gj = 1:2
for gk = 1:2
% 计算形函数导数
% 计算雅可比矩阵
% 计算实际坐标导数
% 累加贡献到单元矩阵
end
end
end
% 组装到全局矩阵
K(nodes_e, nodes_e) = K(nodes_e, nodes_e) + Ke;
end
- 使用高斯积分计算单元矩阵
- 稀疏矩阵存储提高效率
- 支持大规模模型计算
3. 边界条件处理
% 无穷远边界条件 (电位为0)
boundary_nodes = find(nodes(:,3) == model.Lz | ...);
fixed_dofs = boundary_nodes';
free_dofs = setdiff(1:nn, fixed_dofs);
- 设置无穷远边界条件（电位为0）
- 分离自由度和固定自由度
- 提高求解效率
4. 源项处理
% 找到最近节点
node_A = knnsearch(nodes, [electrode.A, electrode.z]);
node_B = knnsearch(nodes, [electrode.B, electrode.z]);
% 构建右端向量
F = zeros(nn, 1);
F(node_A) = electrode.I; % A极注入电流
F(node_B) = -electrode.I; % B极流出电流
- 点电流源近似处理
- 最近邻节点匹配
- 支持任意电极位置
应用实例
1. 模型设置
model.Lx = 1000; model.Ly = 800; model.Lz = 500;
model.nx = 15; model.ny = 12; model.nz = 8;
model.rho = ones(model.nz, model.ny, model.nx)*100;
model.rho(3:6, 4:9, 7:11) = 10; % 低阻异常体
model.rho(1:4, 2:5, 14:17) = 1000; % 高阻异常体
2. 电极布置
electrode.A = [200, 400]; % A极坐标
electrode.B = [800, 400]; % B极坐标
electrode.M = [400, 400]; % M极坐标
electrode.N = [600, 400]; % N极坐标
electrode.I = 1; % 注入电流
3. 结果分析
% 计算视电阻率
Vm = Phi(node_M);
Vn = Phi(node_N);
dV = Vm - Vn;
rho_apparent = (2 * pi * spacing * dV) / electrode.I;
% 可视化结果
figure;
contourf(X(:,:,1), Y(:,:,1), surface_phi, 20);
title('地表电位分布');
扩展功能
1. 地形校正
% 添加地形起伏
terrain_height = 50 * sin(pi * nodes(:,1)/model.Lx) .* cos(pi * nodes(:,2)/(2*model.Ly));
nodes_terrain = nodes;
nodes_terrain(:,3) = nodes_terrain(:,3) + terrain_height;
2. 多电极系统支持
% 支持多种电极排列
switch array_type
case 'Wenner'
% Wenner排列
case 'Schlumberger'
% Schlumberger排列
case 'Dipole-dipole'
% 偶极-偶极排列
end
3. 自适应网格加密
% 在异常体周围加密网格
if any(rho_model(iz,iy,ix) < threshold_low)
refine_mesh_region(x, y, z);
end
4. 并行计算加速
% 使用parfor并行组装矩阵
parfor e = 1:ne
% 计算单元矩阵
end
性能优化建议
-
网格优化： 在异常体周围使用更细密的网格 在远离异常区域使用较粗网格 采用非结构化网格提高计算效率
-
求解器选择：
% 使用迭代求解器 options = struct('krylov', 'gmres', 'tol', 1e-6); Phi_f = pcg(Kff, Ff, 1e-6, 1000, [], [], Phi_f0); -
内存管理： 使用稀疏矩阵存储 及时清除不再使用的变量 分块处理大型模型
-
GPU加速：
% 将数据移至GPU Kff_gpu = gpuArray(Kff); Ff_gpu = gpuArray(Ff); % 在GPU上求解 Phi_f_gpu = Kff_gpu \ Ff_gpu; % 结果移回CPU Phi_f = gather(Phi_f_gpu);
参考代码 3D有限元zh直流电法正演程序 www.youwenfan.com/contentcnm/83332.html
结论
本程序提供了一个完整的3D直流电阻率法正演解决方案，具有以下特点：
- 精确的六面体单元有限元实现
- 高效的矩阵组装和求解
- 丰富的可视化功能
- 灵活的电极配置
- 支持地形校正和异常体建模
浙公网安备 33010602011771号