带热源的二维瞬态热传导问题

问：

有一个关于以热源为中心的二维热传导的问题。 它的形象是这样的。给定：2cm*2cm的热源可以看作是一个节点。t=0s时所有节点的温度是已知的，产生热通量q_generation也是已知的。 问题中的数字都不重要。 求：二维平面随时间变化的温度分布。

我已经知道没有热源的二维（逐节点）数值方法，但是如果有产生热量的热源，我不知道如何解决这个问题。 我应该从哪里开始计算？我编写了一个MATLAB代码，并得到了明显错误的答案。

format longG;
clc;clear;
%% Define Parameters and Initial Conditions
Lx = 0.3; % Length of the region
Ly = 0.2; % Width of the region
Nx = 31; % Number of points in each row/column
Ny = 21; % Number of points in each column/row
Tmax = 1; % Total time
dt = 10^(-2); % Time step
num_time_steps = Tmax/dt; % Total number of time steps
k = 14.4; % Thermal conductivity
h = 100; % Heat transfer coefficient
alpha = 0.387*10^(-5); % Thermal diffusivity
Cp = 461; % Specific heat capacity
den = 7817; % Density
q_gen = 50*10^(-6)*0.01*1;  % Heat source (50 watts per cubic centimeter)
num_i = 0;
num_j = 2;

dx = Lx / (Nx-1);
dy = Ly / (Ny-1);

%% Create an array to store temperature distribution at all time points
T = zeros(Ny, Nx, num_time_steps);
% Initial temperature distribution at t=0 (k=1)
T(:, :, 1) = 25*ones(Ny, Nx);

%% Time Stepping
% Loop to calculate temperature distribution at each time step
for t = 2:num_time_steps
    for i = (Ny+1)/2:Ny  % Horizontal from 11-21
        for j = (Nx+1)/2:Nx % Vertical from 16-31
            % Use the previous time step array to compute the next time step array
            if i ~= 11
                num_j = 2;
            end
            if i == (Ny+1)/2 && j == (Nx+1)/2
                % Calculate temperature at the corner
                T(i, j, t) = ((T(i-1, j, t-1)-2*T(i,j,t-1)+T(i+1,j,t-1))/dx^2 + (T(i, j+1, t-1)-2*T(i,j,t-1)+T(i,j-1,t-1))/dy^2 + q_gen/k)*alpha*dt+T(i,j,t-1);
            elseif i == Ny && j == Nx % Corner
                T(i, j, t) = ((0.5*k/dy^2+0.5*k/dx^2+h/(dy*2)+h/(dx*2))*-T(i, j, t-1)+...
                    (0.5*k/dy^2*T(i,j-1,t-1)+0.5*k/dx^2*T(i-1,j,t-1)+h/(dy*2)*25+...
                    h/(dx*2)*25))*dt+T(i,j,t-1);
                T(i-num_i,j,t)=T(i,j,t);
                T(i,j-num_j,t)=T(i,j,t);
                T(i-num_i,j-num_j,t)=T(i,j,t);
                message = sprintf('Corner has been processed, %d, %d', i, j);
                disp(message);
            elseif i == Ny && j ~= Nx % Rightmost column
                T(i, j, t) = ((alpha/dx^2+0.5*alpha/dy^2+0.5*alpha/dy^2+h/(dx*den*Cp))*-T(i, j, t-1)+...
                    (alpha/dx^2*T(i-1,j,t-1)+0.5*alpha/dy^2*T(i,j+1,t-1)+0.5*alpha/dy^2*T(i,j-1,t-1)+...
                    h/(dx*den*Cp)*25))*dt+T(i,j,t-1); 
                T(i,j-num_j,t)=T(i,j,t);
                message = sprintf('Right column has been processed, %d, %d', i, j);
                disp(message);
            elseif j == Nx && i ~= Ny % Bottommost row
                T(i, j, t) = ((alpha/dx^2+alpha/dy^2+h/(dy*den*Cp))*-T(i, j, t-1)+...
                    (0.5*alpha/dx^2*T(i-1,j,t-1)+0.5*alpha/dx^2*T(i+1,j,t-1)+alpha/dy^2*T(i,j-1,t-1)+...
                    h/(dy*den*Cp)*25))*dt+T(i,j,t-1);
                T(i-num_i, j, t) = T(i,j,t);
                message = sprintf('Bottom row has been processed, %d, %d', i, j);
                disp(message);
            else
                T(i, j, t) = ((T(i-1, j, t-1)-2*T(i,j,t-1)+T(i+1,j,t-1))/dx^2 + (T(i, j+1, t-1)-2*T(i,j,t-1)+T(i,j-1,t-1))/dy^2)*alpha*dt+T(i,j,t-1); % Update temperature with arbitrary values
                if i == 11
                    T(i,j-num_j,t)=T(i, j, t);
                    num_j = num_j + 2;
                else
                    T(i-num_i,j,t)=T(i, j, t);
                    T(i-num_i,j-num_j,t)=T(i,j,t);
                    T(i,j-num_j,t)=T(i,j,t);
                    num_j = num_j + 2;
                end
                message = sprintf('Inner region has been processed, %d, %d', i, j);
                disp(message);
            end
        end % Calculation for each row
        num_i = num_i + 2;
    end % Calculation for one time step

% Plot Temperature Distribution
    % Create a figure window
    figure(1);
    imagesc(T(:,:,t));
    colorbar;
    title(['Time = ' num2str(t*dt)]);
    xlabel('X');
    ylabel('Y');
    drawnow;

    num_i = 0;
    num_j = 2;
end