【无标题】

二维热传导方程–有限差分FAST格式解法

matlab 程序

clear, clc

%% 计算域与网格节点
Lx = 1;Ly = 1;             %指定计算域在二维平面xy上的范围
Nx = 20;Ny = Nx;           %指定单元网格数
dx = Lx/Nx;dy = Ly/Ny;     %x、y方向上的节点间距（m）
x = 0:dx:Lx; y = 0:dy:Ly;  %网格节点的位置
deltaT = 5;                %离散的时间步长（s）

%% 指定材料属性 
% K-热传导系数(W/mK)；rho-材料密度(kg/m^3);Cp-比热容(J/kgK)；
% alpha-热扩散系数(m^2/s);F-无量纲迭代系数。
k = 52.34;
rho = 7850;
Cp = 502; 
alpha = k/(rho*Cp);
F=alpha*deltaT/dx^2;

%% 边界条件
% T1-Neumann边界，隔热条件；T2-Dirichlet边界，指定顶部断面温度
T1 = 0;
T2 =3000;

%% 判敛准则
% ERR-指定两次迭代容许的最大误差；Lmax-最大迭代步长；
ERR=1e-4;
Lmax = 10^4;

%% 初始化温度场
T = T1*ones(Nx+1,Ny+1); % Initialize T array to T1 everywhere
T(2:Nx,Ny+1) = T2;      % Initialize top row to T2 boundary condition
T(1,Ny+1)    = T1;      % Initialize top left
T(Nx+1,Ny+1) = T1;      % Initialize top right
Tplot = ones(Lmax,1);   % Initialize Tplot to allocate memory

%% 求解器-有限差分法
for L = 1:Lmax
    Told = T;
    for j = 2:Ny 
        for i = 2:Nx
            T(i,j) = F*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1)...
                + Told(i,j+1)) + (1-4*F)*Told(i,j);
        end
    end
 % 设置监测点
    Tplot(L) = T(Nx/2,Ny/2);
    err= abs((T(Nx/2,Ny/2) - Told(Nx/2,Ny/2))/(T(Nx/2,Ny/2)));
    fprintf('Time step = %8.0f - err. = %10.6f percent\n', L, err);
    if (err < ERR)
        break;
    end
 % 设置计算云图
    Nc = 50; % Number of contours for plot
    dT = (T2 - T1)/Nc; % Temperature step between contours
    v = T1:dT:T2; % Sets temperature levels for contours
    colormap(jet) % Sets colors used for contour plot
    contourf(x, y, T',v, 'LineStyle', 'none') 
    colorbar                    % Adds a scale to the plot
    axis equal tight            % Makes the axes have equal length
    title(['Contour Plot of Temperature in deg. C at time = ',...
        num2str(deltaT*L/3600),' h']) 
    xlabel('x (m)') 
    ylabel('y (m)')        
    set(gca,'XTick',0:.1:Lx)    % Sets the x-axis tick mark locations
    set(gca,'YTick',0:.1:Ly)    % Sets the y-axis tick mark locations
    pause(0.001)                % Pause between time steps to display graph
    if L ==1||L==500||L==1000||L==2000||L==3000||L==4000 % Chosen time steps to save plot
       saveas(gcf, ['Transient_Plot_Unstable_',num2str(L)], 'jpg'); 
       save plot
    end
end

%% 迭代溢出提示
fprintf('Number of time steps = \t %8.0f \n\n', L)
if (L == Lmax)      % Warn if number of iterations exceeds maximum
    disp('Warning: Maximum time steps exceeded')   
    fprintf('\n') 
end

%% 创建图像
figure
plot(Tplot(1:L))
xlabel('Timestep')
ylabel('Temperature (C)') 
 for c=1:1:20
        disp(T(c,c));
    end