三维热传导方程和泊松方程的有限元方法MATLAB程序

1. 三维热传导方程有限元解法

数学形式：

\(\frac{\partial u}{\partial t} = \alpha \nabla^2 u + f(x,y,z,t)\)

其中 \(\alpha\) 是热扩散系数。

%% 三维热传导方程有限元求解程序
% 使用线性四面体单元
clear; clc; close all;

%% 参数设置
alpha = 1.0;        % 热扩散系数
t_end = 1.0;        % 终止时间
dt = 0.01;          % 时间步长
nt = ceil(t_end/dt); % 时间步数

%% 生成简单的四面体网格（立方体域）
% 定义立方体范围
x_min = 0; x_max = 1;
y_min = 0; y_max = 1;
z_min = 0; z_max = 1;
nx = 5; ny = 5; nz = 5; % 每个方向的节点数

% 生成节点坐标
[x, y, z] = meshgrid(linspace(x_min, x_max, nx), ...
                     linspace(y_min, y_max, ny), ...
                     linspace(z_min, z_max, nz));
nodes = [x(:), y(:), z(:)];

% 生成四面体单元（简化方法：将每个立方体分成5个四面体）
tetra = [];
for i = 1:nx-1
    for j = 1:ny-1
        for k = 1:nz-1
            % 立方体的8个节点
            p1 = (i-1)*ny*nz + (j-1)*nz + k;
            p2 = (i-1)*ny*nz + (j-1)*nz + k+1;
            p3 = (i-1)*ny*nz + j*nz + k;
            p4 = (i-1)*ny*nz + j*nz + k+1;
            p5 = i*ny*nz + (j-1)*nz + k;
            p6 = i*ny*nz + (j-1)*nz + k+1;
            p7 = i*ny*nz + j*nz + k;
            p8 = i*ny*nz + j*nz + k+1;
            
            % 分割成立方体为5个四面体
            tetra = [tetra; p1 p2 p3 p6];
            tetra = [tetra; p1 p3 p4 p6];
            tetra = [tetra; p3 p4 p6 p7];
            tetra = [tetra; p4 p6 p7 p8];
            tetra = [tetra; p3 p6 p7 p5];
        end
    end
end

n_nodes = size(nodes, 1);
n_elem = size(tetra, 1);

%% 初始条件
u0 = zeros(n_nodes, 1);
for i = 1:n_nodes
    % 初始温度分布
    u0(i) = sin(pi*nodes(i,1)) * sin(pi*nodes(i,2)) * sin(pi*nodes(i,3));
end

%% 边界条件（Dirichlet边界）
bc_nodes = find(nodes(:,1)==x_min | nodes(:,1)==x_max | ...
                nodes(:,2)==y_min | nodes(:,2)==y_max | ...
                nodes(:,3)==z_min | nodes(:,3)==z_max);
bc_values = zeros(length(bc_nodes), 1);

%% 组装刚度矩阵和质量矩阵
K = sparse(n_nodes, n_nodes); % 刚度矩阵
M = sparse(n_nodes, n_nodes); % 质量矩阵
F = zeros(n_nodes, 1);        % 载荷向量

for e = 1:n_elem
    % 单元节点
    nodes_e = tetra(e, :);
    coords = nodes(nodes_e, :);
    
    % 计算单元刚度矩阵和质量矩阵
    [Ke, Me, Fe] = tetrahedral_element_matrices(coords, alpha);
    
    % 组装到全局矩阵
    K(nodes_e, nodes_e) = K(nodes_e, nodes_e) + Ke;
    M(nodes_e, nodes_e) = M(nodes_e, nodes_e) + Me;
    F(nodes_e) = F(nodes_e) + Fe;
end

%% 时间步进求解（使用θ方法，θ=0.5为Crank-Nicolson）
theta = 0.5;
u = u0;

% 预计算系数矩阵
A = M + theta * dt * K;
B = M - (1 - theta) * dt * K;

