Finite method for a heat generating furnace

Hello! I'm currently working on a heat transfer problem, I don't consider myself to be an expert, therefore I'm still learning how to use this tool.

Here's the thing: I have pretty much everything working, but there's an issue. The ELSE equation is not giving any results and I don't know exactly why. I would appreciate if someone would be interested in helping me learn about this because I'm at a loss righ now

%% Parámetros geométricos del dominio
A=3;                        %Anchura (en metros)
H=2.5;                      %Altura (en metros)                        
Nx= 7;                      %Número de nodos en la dirección x
Ny=6;                       %Número de nodos en la dirección y
dx=A/(Nx-1);                %Distancia de elemento finito en dirección x (metros)
dy=H/(Ny-1);                %Distancia de elemento finito en direccion y (metros)

%% Condiciones iniciales y de frontera
To=1000;            %Temperatura interna del horno (en C)
T = ones(Nx,Ny);


Tinf=20;            %Temperatura externa (en C)

ho=100;             %Coeficiente de transferencia convectiva interior (en W/m^2K)
hinf=10;            %Coeficiente de transferencia convectiva exterior (en W/m^2K)

q_dot=8000;         %Calor generado por unidad volumétrica
K=1;                %Conductividad del material (en W/mK)

epsilon=1e-5;
Error=5;
%% Computación

iter=0;

while (Error>epsilon)
    iter=iter+1;
    disp(iter);

    Told=T;
   for j=1:Ny
        for i=1:Nx
            if i==Nx
                for j=2:Ny-1
                    T(i,j)=(2*T(i-1,j)+T(i,j-1)+T(i,j+1)+(q_dot*dx^2)/K)/4;
                end
                for j=Ny
                    T(i,j)=T(i,j-1)+T(i-1,j)+(q_dot*dx^2)/(2*K)+(hinf*dx*Tinf)/K;
                end
                for j=1
                    T(i,j)=((q_dot*dx^2)/(2*K)+T(i-1,j)+T(i,j+1));
                end
            elseif i==1
                for j=Ny*4/6:Ny-1
                    T(i,j)=(2*T(i+1,j)+T(i,j-1)+T(i,j+1)+(q_dot*dx^2)/K)/4;
                end
                for j=Ny*0.5
                    T(i,j)=T(i,j+1)+T(i+1,j)+(q_dot*dx^2)/(2*K)+(hinf*dx*Tinf)/K;
                end
                for j=Ny
                    T(i,j)=T(i,j-1)+T(i+1,j)+(q_dot*dx^2)/(2*K)+(ho*dx*To)/K;
                end
            elseif i>1 && i<Nx
                for j=Ny
                    T(i,j)=(T(i-1,j)+T(i+1,j)+2*T(i,j-1)+(q_dot*dx^2)/K+(2*hinf*dx*Tinf)/K)/(4+2*hinf*dx/K);
                end
            elseif 2<i && i<Nx*3/7
                for j=Ny*0.5
                    T(i,j)=((T(i-1,j)+T(i+1,j)+2*T(i,j+1)+(q_dot*dx^2)/K)+(2*ho*dx*To)/K)/(4+2*hinf*dx/K);
                end
            elseif i==Nx*4/7
                for j=2:Ny*2/6
                    T(i,j)=((T(i,j-1)+T(i-j+1)+2*T(i+1,j)+(q_dot*dx^2)/K)+(2*ho*dx*To)/K)/(4+2*hinf*dx/K);
                end            
                for j=1
                    T(i,j)=T(i,j-1)+T(i+1,j)+(q_dot*dx^2)/(2*K)+(ho*dx*To)/K;
                end
                for j=Ny*0.5
                    T(i,j)=(T(i,j+1)+2*T(i+1,j)+2*T(i,j+1)+T(i-1,j)+(3*q_dot*dx^2)/(2*K)+(2*ho*dx*To)/K)/(6+(2*ho*dx)/K);
                end
            else               
                T(i,j)=(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1)+(q_dot*dx^2)/K)/4;
            end
        end
    end
    
Error=sqrt(sumsqr(T-Told));
disp(Error);

end
T(1:3,1) = To; 
T(1:3,2)=To;                 
                
%% Graficando el resultado
x=0:dx:A;
y=0:dy:H;
% Definir los límites de color para resaltar variaciones
min_temp = min(min(T)); % Obtener el valor mínimo de temperatura
max_temp = max(max(T)); % Obtener el valor máximo de temperatura

