For loop doesn't work well.

I am trying to size a tube for a two-phase heat exchange problem.

I calculate the local pressure drop, I update the thermophysical properties via refrprop and then I calculate the heat exchange coefficient, that permits me to find finally the tube lenght.

The problem is that, a change on the steps number occurs in a change of results; it seems like the pressure drop is proportional to the numner of steps.

[SL: reformatted code as code not italicized text]

n = 500;
dx = (X_outlet_eva-X_inlet_eva)/(n);
x = zeros(n,1);
x(1) = X_inlet_eva;

P = zeros(n,1);
P(1) = Lowpressure;

for i = 2 : n+1
    
    x(i) = x(i-1) + dx;
    
    Viscosity_liquid(i) = refpropm('V','P',P(i-1),'Q',0,fluid); %[Pa s]
    Viscosity_vapor(i) = refpropm('V','P',P(i-1),'Q',1,fluid); %[Pa s]
    Density_liquid(i) = refpropm('D','P',P(i-1),'Q',0,fluid); %[kg/mc]
    Density_vapor(i) = refpropm('D','P',P(i-1),'Q',1,fluid); %[kg/mc]
    Enthalpy_liquid(i) = refpropm('H','P',P(i-1),'Q',0,fluid); %[J/kg]
    Enthalpy_vapor(i) = refpropm('H','P',P(i-1),'Q',1,fluid); %[J/kg]
    Latent(i) = (Enthalpy_vapor(i) - Enthalpy_liquid(i)); %[J/kg]
    Prandtl_liquid(i) = refpropm('^','P',P(i-1),'Q',0,fluid); %[-]
    Surfacetension_liquid(i) = refpropm('I','P',P(i-1),'Q',0,fluid); %[N/m]
    Thermalconductivity_liquid(i) = refpropm('L','P',P(i-1),'Q',0,fluid); %[W/mK]
    Specificheat_liquid(i) = refpropm('C','P',P(i-1),'Q',0,fluid); %[J/kgK]
    
    Ki_tt(i) = (((1-x(i))/x(i))^(0.9))*((Density_vapor(i)/Density_liquid(i))^(0.5))*((Viscosity_liquid(i)/Viscosity_vapor(i))^(0.1));
    
    Reynolds_liquid(i) = (G*(1-x(i))*Internalradius_eva*2)/Viscosity_liquid(i);
    Moody(i) = 0.316/(Reynolds_liquid(i)^(0.25));
    
    alfa_martinelli(i) = (1+0.28*Ki_tt(i)^0.71)^(-1);
    phi_nonhomo(i) =  (((1 - x(i))^2)/(1-alfa_martinelli(i)))+(Density_liquid(i)/Density_vapor(i))*((x(i)^2)/alfa_martinelli(i));
    dPD_acc(i) = (((x(i)^2)*G^2)/(alfa_martinelli(i)*Density_vapor(i))) + (((1-x(i))^2)*G^2)/(((1-alfa_martinelli(i))*Density_liquid(i))); %[Pa]
    dPD_fric(i) = phi_nonhomo(i)*Moody(i)*(L_tube(i-1)/(Internalradius_eva*2))*((G^2)/(2*Density_liquid(i))); %[Pa]
    dPD_tot(i) = (dPD_acc(i)+dPD_fric(i))*10^(-3); %[kPa]
    
    P(i) = P(i-1) - dPD_tot(i);
    
    T(i) = refpropm('T','P',P(i),'Q',0,fluid);
    
    % Convective heat exchangeo
    h_cb(i) = 0.023*(Thermalconductivity_liquid(i)/(Internalradius_eva*2))*(Reynolds_liquid(i)^(0.8))*(Prandtl_liquid(i)^(0.4)); %[W/m^2K]
    
    %Nucleate boiling heat transfer
    Twall_eva(i) = (Convectiveexchange_air*Temperature_air_inlet_eva*Externalradius_eva + h_chen(i-1)*T(i)*Internalradius_eva)/(Convectiveexchange_air*Externalradius_eva + h_chen(i-1)*Internalradius_eva);
    Psat_twall(i) = refpropm('P','T',Twall_eva(i),'Q',x(i),fluid);
    A_NB(i) = 0.00122*(Thermalconductivity_liquid(i)^(0.79))*(Specificheat_liquid(i)^(0.45))*(Density_liquid(i)^(0.49));
    B_NB(i) = (Twall_eva(i)-T(i))^(0.24);
    C_NB(i) = ((Psat_twall(i)-refpropm('P','T',T(i),'Q',0,fluid))*10^3)^(0.75);
    D_NB(i) = (Surfacetension_liquid(i)^(0.5))*(Viscosity_liquid(i)^(0.29))*(Latent(i)^(0.24))*(Density_vapor(i)^(0.24));
    h_nb(i) = (A_NB(i)/D_NB(i))*B_NB(i)*C_NB(i); %[W/m^2K]
    
    if (1/Ki_tt(i)) > 0.1
        F_chen(i) = 2.35*((0.213+(1/Ki_tt(i)))^(0.736));
    elseif (1/Ki_tt(i)) < 0.1
        F_chen(i) = 1;
    end
    
    S_chen(i) = 1/(1+((2.53*10^(-6))*(Reynolds_liquid(i)*(F_chen(i)^(1.25))^(1.17))));
    h_chen(i) = F_chen(i)*h_cb(i) + S_chen(i)*h_nb(i); %[W/m^2K]
    deltaT(i) = (Temperature_air_inlet_eva-T(i));
    U(i) = 1/((1/Convectiveexchange_air)+(1/h_chen(i)));
    Qv(i) = (Massflow_R32*Latent(i)*(x(i)-x(i-1))); %[W]
    L_tube(i) = Qv(i)/(U(i)*deltaT(i)*2*pi*Externalradius_eva);

[Steven Lord]

What does "doesn't work well" mean in this context?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support so we can investigate.

Answer 1

[Torsten]

Somewhere in the line

P(i) = P(i-1) - dPD_tot(i);

the "dx" must come into play because you usually compute the pressure gradient (P(i)-P(i-1))/dx in a point.

[Giovanni Karahoxha]

Thank you so much forte the answer, problem solved! I misread the equation and not taken care on the dimensions… I will take it in mind for future!