二维热传导方程–有限差分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