Modelling Heat Transfer of water flow through a concrete duct, Can't get my water value to stay at 5 degrees

clear;

tic
k = 5; %Thermal conductivity of concrete - W/mK
T_w = 5; %Temperature of water - DegC
 %Temperature surface 
h = 100; %Convective heat transfer of water - W/m^2K
q_gen = 500; %Heat generation from top surface W/m^2
T_ts = 100; %Intial guess of top surface temperature - DegC
R = 5; %Grid resolution - mm
dx = 0.001*R; %Distance between nodes in x direction - m
dy = 0.001*R; %Distance between nodes in y direction - m
A = (100/R)+1;
B = (60/R)+1;
C = (30/R)+1;
D = (20/R)+1;
E = 5/R+1;
F = 35/R+1;
G = 50/R+1;
H = 25/R+1;
I = 75/R+1;
AA = 80/R +1;

T = ones(A,B)*T_ts;
Tp = T;
%x = 1;
i = 0;
error = 1;
max_error = 1e-6;
max_iter = 100000;

q = zeros(1,A);
Ts = zeros(1,A);
L = (0:R:100);
N = 1;
n=1;
m=1;

 



     while error>max_error && i<max_iter
         i = i + 1;
         Tp = T;

         for n = 1:A
            for m = 1:B
                
                
     %Water temperature 
         if (D<n)&&(n<AA)&&(C<m)&&(m<G)
                    
                   Tp(n,m) = T_w;
                    
 %Top left corner Node ("Node 1")
         elseif (m==1)&&(n==1)

    % T(n,m) = (Tp(m,n+1)+T_ts+Tp(m+1,n))/3 + (2*q_gen*dx*dx)/(3*k);
    T(n,m) = (Tp(n+1,m)+T_ts+Tp(n,m+1))/3 + (2*q_gen*dx*dx)/(3*k);
 
 %Top horizontal edge of concrete ("Node 2")
 elseif ((m==E)&&(m==B))&&(1==n)
 
     T(n,m) = (Tp(m-1,n))/6 + (T_ts+Tp(m,n-1)+Tp(m+1,n))/3 + (q_gen*dx^2)/(3*k);
   
 
 %Fully Enclosed Concrete Nodes ("Node ")
 elseif (1<m)&&(m<B)&&(1<n)&&(n<D) || (1<m)&&(m<B)&&(n<17)&&(n==A) || (1<n)&&(n<A)&&(1<m)&&(m<C) || (G<m)&&(m<B)&&(1<n)&&(n<A)


     %T(n,m) = (Tp(m-1,n)+ Tp(m+1,n)+Tp(m,n+1)+Tp(m,n-1))/4;
    T(n,m) = (Tp(n,m-1)+ Tp(n,m+1)+Tp(n+1,m)+Tp(n-1,m))/4;
 
 %Left vertical edge  %("Node ")
 elseif (n==1) && (1>m) && (m<A)

     %T(n,m) = (Tp(m,n+1)+Tp(m,n-1))/4 + (Tp(m+1,n))/2;
     T(n,m) = (Tp(n+1,m)+Tp(n-1,m))/4 + (Tp(n,m+1))/2;

 
 %Top left Corner of Water Node 59
    elseif (m==D)&&(n==C)
       
        %T(n,m)= (Tp(m-1,n)+Tp(m,n+1))/(3+(h*dx)/k)+ (Tp(m+1,n)+Tp(m,n-1))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
       T(n,m)= (Tp(n,m-1)+Tp(n+1,m))/(3+(h*dx)/k)+ (Tp(n,m+1)+Tp(n-1,m))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
 
       %Top Right Corner of Water
 %
     elseif (m==D)&&(n==G)
 
    % T(n,m) = (Tp(m+1,n)+Tp(m,n+1))/(3+(h*dx)/k)+ (Tp(m-1,n)+Tp(m,n-1))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
     T(n,m) = (Tp(n,m+1)+Tp(n+1,m))/(3+(h*dx)/k)+ (Tp(n,m-1)+Tp(n-1,m))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
 
 %Bottom Left Corner of Water

                 elseif (m==AA)&&(n==C)
 
     % T(n,m) = (Tp(m-1,n)+Tp(m,n-1))/(3+(h*dx)/k)+ (Tp(m+1,n)+Tp(m,n+1))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
    T(n,m) = (Tp(n,m-1)+Tp(n-1,m))/(3+(h*dx)/k)+ (Tp(n,m+1)+Tp(n+1,m))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);

 %Bottom Right Corner of Water Node 219
 
                 elseif (m==AA)&&(n==G)

                    % T(n,m) = (Tp(m+1,n)+Tp(m,n-1))/(3+(h*dx)/k)+ (Tp(m-1,n)+Tp(m,n+1))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
                     T(n,m) = (Tp(n,m+1)+Tp(n-1,m))/(3+(h*dx)/k)+ (Tp(n,m-1)+Tp(n+1,m))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);

   %Top of Water

                 elseif (n==D)&&(C<m)&&(m<G)

                         %T(n,m) = (Tp(m,n+1)*k + (Tp(m-1,n)+Tp(m+1,n))*(k/2-(h*dx)/2)+h*dx*(Tp(m,n-1)+T_w))/(2*k+h*dx);
                          T(n,m) = (Tp(n+1,m)*k + (Tp(n,m-1)+Tp(n,m+1))*(k/2-(h*dx)/2)+h*dx*(Tp(n-1,m)+T_w))/(2*k+h*dx);


 %Bottom of Water
 
                 elseif (n==AA)&&(C<m)&&(m<G) %(17, 8-10)
 
  %T(n,m) = (Tp(m,n-1)*k + (Tp(m-1,n)+Tp(m+1,n))*(k/2-(h*dx)/2)+h*dx*(Tp(m,n+1)+T_w))/(2*k+h*dx);

  T(n,m) = (Tp(n-1,m)*k + (Tp(n,m-1)+Tp(n,m+1))*(k/2-(h*dx)/2)+h*dx*(Tp(n+1,m)+T_w))/(2*k+h*dx);

 %Left side of Water

                 elseif (m==C)&&(D<n)&&(n<AA)

    
    T(n,m) = (Tp(n,m-1)*k + (Tp(n+1,m)+Tp(n-1,m))*(k/2-(h*dx)/2)+h*dx*(Tp(n,m+1)+T_w))/(2*k+h*dx);
