Boundary conditions for 2D heat transfer in a furnace.
1 view (last 30 days)
Show older comments
I've been working on a 2D heat transfer problem, I asked here and got a pretty good answer and actually managed to improve my code yesterday, but there's still an issue: There are 3 temperatures from T that are not updating at all. I think this is due to the internal and upper wall nodes that are somehow affecting the Vertical interior nodes and the bottom ones. I can say this because upon swapping the internal nodes' position in the code with the Bottom nodes I get a different result.
%%Parámetros geométricos del dominio
W=2; %Grosor (metros)
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-4;
Error=5;
max_iter=15000;
%%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 >= 2 && i <= Nx-1 %Pared superior (hinf y Tinf)
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 >= 2 && j <= Ny-1 %Nodos internos (inner nodes)
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>round(Nx*(4/7)) && i<Nx %Pared inferior (bottom)
if j==1
T(i,j)=(2*T(i,j+1)+T(i-1,j)+T(i+1,j)+(q_dot*dx^2)/K)/4;
end
elseif i == round(Nx * (4 / 7)) && j>=1 && j<= round(Ny*0.5) %Pared vertical interior (ho 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) + 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 > 1 && i < round(Nx * 4 / 7) %Pared horizontal interior (ho 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 == 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
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);
3 Comments
Torsten
on 17 Nov 2023
Edited: Torsten
on 17 Nov 2023
So 1/2 of the length in x-direction and 2/5 of the length in y-direction is cut off from the rectangular domain and boundary conditions are set at the vertical and horizontal walls ? What are these boundary conditions ? What are the boundary conditions at the other outer surfaces ? Don't you have a summary of your problem with geometric lengths, physical data, equation and boundary conditions ?
Answers (0)
See Also
Categories
Find more on Instrument Control Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!