Array Values Not Changing Within Loop

Hi everyone,

I am working on a Runge Kutta code for modeling annular flow within a pipe.

I am able to plot graphs for my uR,T and Z as intended in the code, however my arrays for these values do not change from the initial values in the for loops so my graphs are simply straight lines at these values.

I believe there is something wrong with my loops.

Any help would be appreciated,

thanks.

 % Parameters & Boundary Conditions
dTdz=30.12; Ri=0.05; Ro=0.075; Rm=0.0621; hcoeff=100;  To=300; Ti=0; Zo=0; Zi=1.3*Ri*hcoeff*(Ti-To); uRi=0; uRo=0; tawi=4.3465; tawo=3.792; npoints=1000;
h=(Ro-Ri)/npoints;
r=(Ri+h:h:Ro)';
uR= zeros(npoints,1);
T= zeros(npoints,1);
Z = zeros(npoints,1);


%  Runge Kutta Function Fourth Order

for i=1:1:(npoints-1)
     T(1)=Ti;
     Z(1)=Zi;
     uR(1)= uRi;
        
     if r<Rm
         
         k1=h*dUdrI(r(i),uR(i));
         w1=h*dTdr(r(i),T(i));
         f1=h*dZdr(r(i),Z(i));
         
         k2=h*dUdrI((r(i)+h/2),(uR(i)+h*k1/2));
         w2=h*dTdr((r(i)+h/2),(T(i)+h*w1/2));
         f2=h*dZdr((r(i)+h/2),(Z(i)+h*f1/2));
         
         k3=h*dUdrI((r(i)+h/2),(uR(i)+h*k2/2));
         w3=h*dTdr((r(i)+h/2),(T(i)+h*w2/2));
         f3=h*dZdr((r(i)+h/2),(Z(i)+h*f2/2));
         
         k4=h*dUdrI((r(i)+h),(uR(i)+h*k3));
         w4=h*dTdr((r(i)+h),(T(i)+h*w3));
         f4=h*dZdr((r(i)+h),(Z(i)+h*f3));
         
         uR(i+1) = uR(i) + (k1+2*k2+2*k3+k4)/6;
         T(i+1) = T(i) + (w1+2*w2+2*w3+w4)/6;
         Z(i+1) = Z(i) + (f1+2*f2+2*f3+f4)/6;
      end
      
end
 
 for i=(npoints-1):1:1
        T(1)=To;
        Z(1)=Zo;
        uR(1)=uRo;
        if r>=Rm
            k1=h*dUdrO(r(i),uR(i));
            w1=h*dTdr(r(i),T(i));
            f1=h*dZdr(r(i),Z(i));
            
            k2=h*dUdrO((r(i)-h/2),(uR(i)-h*k1/2));
            w2=h*dTdr((r(i)-h/2),(T(i)-h*w1/2));
            f2=h*dZdr((r(i)-h/2),(Z(i)-h*f1/2));
            
            k3=h*dUdrO((r(i)-h/2),(uR(i)-h*k2/2));
            w3=h*dTdr((r(i)-h/2),(T(i)-h*w2/2));
            f3=h*dZdr((r(i)-h/2),(Z(i)-h*f2/2));
            
            k4=h*dUdrO((r(i)-h),(uR(i)-h*k3));
            w4=h*dTdr((r(i)-h),(T(i)-h*w3));
            f4=h*dZdr((r(i)-h),(Z(i)-h*f3));
            
            uR(i+1)=uR(i)+(k1+2*k2+2*k3+k4)/6;
            T(i+1)=T(i)+(w1+2*w2+2*w3+w4)/6;
            Z(i+1)=Z(i)+(f1+2*f2+2*f3+f4)/6;
            
        end
 end

 
figure, plot (r,uR)
xlabel('Radius')
ylabel('Velocity')
title('Velocity vs Radius')

figure, plot (r,T)
xlabel('Radius')
ylabel('Temperature')
title('Temperature vs Radius')

figure, plot (r,Z)
xlabel('Radius')
ylabel('Heat transfer per unit length')
title('Heat transfer per unit length vs Radius')

