Явная схема трехмерной теплопроводности

[trofimovep]

Здравствуйте. Написал явную схему трехмерной теплопроводности самым простым методом сеток, с ГУ третьего рода.. Но решение находит неправильно, а именно: вывожу температурное поле в сечениях по х, у, 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
  

Кто может заметил, или знает, что не так, подскажите, пожалуйста.

[Karan]

Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Computer Science»

[GAA]

1. trofimovep, запишите исходную задачу (формулы с описаниями).
2. Обычно задают число узлов сетки (вдоль направления), и уже потом по этому числу находят длину шага (вдоль направления).
3. i, j, k изменяются в цикле от 1 до Nx, Ny, Nz. Проверка if i == A вызывает удивление.
4. Я бы не вставлял в основной цикл вычисление значений в граничных точках, а выполнил бы это отдельными циклами. (Для скорости вычислений.)
5. Я бы задавал размер массива при помощи zeros . Это позволит избежать изменения размера массива в цикле.