% 时间迭代
for t_step = 1:nt
    t = t_step * dt;
    
    % 更新载荷向量（考虑时间相关的源项）
    F_time = F * exp(-t); % 示例：指数衰减的源项
    
    % 形成右端项
    rhs = B * u + dt * ((1 - theta) * F_time + theta * F_time);
    
    % 处理Dirichlet边界条件
    A_fixed = A;
    rhs_fixed = rhs;
    
    for i = 1:length(bc_nodes)
        node_id = bc_nodes(i);
        A_fixed(node_id, :) = 0;
        A_fixed(node_id, node_id) = 1;
        rhs_fixed(node_id) = bc_values(i); % 边界值
    end
    
    % 求解线性系统
    u = A_fixed \ rhs_fixed;
    
    % 每10步输出一次
    if mod(t_step, 10) == 0
        fprintf('Time step: %d, Time: %.3f\n', t_step, t);
    end
end

%% 后处理和可视化
% 提取截面上的数据
z_slice = 0.5;
slice_nodes = find(abs(nodes(:,3) - z_slice) < 0.05);

figure;
scatter3(nodes(slice_nodes,1), nodes(slice_nodes,2), ...
         nodes(slice_nodes,3), 40, u(slice_nodes), 'filled');
colorbar;
title('三维热传导方程解（z=0.5截面）');
xlabel('X'); ylabel('Y'); zlabel('Z');