% Heat Transfer per Unit Length
function [dZdr_] = dZdr (r,rho,Cp,uR,dTdz)
dZdr_ = -r*rho(T)*Cp(T)*uR*dTdz;
end
 % Temperature Profile
function [dTdr_] = dTdr (Z,r,k,rho,epsilon_m)
dTdr_ = Z/(r*(k(T)+rho(T)*Cp(T)*epsilon_m(r,T,Rm,Ri,Ro,tawi,tawo)));
end 
% Velocity Profile
function [dUdrI_] = dUdrI (tawi,mu,rho,epsilon_m,Rm,r,Ri)
dUdrI_= tawi./(mu(T)+rho(T)*epsilon_m(r,T,Rm,Ri,Ro,tawi,tawo))*((Rm^2-r^2)/(Rm^2-Ri^2))*(Ri/r);
end

function [dUdrO_] = dUdrO (tawo,mu,rho,epsilon_m,r,Rm,Ro)
dUdrO_= tawo./(mu(T)+rho(T)*epsilon_m(r,T,Rm,Ri,Ro,tawi,tawo))*((r^2-Rm^2)/(Ro^2-Rm^2))*(Ro/r);
end

%eddy diffusivities of momentum
function [epsilon_m_] = epsilon_m (r,T,rm,ri,ro,tawi,tawo)
      if r<rm
        radratio=(r-ri)/(rm-ri);
        rEpsI=6.516e-4+3.9225e-1*radratio+(-6.0949e-1)*(radratio).^2+(2.3391e-1)*(radratio).^3+(1.410e-1)*(radratio).^4+(-9.6098e-2)*(radratio).^5;
        epsilon_m_=rEpsI*sqrt(tawi./rho(T))*(rm-ri);
      else
        radratio=(ro-r)/(ro-rm);
        rEpsO=6.516e-4+3.9225e-1*radratio+(-6.0949e-1)*(radratio).^2+(2.3391e-1)*(radratio).^3+(1.410e-1)*(radratio).^4+(-9.6098e-2)*(radratio).^5;
        epsilon_m_=rEpsO*sqrt(tawo./rho(T))*(ro-rm);
      end
end
     
function [Cp_]=Cp(T)
Cp_=1010.4755 + 0.1151*(T) + 4.00e-5*(T)^2;
end

function [k_]=k(T)
k_=0.0243+(6.548e-5)*(T) - (1.65e-8)*(T)^2;
end

%Dynamic Viscosity
function [mu_]=mu(T)
mu_=1.747e-5 + 4.404e-8*(T) - 1.645e-11*((T).^2);
end

function [Pr_]=Pr(T)
Pr_=0.7057*10^(2.06e-5*(T));
end

function [rho_]=rho (T)
rho_ =1e5/(287*(T));
end

function [vis_]=vis(T)
vis_=1.380e-5 + (8.955e-8)*(T) - (1.018e-10)*(T)^2;
end

[Sriram Tadavarty]

Hi Thomas,

This doesn't actually run in the for loops. so, the result of zeros all the time.

The if condition in the for loop must be replaced with if r(i) < Rm.

Then, the functions written dUdRi have different number of input arguments in the function call and the usage. Make sure all the calls to the different functions have same number of input arguments. That will help to see what other issues are in the code.

Regards,

Sriram

[Thomas Faulkner]

Hi Sriram,

thanks for your response.

I have tried changing the functions many times but keep getting 'not enough input arguments' for my dUdrI functoin.

i have attached a copy of my new code.

Kind regards,

Thomas

[Sriram Tadavarty]

It does complain because, function dUdrl has the following function signature

dUdrI (tawi,mu,rho,epsilon_m,Rm,r,Ri) % Implies there are 7 inputs

But you are passing only two inputs in the for loop code

k1=h*dUdrI(r(i),uR(i)); % Passing only r and uR

This is the issue for your case. You need to either update the function signature for velocity profile or send those many inputs to the function. I am not sure what is intended output with those functions.