Здравствуйте. Написал явную схему трехмерной теплопроводности самым простым методом сеток, с ГУ третьего рода.. Но решение находит неправильно, а именно: вывожу температурное поле в сечениях по х, у, z. И нигде нет совпадения температуры. Вот код в MATLAB:
Код:
clear all; clc
% Matetial
rho = 7800;
c = 500;
lambda = 46;
a = lambda / (rho * c); % thermal diffusitivy
% Parameters
A = 1.5; % i % height
B = 0.25; % j % width
C = 0.5 ; % k % length
t = 6; % p
% Bordering
T0 = 20;
hx = 0.01;
hy = 0.01;
hz = 0.1;
tau = 1 / (a * 30 * (1 / (hx * hx) + 1 / (hy * hy)) + 1 / (hz * hz)) ;
%hy = (A / B) * hx;
Nx = round(A / (hx));
Ny = round(B / (hy));
Nz = round(C / (hz));
Nt = round(t / tau);
%Tempearture of environment
Tl = 1300; % left
Tr = 1300; % right
Tu = 1300; % up
Td = 1300; % down
Tf = 1300; % front
Tb = 1300; % back
% flow
alpha_u = 1900;
alpha_l = 850;
alpha_r = 850;
alpha_d = 1900/5;
alpha_f = 850;
alpha_b = 850;
%
tic
for i = 1:Nx+1
for j = 1:Ny+1
for k = 1:Nz+1
T(i, j, k, 1) = T0;
end
end
end
for p = 1:Nt-1
for i = 1:1:Nx
for j = 1:1:Ny
for k =1:1:Nz
% left side
if i == 1
T(i, j, k, p+1) = T(i, j, k, p) + (alpha_l * (Tl - T(i, j, k, p))...
+ (lambda * (T(i+1, j, k, p) - T(i, j, k, p)) / hx)) * ...
(2 * tau / (c * rho * hx));
% right side
else if i == A
T(i, j, k, p+1) = T(i, j, k, p) + (alpha_r * (Tr - T(i, j, k, p))...
+ (lambda * (T(i-1, j, k, p) - T(i, j, k, p)) / hx)) * ...
(2 * tau / (c * rho * hx));
% down side
else if j == 1
T(i, j, k, p+1) = T(i, j, k, p) + (alpha_d * (Td - T(i, j, k, p))...
+(lambda * (T(i, j+1, k, p) - T(i, j, k, p)) / hx)) * ...
(2 * tau / (c * rho * hy));
% up side
else if j == B
T(i, j, k, p+1) = T(i, j, k, p) + (alpha_u * (Tu - T(i, j, k, p))...
+(lambda * (T(i, j-1, k, p) - T(i, j, k, p)) / hx)) * ...
(2 * tau / (c * rho * hy));
% front side
else if k == 1
T(i, j, k, p+1) = T(i, j, k, p) + (alpha_f * (Tf - T(i, j, k, p))...
+(lambda * (T(i, j, k+1, p) - T(i, j, k, p)) / hx)) * ...
(2 * tau / (c * rho * hz));
% back side
else if k == C
T(i, j, k, p+1) = T(i, j, k, p) + (alpha_b * (Tb - T(i, j, k, p))...
+(lambda * (T(i, j, k-1, p) - T(i, j, k, p)) / hz)) * ...
(2 * tau / (c * rho * hz));
% implement points
else
T(i , j, k, p+1) = T(i, j, k, p) +...
a * tau * (T(i+1, j, k, p) - 2 * T(i, j, k, p) + T(i-1, j, k, p))/hx^2 +...
a * tau * (T(i, j+1, k, p) - 2 * T(i, j, k, p) + T(i, j-1, k, p))/hy^2 +...
a * tau * (T(i, j, k+1, p) - 2 * T(i, j, k, p) + T(i, j, k-1, p))/hz^2;
end
end
end
end
end
end
end
end
end
end
% section in a yz plane
x = round(A / 2);
% section in a xz plane
y = 125;
% section in a xy plane
z = 3;
for i = 1:Nx+1
for j = 1:Ny
for k = 1:Nz
Txyend(i, j) = T(i, j, z, Nt);
Tyzend(j, k) = T(x, j, k, Nt);
Txzend(i, k) = T(i, 1, k, Nt);
end
end
end
pcolor(Txyend)
colorbar
grid on
s1 = ['section with z = ', z];
title(s1)
xlabel('x')
ylabel('y')
figure;
pcolor(Tyzend)
colorbar
grid on
s2 = ['section with x = ', x];
title(s2)
xlabel('y')
ylabel('z')
figure;
pcolor(Txzend)
colorbar
grid on
s3 = ['section with y = ', y];
title(s3)
xlabel('x')
ylabel('z')
toc
Кто может заметил, или знает, что не так, подскажите, пожалуйста.
