二维非稳态热传导问题

非稳态二维热传导方程为

$$\frac{\partial T}{\partial t}=\alpha \left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right)$$
上述方程是针对边缘保持恒定温度的方形板求解的。换句话说，指定了狄尔奇雷边界条件。使用有限差分法对二阶偏微分方程进行离散化。使用显式和 ADI 方法求解离散化的 PDE。

% Parameters
thermal_diffusivity = 0.835;  % Thermal diffusivity (alpha)
L = 10;                       % Domain length (m)
t_start = 0;                  % Start time
t_end = 100;                  % End time
x_start = 0;                  % x start
x_end = 40;                   % x end
dt = 10;                      % Time step
dx = 10;                      % Spatial step (uniform in x and y)
lambda = (thermal_diffusivity * dt) / dx^2;  % Stability parameter

% Grid initialization
nx = (x_end - x_start)/dx + 1;  % Number of x nodes
ny = nx;                        % Assume square domain
nt = (t_end - t_start)/dt + 1;  % Number of time steps
T = zeros(nx, ny, nt);          % Temperature field

% Boundary conditions
T(:, 1, :) = 75;      % Left boundary
T(:, end, :) = 50;    % Right boundary
T(1, :, :) = 0;       % Top boundary
T(end, :, :) = 100;   % Bottom boundary

% ADI Method
for n = 2:nt
    % Step 1: X-direction implicit (solve for T_half)
    T_half = T(:, :, n-1);
    for i = 2:nx-1
        % Construct tridiagonal system
        e = -lambda * ones(nx-2, 1);  % Lower diagonal
        f = 2*(1 + lambda) * ones(nx-1, 1);  % Main diagonal
        g = -lambda * ones(nx-2, 1);  % Upper diagonal
        R = zeros(nx-1, 1);
        
        % Fill R vector
        for j = 2:ny-1
            R(j-1) = lambda*T_half(j, i-1) + 2*(1-lambda)*T_half(j, i) + lambda*T_half(j, i+1);
        end
        R(1) = R(1) + lambda * T(1, i, n-1);      % Top BC
        R(end) = R(end) + lambda * T(end, i, n-1); % Bottom BC
        
        % Solve and update
        T_half(2:end-1, i) = ThomasAlgorithm(e, f, g, R);
    end
    T(:, :, n) = T_half;
    
    % Step 2: Y-direction implicit (solve for full step)
    for j = ny-1:-1:2
        e = -lambda * ones(ny-2, 1);
        f = 2*(1 + lambda) * ones(ny-1, 1);
        g = -lambda * ones(ny-2, 1);
        R = zeros(ny-1, 1);
        
        for i = 2:nx-1
            R(i-1) = lambda*T(j-1, i, n) + 2*(1-lambda)*T(j, i, n) + lambda*T(j+1, i, n);
        end
        R(1) = R(1) + lambda * T_half(j, 1);      % Left BC
        R(end) = R(end) + lambda * T_half(j, end); % Right BC
        
        T(j, 2:end-1, n) = transpose(ThomasAlgorithm(e, f, g, R));
    end
end

% Visualization
[x, y] = meshgrid(x_start:dx:x_end, x_start:dx:x_end);
T_last = T(:, :, end);
figure;
mesh(x, y, T_last);
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Temperature (°C)');
title('ADI Solution for 2D Heat Equation');
colorbar;

% Thomas Algorithm (optimized)
function X = ThomasAlgorithm(e, f, g, R)
    n = length(f);
    % Decomposition
    for i = 2:n
        e(i-1) = e(i-1) / f(i-1);
        f(i) = f(i) - e(i-1) * g(i-1);
    end
    % Forward substitution
    for i = 2:n
        R(i) = R(i) - R(i-1) * e(i-1);
    end
    % Backward substitution
    X = zeros(n, 1);
    X(n) = R(n) / f(n);
    for i = n-1:-1:1
        X(i) = (R(i) - g(i) * X(i+1)) / f(i);
    end
end