%% 四面体单元矩阵计算函数
function [Ke, Me, Fe] = tetrahedral_element_matrices(coords, alpha)
    % 计算四面体单元的刚度矩阵、质量矩阵和载荷向量
    % coords: 4×3矩阵，单元节点坐标
    
    % 单元体积计算
    A = [ones(4,1), coords];
    vol = abs(det(A)) / 6;
    
    % 自然坐标梯度
    dN_dxi = [-1, 1, 0, 0; 
              -1, 0, 1, 0;
              -1, 0, 0, 1]';
    
    % Jacobian矩阵
    J = dN_dxi' * coords;
    invJ = inv(J);
    
    % 形函数对物理坐标的梯度
    dN_dx = invJ' * dN_dxi';
    
    % 刚度矩阵 (∇N·∇N^T)
    Ke = alpha * vol * (dN_dx * dN_dx');
    
    % 质量矩阵 (N·N^T)
    % 对于线性四面体，质量矩阵的闭式表达式
    Me = vol/20 * (ones(4,4) + eye(4));
    
    % 载荷向量（假设单位源项）
    Fe = vol/4 * ones(4,1);
end

2. 三维泊松方程有限元解法

数学形式：

\(-\nabla^2 u = f(x,y,z)\)

带有边界条件。

%% 三维泊松方程有限元求解程序
clear; clc; close all;

%% 问题定义
% 泊松方程: -∇²u = f
% 定义源项函数
f = @(x,y,z) 3*pi^2 * sin(pi*x) .* sin(pi*y) .* sin(pi*z);

%% 生成网格（使用更规则的立方体网格）
Lx = 1; Ly = 1; Lz = 1;  % 区域尺寸
nx = 10; ny = 10; nz = 10; % 每个方向的单元数

% 创建六面体网格并转换为四面体
[x, y, z] = meshgrid(linspace(0, Lx, nx+1), ...
                     linspace(0, Ly, ny+1), ...
                     linspace(0, Lz, nz+1));

nodes = [x(:), y(:), z(:)];
n_nodes = size(nodes, 1);

% 生成四面体网格（简单方法）
tetra = [];
for i = 1:nx
    for j = 1:ny
        for k = 1:nz
            % 当前立方体的8个节点
            n1 = (i-1)*(ny+1)*(nz+1) + (j-1)*(nz+1) + k;
            n2 = n1 + 1;
            n3 = n1 + (nz+1);
            n4 = n3 + 1;
            n5 = n1 + (ny+1)*(nz+1);
            n6 = n2 + (ny+1)*(nz+1);
            n7 = n3 + (ny+1)*(nz+1);
            n8 = n4 + (ny+1)*(nz+1);
            
            % 将立方体分成6个四面体
            tetra = [tetra; n1 n2 n3 n6];
            tetra = [tetra; n1 n3 n4 n6];
            tetra = [tetra; n1 n5 n6 n3];
            tetra = [tetra; n3 n6 n7 n5];
            tetra = [tetra; n3 n4 n6 n7];
            tetra = [tetra; n4 n6 n7 n8];
        end
    end
end

n_elem = size(tetra, 1);

%% 边界条件（Dirichlet边界）
% 在边界上施加零边界条件
tol = 1e-10;
bc_nodes = find(abs(nodes(:,1)) < tol | abs(nodes(:,1)-Lx) < tol | ...
                abs(nodes(:,2)) < tol | abs(nodes(:,2)-Ly) < tol | ...
                abs(nodes(:,3)) < tol | abs(nodes(:,3)-Lz) < tol);
bc_values = zeros(size(bc_nodes));

%% 组装全局刚度矩阵和载荷向量
K_global = sparse(n_nodes, n_nodes);
F_global = zeros(n_nodes, 1);

% 计算每个单元的贡献
for e = 1:n_elem
    elem_nodes = tetra(e, :);
    coords = nodes(elem_nodes, :);
    
    [Ke, Fe] = poisson_element_matrices(coords, f);
    
    % 组装
    K_global(elem_nodes, elem_nodes) = K_global(elem_nodes, elem_nodes) + Ke;
    F_global(elem_nodes) = F_global(elem_nodes) + Fe;
end

%% 施加边界条件
K_fixed = K_global;
F_fixed = F_global;

% Dirichlet边界条件处理（置大数法）
for i = 1:length(bc_nodes)
    node_id = bc_nodes(i);
    K_fixed(node_id, :) = 0;
    K_fixed(node_id, node_id) = 1;
    F_fixed(node_id) = bc_values(i);
end

%% 求解线性系统
u = K_fixed \ F_fixed;

%% 计算精确解（用于验证）
u_exact = sin(pi*nodes(:,1)) .* sin(pi*nodes(:,2)) .* sin(pi*nodes(:,3));

%% 误差分析
error_L2 = sqrt(sum((u - u_exact).^2)/n_nodes);
fprintf('L2误差: %.6f\n', error_L2);
fprintf('最大误差: %.6f\n', max(abs(u - u_exact)));

%% 可视化结果
figure('Position', [100, 100, 1200, 400]);

% 子图1：数值解
subplot(1,3,1);
scatter3(nodes(:,1), nodes(:,2), nodes(:,3), 20, u, 'filled');
colorbar; title('数值解'); xlabel('X'); ylabel('Y'); zlabel('Z');

% 子图2：精确解
subplot(1,3,2);
scatter3(nodes(:,1), nodes(:,2), nodes(:,3), 20, u_exact, 'filled');
colorbar; title('精确解'); xlabel('X'); ylabel('Y'); zlabel('Z');

% 子图3：误差
subplot(1,3,3);
error = abs(u - u_exact);
scatter3(nodes(:,1), nodes(:,2), nodes(:,3), 20, error, 'filled');
colorbar; title('绝对误差'); xlabel('X'); ylabel('Y'); zlabel('Z');

%% 绘制截面图
figure;
z_idx = find(abs(nodes(:,3) - 0.5) < 0.01);
[~, sort_idx] = sort(nodes(z_idx,1));
plot(nodes(z_idx(sort_idx),1), u(z_idx(sort_idx)), 'b-', 'LineWidth', 2);
hold on;
plot(nodes(z_idx(sort_idx),1), u_exact(z_idx(sort_idx)), 'r--', 'LineWidth', 2);
legend('数值解', '精确解');
title('z=0.5截面上的解');
xlabel('X'); ylabel('u');

%% 泊松方程单元矩阵计算函数
function [Ke, Fe] = poisson_element_matrices(coords, f)
    % 计算四面体单元的刚度矩阵和载荷向量
    % coords: 4×3矩阵，单元节点坐标
    
    % 单元体积
    A = [ones(4,1), coords];
    vol = abs(det(A)) / 6;
    
    % 自然坐标梯度
    dN_dxi = [-1, 1, 0, 0; 
              -1, 0, 1, 0;
              -1, 0, 0, 1]';
    
    % Jacobian矩阵
    J = dN_dxi' * coords;
    invJ = inv(J);
    
    % 形函数对物理坐标的梯度
    dN_dx = invJ' * dN_dxi';
    
    % 刚度矩阵
    Ke = vol * (dN_dx * dN_dx');
    
    % 载荷向量（数值积分，使用中点规则）
    % 单元中心坐标
    xc = mean(coords(:,1));
    yc = mean(coords(:,2));
    zc = mean(coords(:,3));
    
    % 单元中心处的源项值
    f_val = f(xc, yc, zc);
    
    % 均匀分布载荷
    Fe = vol/4 * f_val * ones(4,1);
end

3. 高级功能：使用MATLAB PDE工具箱

%% 使用MATLAB PDE工具箱求解三维泊松方程
clear; clc; close all;

% 创建三维几何（单位球）
model = createpde();
importGeometry(model, 'Sphere.stl'); % 需要STL文件
% 或者使用内置几何体
% gm = multisphere(1);
% model.Geometry = gm;

% 生成网格
generateMesh(model, 'Hmax', 0.2);

% 查看网格
pdeplot3D(model);
title('三维有限元网格');

% 定义方程系数
% 对于泊松方程：-∇·(c∇u) + au = f
specifyCoefficients(model, 'm', 0, 'd', 0, 'c', 1, 'a', 0, 'f', 1);

% 边界条件（球表面u=0）
applyBoundaryCondition(model, 'dirichlet', 'Face', 1:model.Geometry.NumFaces, 'u', 0);

% 求解
results = solvepde(model);
u = results.NodalSolution;

% 可视化
figure;
pdeplot3D(model, 'ColorMapData', u);
title('三维泊松方程解（球域）');
colorbar;

4. 增强功能：自适应网格细化

%% 自适应网格细化求解泊松方程
function adaptive_poisson_3d()
    % 初始粗网格
    [nodes, tetra] = generate_initial_mesh();
    
    max_levels = 4;
    tolerance = 1e-3;
    
    for level = 1:max_levels
        % 求解当前网格
        [u, error_est] = solve_on_mesh(nodes, tetra);
        
        % 计算误差估计
        elem_error = estimate_error(nodes, tetra, u);
        
        % 标记需要细化的单元
        refine_flag = elem_error > tolerance * max(elem_error);
        
        if sum(refine_flag) == 0
            break;
        end
        
        % 细化网格
        [nodes, tetra] = refine_mesh(nodes, tetra, refine_flag);
        
        fprintf('Level %d: %d elements, max error: %.6f\n', ...
                level, size(tetra,1), max(elem_error));
    end
    
    % 最终解的可视化
    visualize_solution(nodes, tetra, u);
end

参考代码 对三维热传导方程和泊松方程的有限元方法解法matlab程序 www.3dddown.com/cnb/46697.html

说明：

1. 网格生成：三维有限元的关键是网格生成。上述代码使用了简单的结构化网格，实际应用中可能需要更复杂的网格生成器。

2. 单元类型：使用了线性四面体单元，这是三维有限元最常用的单元类型。

3. 边界条件：实现了Dirichlet边界条件处理。

4. 矩阵组装：利用MATLAB的稀疏矩阵存储提高效率。

5. 时间离散：对于热传导方程，使用了θ方法进行时间离散。

建议：

1. 对于复杂几何，建议使用专业的网格生成工具（如Gmsh）
2. 对于大规模问题，考虑使用迭代求解器（如PCG）代替直接求解器
3. 可以使用MATLAB的并行计算功能加速矩阵组装
4. 考虑使用pdetoolbox进行更复杂的有限元分析