% T(n,m) = (Tp(m-1,n)*k + (Tp(m,n+1)+Tp(m,n-1))*(k/2-(h*dx)/2)+h*dx*(Tp(m+1,n)+T_w))/(2*k+h*dx);
 %Right Side of Water
 
                 elseif (m==G)&&(D<n)&&(n<AA)

                   %T(n,m) = (Tp(m+1,n)*k + (Tp(m,n+1)+Tp(m,n-1))*(k/2-(h*dx)/2)+h*dx*(Tp(m-1,n)+T_w))/(2*k+h*dx);
                  T(n,m) = (Tp(n,m+1)*k + (Tp(n+1,m)+Tp(n-1,m))*(k/2-(h*dx)/2)+h*dx*(Tp(n,m-1)+T_w))/(2*k+h*dx);

 %Bottom Nodes Insulated

                 elseif (m==A)&&(1<n)&&(n==B)

                    %T(n,m) = (Tp(m-1,n)+Tp(m+1,n))/4 + (Tp(m,n+1))/2;
                    T(n,m) = (Tp(n,m-1)+Tp(n,m+1))/4 + (Tp(n+1,m))/2;

 %Bottom Left Corner Node
 elseif (m<1)&&(n==A)

    % T(n,m) = (Tp(m+1,n)+Tp(m,n+1))/2;
     T(n,m) = (Tp(n,m+1)+Tp(n+1,m))/2;

   
                   
                   end
                end
         end
%Temp(i, :, :) = T;
 %Error Calculation
error = max(abs(T-Tp),[],'all');
 it_error(i) = error;
 Tp = T;
 end

 %T = T';
% Tp = Tp';

 for N = 1:B
 q(N) = h*dx*(T(1,N)-T_w);
 
 Ts(N) = T(1,N);


end

q_total = sum(q);
Q = q_total/0.1;

For the water temp part I need the internal water nodes to be constant at 5 degrees, however it keeps reverting to the surface temp guess

Answer 1

[Udit06]

Hi Lee,

After analysing your code, I found out that you are updating the variable Tp instead of updating T due to which you are facing this issue.

You can make the following change to resolve the issue.

%Water temperature 
if (D<n)&&(n<AA)&&(C<m)&&(m<G)
    T(n,m) = T_w;
    
    %Rest of the code    
end

I hope this helps.