Hi I am working on a program to calculate the heat lost by a pipe while water is traveling through it.
Here is the code
clearWaterTemperatureCelcius = [0.01 10:10:100];WaterTemperature = WaterTemperatureCelcius+273.15;WaterDensityMatrix= [916.8 999.8 998.3 995.7 992.3 988 983 978 972 965 958];WaterDynamicViscosityMatrix = [1.78 1.31 1.00 0.798 0.653 0.547 0.467 0.404 0.355 0.314 0.281]; WaterKinematicViscosityMatrix = [1.792 1.304 1.004 0.801 0.658 0.553 0.474 0.413 0.365 0.326 0.295]; KWaterMatrix = [-0.07 0.088 0.207 0.303 0.385 0.457 0.523 0.585 0.643 0.665 0.752];PrWaterMatrix = [13.67 9.47 7.01 5.43 4.34 3.56 2.99 2.56 2.23 1.96 1.75];CWaterMatrix = [4.21 4.193 4.183 4.179 4.179 4.182 4.185 4.191 4.198 4.208 4.219]; AirTemperature = 100:50:500;AirDynamicViscosityMatrix = [0.6924 1.0283 1.3289 1.488 1.893 2.075 2.286 2.484 2.671]; AirKinematicViscosityMatrix = [1.923 4.343 4.49 9.49 15.68 20.76 25.9 28.86 37.9];AirDensityMatrix = [3.601 2.3675 1.7684 1.4128 1.774 0.998 0.8826 0.7833 0.7048];KAirMatrix = [0.009246 0.013735 0.01809 0.02227 0.02624 0.03003 0.03365 0.03707 0.04038]; PrAirMatrix = [0.770 0.753 0.739 0.722 0.708 0.697 0.689 0.683 0.680];PipeTemperatureMatrixCelcius = [-100 0 100 200 300 400 600 800 1000 1200];PipeTemperatureMatrix = PipeTemperatureMatrixCelcius+273.15;KPipeMatrix = [87 73 67 62 55 48 40 36 35 36];Two = input ('Guess the outside temperature of the pipe wall (K): ');Twi = input ('Guess the inside temperature of the pipe wall (K): ');Tm = input ('Guess the temperature in the middle of the flow (K): ');Do = input ('Please enter the outside diameter of the pipe (m): ');Di = input ('Please enter the inside diameter of the pupe (m): ');Flow = input ('Please enter 1 for laminar flow and 2 for turbulent flow: ');Tinf = 266.48; Vinf = 6.7056; Tin = 372.04; L = 30.48; dx = 0.1; ivalue = L/dx;Pi = 3.141593;TTest1 = 1;TTest2 = 1;TTest3 = 1;Tx = Tin;Ai = Pi*Di*dx;Ao = Pi*Do*dx;for i=1:ivalue, if i == 1 Tx=Tin; end if i > 1 Tx = Tout; end while (abs(TTest1-Two) > 0.001) && (abs(TTest2-Twi) > 0.001) && (abs(TTest3-Tm) > 0.001) disp(['Two is ' num2str(Two)]) disp(['Twi is ' num2str(Twi)]) disp(['Tm is ' num2str(Tm)]) Tfo = (Two+Tinf)/2; Tfi = (Twi+Tin)/2; PipeTemperature = Two-Twi; WaterKinematicViscosity = (interp1(WaterTemperature, WaterKinematicViscosityMatrix, Tfi))*(10^(-6)); WaterDynamicViscosity = (interp1(WaterTemperature, WaterDynamicViscosityMatrix, Tfi))*(10^(-2)); WaterDensity = interp1(WaterTemperature, WaterDensityMatrix, Tfi); KWater = interp1(WaterTemperature, KWaterMatrix, Tfi); PrWater = interp1(WaterTemperature, PrWaterMatrix, Tfi); CWater = (interp1(WaterTemperature, CWaterMatrix, Tfi)*1000); AirKinematicViscosity = (interp1(AirTemperature, AirKinematicViscosityMatrix, Tfo))*(10^(-6)); AirDynamicViscosity = (interp1(AirTemperature, AirDynamicViscosityMatrix, Tfo))*(10^(-5)); AirDensity = interp1(AirTemperature, AirDensityMatrix, Tfo); KAir = interp1(AirTemperature, KAirMatrix, Tfo); PrAir = interp1(AirTemperature, PrAirMatrix, Tfo); KPipe = interp1(PipeTemperatureMatrix, KPipeMatrix, PipeTemperature); if Flow==1 ReDi = 2000; hi = (4.363*KWater/Di); else ReDi = 5000; hi= (KWater*0.023*(ReDi^0.8)*(PrWater^0.3)); end Vm = (ReDi*WaterKinematicViscosity)/Di; ReDf = Vinf*Do/AirKinematicViscosity; disp(['ReDf is ' num2str(ReDf)]) disp(['ReDi is ' num2str(ReDi)]) disp(['hi is ' num2str(hi)]) if 0.4 < ReDf <= 4 c=0.989; n=0.33; elseif 4 < ReDf <= 40 c=0.911; n=0.385; elseif 40 < ReDf <= 4000 c=0.683; n=0.466; elseif 4000 < ReDf <= 40000 c=0.193; n=0.618; elseif 40000 < ReDf <= 400000 c=0.0266; n=0.805; end disp(['c is ' num2str(c)]) disp(['n is ' num2str(n)]) NuDf = c*(ReDf^n)*(PrAir^(1/3)); disp(['NuDf is ' num2str(NuDf)]) ho = NuDf*KAir/Do; MassFlowRate = WaterDensity*Vm*Pi*Di*Di/4; disp(['Mass flow rate is ' num2str(MassFlowRate) ' kg/s']) q = (Tx-Tinf)/(((1/(ho*Ao))+((1/(2*Pi*KPipe*dx))*log(Do/Di))+(1/(hi*Ai))+(1/(2*MassFlowRate*CWater)))); TwoNew = ((q/(ho*Ao))+Tinf); TwiNew = (q*(log(Do/Di))/(2*Pi*KPipe*dx))+TwoNew; TmNew = Tx-(q/(2*MassFlowRate*CWater)); disp(['TwoNew is ' num2str(TwoNew)]) disp(['TwiNew is ' num2str(TwiNew)]) disp(['TmNew is ' num2str(TmNew)]) TTest1 = Two; TTest2 = Twi; TTest3 = Tm; Two = TwoNew; Twi = TwiNew; Tm = TmNew; end Txdx = (2*TmNew) - Tx; Tout = Tx + Txdx; endq1 = ho*(Twi-Tinf)*Ao;q2 = (Twi-Two)/((1/(2*Pi*KPipe*dx))*log(Do/Di));q3 = hi*(Tm-Twi)*Ai;q4 = MassFlowRate*CWater*(Tx-((2*Tm)-Tx));disp (['The value of q1 is ' num2str(q1) ' W'])disp (['The value of q2 is ' num2str(q2) ' W'])disp (['The value of q3 is ' num2str(q3) ' W'])disp (['The value of q4 is ' num2str(q4) ' W'])disp (['The value of q5 is ' num2str(q) ' W'])I keep getting NaN output for the temperature, which stops the loop
Guess the outside temperature of the pipe wall (K): 320 Guess the inside temperature of the pipe wall (K): 330 Guess the temperature in the middle of the flow (K): 340 Please enter the outside diameter of the pipe (m): 1 Please enter the inside diameter of the pupe (m): 0.9 Please enter 1 for laminar flow and 2 for turbulent flow: 1 Two is 320 Twi is 330 Tm is 340 ReDf is 451765.1016 ReDi is 2000 hi is 3.0572 c is 0.989 n is 0.33 NuDf is 64.8209 Mass flow rate is 0.51629 kg/s TwoNew is NaN TwiNew is NaN TmNew is NaN The value of q1 is NaN W The value of q2 is NaN W The value of q3 is NaN W The value of q4 is NaN W The value of q5 is NaN WAny advice on solving this would help.
Thanks.
Best Answer