colormap(hot); % Cambiar el colormap (puedes probar otros)
contourf(x, y, T', 'LineWidth', 1); % Añadir un grosor de línea para los contornos
colorbar;
title('Distribución de Temperatura')
xlabel('Dirección en x')
ylabel('Dirección en y')

% Establecer los límites de color para resaltar variaciones
clim([min_temp, max_temp]); % Establecer los límites de color según los valores mínimos y máximos de temperatura

Answer 1

[Torsten]

You use "j" twice as loop index in a nested loop. This will give problems.

for j=1:Ny

for i=1:Nx

if i==Nx

for j=2:Ny-1

[Santiago]

I just fixed it, and even MATLAB agrees with the decision, but I still won't let me get the code through. When I check the T matrix I still have some ones that haven't updated and apparently the index in position 1 of the 78th line is invalid. I think the main code is fine though.

%% Parámetros geométricos del dominio
A=3;                        %Anchura (en metros)
H=2.5;                      %Altura (en metros)                        
Nx= 7;                      %Número de nodos en la dirección x
Ny=6;                       %Número de nodos en la dirección y
dx=A/(Nx-1);                %Distancia de elemento finito en dirección x (metros)
dy=H/(Ny-1);                %Distancia de elemento finito en direccion y (metros)

%% Condiciones iniciales y de frontera
To=1000;            %Temperatura interna del horno (en C)
T = ones(Nx,Ny);


Tinf=20;            %Temperatura externa (en C)

ho=100;             %Coeficiente de transferencia convectiva interior (en W/m^2K)
hinf=10;            %Coeficiente de transferencia convectiva exterior (en W/m^2K)

q_dot=8000;         %Calor generado por unidad volumétrica
K=1;                %Conductividad del material (en W/mK)

epsilon=1e-5;
Error=5;
%% Computación

iter=0;

while (Error > epsilon)
    iter = iter + 1;
    disp(iter);

      Told = T;
    for j = 1:Ny
        for i = 1:Nx
            
            if i == Nx && 1<=j && Ny>=j                     %Pared derecha
                xx=1
                if j == Ny
                    T(i, j) = T(i, j - 1) + T(i - 1, j) + (q_dot * dx^2) / (2 * K) + (hinf * dx * Tinf) / K;       
                elseif j == 1
                    T(i, j) = ((q_dot * dx^2) / (2 * K) + T(i - 1, j) + T(i, j + 1));
                else 
                    T(i, j) = (2 * T(i - 1, j) + T(i, j - 1) + T(i, j + 1) + (q_dot * dx^2) / K) / 4;
                end
            elseif i == 1 && j>=0.5*Ny && j<=Ny             %Pared Izquierda
                xx=2
                if j == Ny * 0.5
                    T(i, j) = T(i, j + 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (hinf * dx * Tinf) / K;
                elseif j == Ny
                    T(i, j) = T(i, j - 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (ho * dx * To) / K;
                else 
                    T(i, j) = (2 * T(i + 1, j) + T(i, j - 1) + T(i, j + 1) + (q_dot * dx^2) / K) / 4;
                end
            elseif i > 1 && i < Nx                          %Pared superior (h y t inf)
                xx=3
                if j == Ny
                    T(i, j) = (T(i - 1, j) + T(i + 1, j) + 2 * T(i, j - 1) + (q_dot * dx^2) / K + (2 * hinf * dx * Tinf) / K) / (4 + 2 * hinf * dx / K);
                end
            elseif i > 1 && i < Nx * 3 / 7                  %Pared horizontal interior (h y T o)
                xx=4
                if j == Ny * 0.5
                    T(i, j) = ((T(i - 1, j) + T(i + 1, j) + 2 * T(i, j + 1) + (q_dot * dx^2) / K) + (2 * ho * dx * To) / K) / (4 + 2 * hinf * dx / K);
                end
            elseif i == Nx * (4 / 7) && j>=1 && j<=Ny* (2/6)  %Pared vertical interior (h y To)
                xx=5
                if j == 1
                    T(i, j) = T(i, j - 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (ho * dx * To) / K;
                elseif j == Ny * 0.5
                    T(i, j) = (T(i, j + 1) + 2 * T(i + 1, j) + 2 * T(i, j + 1) + T(i - 1, j) + (3 * q_dot * dx^2) / (2 * K) + (2 * ho * dx * To) / K) / (6 + (2 * ho * dx) / K);
                else   
                    T(i, j) = ((T(i, j - 1) + T(i - j + 1, j) + 2 * T(i + 1, j) + (q_dot * dx^2) / K) + (2 * ho * dx * To) / K) / (4 + 2 * hinf * dx / K);
                end

            elseif i>Nx*(4/7) && i<Nx                         %Pared  inferior
                xx=6
                if j==1
                T(i,j)=2*T(i,j+1)+T(i-1,j)+T(i+1,j)+(q_dot*dx^2)/K;
                end
                     
            else 
                 xx=7
                 i
                 j
                T(i, j) = (T(i-1, j) + T(i+1, j) + T(i, j-1) + T(i, j+1) + (q_dot * dx^2) / K) / 4;
               
            end
            
        end
      
    end     
            
        
    
    
    Error = sqrt(sumsqr(T - Told));
    disp(Error);

end  

     1

xx = 7

i = 1

j = 1

Index in position 1 is invalid. Array indices must be positive integers or logical values.

T(1:3,1) = To; 
T(1:3,2)=To;                 
                
%% Graficando el resultado
x = linspace(0, A, Nx);
y = linspace(0, H, Ny);

% Interpolación para obtener una visualización más suave
[X, Y] = meshgrid(x, y);
T_interp = interp2(X, Y, T', linspace(0, A, 10 * Nx), linspace(0, H, 10 * Ny)');

% Definir los límites de color para resaltar variaciones
min_temp = min(min(T_interp)); % Obtener el valor mínimo de temperatura
max_temp = max(max(T_interp)); % Obtener el valor máximo de temperatura

colormap(jet); % Cambiar el colormap (puedes probar otros)
pcolor(linspace(0, A, 10 * Nx), linspace(0, H, 10 * Ny), T_interp);
shading interp; % Interpolación para un mapa de colores suave
colorbar;
title('Distribución de Temperatura')
xlabel('Dirección en x')
ylabel('Dirección en y')

[Torsten]

If i =1 and j = 1, MATLAB enters the last "else" condition.

But to update

T(i, j) = (T(i-1, j) + T(i+1, j) + T(i, j-1) + T(i, j+1) + (q_dot * dx^2) / K) / 4;

both i and j have to be > 1 because else, T(i-1, j) and T(i, j-1) are not defined (i-1 = 0 and j-1 = 0).

[Santiago]

Okay, I've been trying all this time to change my code a bit and did some little improvements. Now I have a result but Matlab isn't updating the inner values. I've tried adding a for loop for the last equation, as its limits inside the mesh are i=2:Nx-1 and j=2:Ny-1 but matlab just doesn't give any logical answer. I feel like I'm on the right track but I just need to get T(i, j) = (T(i-1, j) + T(i+1, j) + T(i, j-1) + T(i, j+1) + (q_dot * dx^2) / K) / 4; to work. I've tried modifying its domain using an if condition without success.

%% Parámetros geométricos del dominio
A=3;                        %Anchura (en metros)
H=2.5;                      %Altura (en metros)                        
Nx= 7;                      %Número de nodos en la dirección x
Ny=6;                       %Número de nodos en la dirección y
dx=A/(Nx-1);                %Distancia de elemento finito en dirección x (metros)
dy=H/(Ny-1);                %Distancia de elemento finito en direccion y (metros)

%% Condiciones iniciales y de frontera
To=1000;            %Temperatura interna del horno (en C)
T = ones(Nx,Ny);


Tinf=20;            %Temperatura externa (en C)

ho=100;             %Coeficiente de transferencia convectiva interior (en W/m^2K)
hinf=10;            %Coeficiente de transferencia convectiva exterior (en W/m^2K)

q_dot=8000;         %Calor generado por unidad volumétrica
K=1;                %Conductividad del material (en W/mK)

epsilon=1e-5;
Error=epsilon+1;
max_iter=10000;
%% Computación

iter=0;
Told=T;
while (Error > epsilon && iter<max_iter)
    iter = iter + 1;
    disp(iter);
    for j = 1:Ny
        for i = 1:Nx
           
            if i == Nx && 1<=j && Ny>=j                     %Pared derecha
                if j == Ny
                    T(i, j) = T(i, j - 1) + T(i - 1, j) + (q_dot * dx^2) / (2 * K) + (hinf * dx * Tinf) / K;       
                elseif j == 1
                    T(i, j) = ((q_dot * dx^2) / (2 * K) + T(i - 1, j) + T(i, j + 1));
                else 
                    T(i, j) = (2 * T(i - 1, j) + T(i, j - 1) + T(i, j + 1) + (q_dot * dx^2) / K) / 4;
                end
            elseif i == 1 && j>=0.5*Ny && j<=Ny             %Pared Izquierda
                if j == Ny * 0.5
                    T(i, j) = T(i, j + 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (hinf * dx * Tinf) / K;
                elseif j == Ny
                    T(i, j) = T(i, j - 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (ho * dx * To) / K;
                else 
                    T(i, j) = (2 * T(i + 1, j) + T(i, j - 1) + T(i, j + 1) + (q_dot * dx^2) / K) / 4;
                end
            elseif i > 1 && i < Nx                          %Pared superior (h y t inf)
                if j == Ny
                    T(i, j) = (T(i - 1, j) + T(i + 1, j) + 2 * T(i, j - 1) + (q_dot * dx^2) / K + (2 * hinf * dx * Tinf) / K) / (4 + 2 * hinf * dx / K);
                end
            elseif i > 1 && i < Nx * 3 / 7                  %Pared horizontal interior (h y T o)
                if j == Ny * 0.5
                    T(i, j) = ((T(i - 1, j) + T(i + 1, j) + 2 * T(i, j + 1) + (q_dot * dx^2) / K) + (2 * ho * dx * To) / K) / (4 + 2 * hinf * dx / K);
                end
            elseif i == Nx * (4 / 7) && j>=1 && j<=Ny* (2/6)  %Pared vertical interior (h y To)
                if j == 1
                    T(i, j) = T(i, j - 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (ho * dx * To) / K;
                elseif j == Ny * 0.5
                    T(i, j) = (T(i, j + 1) + 2 * T(i + 1, j) + 2 * T(i, j + 1) + T(i - 1, j) + (3 * q_dot * dx^2) / (2 * K) + (2 * ho * dx * To) / K) / (6 + (2 * ho * dx) / K);
                else   
                    T(i, j) = ((T(i, j - 1) + T(i + 1, j) + 2 * T(i + 1, j) + (q_dot * dx^2) / K) + (2 * ho * dx * To) / K) / (4 + 2 * hinf * dx / K);
                end

            elseif i>Nx*(4/7) && i<Nx                         %Pared  inferior
                if j==1
                T(i,j)=2*T(i,j+1)+T(i-1,j)+T(i+1,j)+(q_dot*dx^2)/K;
                end            
            
                     
            else 
                if i>1 && j>1 && i<Nx && j<Ny
                T(i, j) = (T(i-1, j) + T(i+1, j) + T(i, j-1) + T(i, j+1) + (q_dot * dx^2) / K) / 4;
                        
                end
            end
           
            
        end
      
    end  
Error = sqrt(sumsqr(T - Told));
 disp(Error);
T(1:3,1)=To; 
T(1:3,2)=To;
Told = T;
      
end  

[Torsten]

Maybe it's time that you write down in a mathematical way what problem you are trying to solve.

[Santiago]

I just did, and I got a better code, but there's still a problem: Matlab doesn't calculate the values for the last equation. When you run the code you get 3 ones and I don't really know why.

%% Parámetros geométricos del dominio
A=3;                        %Anchura (en metros)
H=2.5;                      %Altura (en metros)                        
Nx= 7;                      %Número de nodos en la dirección x
Ny=6;                       %Número de nodos en la dirección y
dx=A/(Nx-1);                %Distancia de elemento finito en dirección x (metros)
dy=H/(Ny-1);                %Distancia de elemento finito en direccion y (metros)

%% Condiciones iniciales y de frontera
To=1000;            %Temperatura interna del horno (en C)
T = ones(Nx,Ny);


Tinf=20;            %Temperatura externa (en C)

ho=100;             %Coeficiente de transferencia convectiva interior (en W/m^2K)
hinf=10;            %Coeficiente de transferencia convectiva exterior (en W/m^2K)

q_dot=8000;         %Calor generado por unidad volumétrica
K=1;                %Conductividad del material (en W/mK)

epsilon=1e-5;
Error=epsilon+1;
max_iter=10000;
%% Computación

iter=0;
Told=T;
while (Error > epsilon && iter<max_iter)
    iter = iter + 1;
    %disp(iter);
    for j = 1:Ny
        for i = 1:Nx
           
            if i == Nx && 1<=j && Ny>=j                     %Pared derecha
                if j == Ny
                    T(i, j) = (T(i, j - 1) + T(i - 1, j) + (q_dot * dx^2) / (2 * K) + (hinf * dx * Tinf) / K)/(2+hinf*dx/K);       
                elseif j == 1
                    T(i, j) = ((q_dot * dx^2) / (2 * K) + T(i - 1, j) + T(i, j + 1));
                else 
                    T(i, j) = (2 * T(i - 1, j) + T(i, j - 1) + T(i, j + 1) + (q_dot * dx^2) / K) / 4;
                end
            elseif i == 1 && j>=round(0.5*Ny) && j<=Ny             %Pared Izquierda
                if j == round(Ny * 0.5)
                    T(i, j) = (T(i, j + 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (ho * dx * To) / K)/(2+ho*dx/K);
                elseif j == Ny
                    T(i, j) = (T(i, j - 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (hinf * dx * To) / K)/(2+hinf*dx/K);
                else 
                    T(i, j) = (2 * T(i + 1, j) + T(i, j - 1) + T(i, j + 1) + (q_dot * dx^2) / K) / 4;
                end
            elseif i > 1 && i < Nx                          %Pared superior (h y t inf)
                if j == Ny
                    T(i, j) = (T(i - 1, j) + T(i + 1, j) + 2 * T(i, j - 1) + (q_dot * dx^2) / K + (2 * hinf * dx * Tinf) / K) / (4 + 2 * hinf * dx / K);
                elseif j > 1 && j < Ny
                T(i, j) = (T(i - 1, j) + T(i + 1, j) + T(i, j - 1) + T(i, j + 1) - 4 * T(i, j)) * K / q_dot + T(i, j);  
                
                end
            elseif i > 1 && i < round(Nx * 4 / 7)                  %Pared horizontal interior (h y T o)
                if j == Ny * 0.5
                    T(i, j) = ((T(i - 1, j) + T(i + 1, j) + 2 * T(i, j + 1) + (q_dot * dx^2) / K) + (2 * ho * dx * To) / K) / (4 + 2 * ho * dx / K);
                end
            elseif i == round(Nx * (4 / 7)) && j>=1 && j<= round(Ny*0.5)  %Pared vertical interior (h y To)
                if j == 1
                    T(i, j) = (T(i, j + 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (ho * dx * To) / K)/(2+ho*dx/K);
                elseif j == round(Ny * 0.5)
                    T(i, j) = (T(i, j + 1) + 2 * T(i + 1, j) + 2 * T(i, j + 1) + T(i - 1, j) + (3 * q_dot * dx^2) / (2 * K) + (2 * ho * dx * To) / K) / (6 + (2 * ho * dx) / K);
                else   
                    T(i, j) = ((T(i, j - 1) + T(i + 1, j) + 2 * T(i + 1, j) + (q_dot * dx^2) / K) + (2 * ho * dx * To) / K) / (4 + 2 * ho * dx / K);
                end

            elseif i>round(Nx*(4/7)) && i<Nx && j==1                        %Pared  inferior
                
                T(i,j)=(2*T(i,j+1)+T(i-1,j)+T(i+1,j)+(q_dot*dx^2)/K)/4;

                          
  
            end
           
            
        end
      
    end  
Error = max(max(abs(T-Told)));
 %disp(Error);
T(1:3,1)=To; 
T(1:3,2)=To;
Told=T;
      
end  
disp(T);              

   1.0e+03 *

    1.0000    1.0000    1.0253    1.5195    1.5007    1.1218
    1.0000    1.0000    0.7977    0.7760    0.6808    0.3521
    1.0000    1.0000    0.5770    0.4111    0.3361    0.2463
    0.0010    0.3784    0.3298    0.2574    0.2296    0.2244
    0.0010    0.3210    0.3747    0.3435    0.2985    0.2364
    0.0010    0.7487    0.8388    0.7700    0.6045    0.2880
    3.1240    2.1230    1.8708    1.6826    1.3196    0.3868

                
%% Graficando el resultado
x = linspace(0, A, Nx);
y = linspace(0, H, Ny);

% Interpolación para obtener una visualización más suave
[X, Y] = meshgrid(x, y);
T_interp = interp2(X, Y, T', linspace(0, A, 10 * Nx), linspace(0, H, 10 * Ny)');

% Definir los límites de color para resaltar variaciones
min_temp = min(min(T_interp)); % Obtener el valor mínimo de temperatura
max_temp = max(max(T_interp)); % Obtener el valor máximo de temperatura

colormap(jet); % Cambiar el colormap (puedes probar otros)
pcolor(linspace(0, A, 10 * Nx), linspace(0, H, 10 * Ny), T_interp);
shading interp; % Interpolación para un mapa de colores suave
colorbar;
title('Distribución de Temperatura')
xlabel('Dirección en x')
ylabel('Dirección en y')