How to resolve error using erfinv: Input must be real and full?
2 次查看(过去 30 天)
显示 更早的评论
The error that I received from the following code is as follows:
Error using erfinv
Input must be real and full.
Error in InletJetMixing (line 173)
tnew(i,1) = (1/(4*alfa))*(Wth(i,1)/(2*erfinv(1 + 0.001*(Tc/(Tc-TaveS(i+1,1))))))^2;
That would be appreciated if you could help in resolving it.
% Model Inputs
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------
Th = 50; % Hot Temperature (C)
Tc = 20; % Cold Temperature (C)
Tave = (Th+Tc)/2; % Average Temperature (C)
mf = 0.001; % Mass Flow Rate (kg/s)
H = 1; % Tank Height (m)
Dtank = 0.3; % Tank Diameter (m)
Dpipe = 0.01; % Inlet Pipe Diameter (m)
visc = 0.00070057; % Average Fluid Viscosity (m2/s)
rhoH = 987.68; % Hot Water Density (kg/m3)
rhoC = 997.78; % Cold Water Density (kg/m3)
rho = 0.5*(rhoH + rhoC); % Average Fluid Density (kg/m3)
k = 0.62614; % Thermal Conductivity (W/mK)
cp = 4068.5; % Specific Heat Capacity (J/kgK)
beta = 0.00032452; % Thermal Expansion Coeff. (1/K)
alfa = k/(rho*cp); % Diffusivity (m2/s)
N = 100; % Number of Nodes
delt = 1; % Time Step Size (s)
tTotal = 7200; % Total Simulation Time (s)
g = 9.81; % Gravitational Acceleration (m/s2)
% Calculate Model Variables
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------
Q = mf/rho; % Volumetric Flow Rate (m3/s)
Ac = (pi/4)*(Dtank^2); % Tank Cross-section area (m2)
u = mf/(rho*(pi/4)*(Dtank^2)); % Mean Flow Velocity (m/s)
uIN = mf/(rhoH*(pi/4)*(Dpipe^2)); % Inlet Pipe Velocity (m/s)
x(1:N,1) = (H/(2*N)):(H/N):(H-H/(2*N)); % Node Locations
tEND = (tTotal/delt)+1; % Last Time Step Value
t(1:tEND,1) = 0:delt:tTotal; % Simulation time at each time step (s)
% Initial Jet Penetration Depth and Plume Velocity (from Nizami [12])
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------
% NIZAMI CONSTANTS %
dpipe = Dpipe*1000; % Inlet Pipe Diameter (mm)
ahj = -0.0150*(dpipe^2) + 1.40*(dpipe) + 0.510; % Penetration Depth Parameter a
bhj = 0.00535*(dpipe) + 0.448; % Penetration Depth Parameter b
% MODELING CONSTANTS %
Rio = (g*beta*Dpipe*(Th-Tc))/(uIN^2); % Initial Richardson Number
hIo = (ahj/1000)*(Rio^(-bhj)); % Initial Penetration Depth
% INLET PLUME VELOCITY (NIZAMI METHOD) %
mfp = mf*1.062*(Rio^(-0.278)); % Inlet Plume Mass Flow Rate (kg/s)
up = mfp/(rho*(pi/4)*(Dtank^2)); % Inlet Plume Velocity (m/s)
% Time Interval 1 (0<t<to)
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------
% Determine TI and hj for 0<t<to:
% ----------------------------------------------------------------------
A2 = (pi/4)*((Dtank^2) - ((3*Dpipe)^2)); % Cross-sectional Area of Region 2 (m2)
V1 = (pi/4)*((3*Dpipe)^2)*hIo; % Volume of Region 1 (m3)
Ci = 0.00001; % Initial Thermocline Location (m)
Tau1 = -rho*V1/mfp; % Time Constant for Region 1
CT2 = Tc*rho*A2*Ci + mf*(Th-Tc)*Tau1; % Constant Term for Region 2
for i = 1:tEND
T2Top(i,1) = mfp*Tc*t(i,1) + mf*(Th-Tc)*(t(i,1)-Tau1*exp(t(i,1)/Tau1))+CT2;
T2Bot(i,1) = mfp*t(i,1) + rho*A2*Ci;
T2(i,1) = T2Top(i,1)/T2Bot(i,1); % Temperature of Region 2
Co(i,1) = (mfp/(A2*rho))*t(i,1) + Ci; % Location of Thermocline
if i == 1
TIo(i,1) = Tave;
else
TIo(i,1) = TIo(i-1,100);
end
% Iterate to solve for the jet penetration depth, h, and inlet
% mixing region temperature, TIo:
for j = 1:100
h(i,j)=(ahj/1000)*((g*beta*(Th-TIo(i,j))*Dpipe)/(uIN^2))^(-bhj);
TIo(i,j+1) = (Co(i,1)*T2(i,1) + (h(i,j) - Co(i,1))*Tc)/(h(i,j));
end
Ri(i,1) = (g*beta*(Th-TIo(i,100))*Dpipe)/(uIN^2);
hj(i,1) = h(i,100);
end
% Calculate end of time interval 1, to
% ----------------------------------------------------------------------
for i = 2:tEND
if hj(i,1) <= Co(i,1) % Determine when thermocline passes jet penetration depth
n = i-1;
break
end
end
% Determine time at end of time interval 1, to:
to = t(n,1) + ((hj(n,1) - Co(n,1))*(t(n+1,1) - t(n,1)))/((Co(n+1,1) - Co(n,1)) - (hj(n+1,1) - hj(n,1)));
TIto = TIo(n,100); % Temperature at the end of time interval 1 (degC)
hto = hj(n,1); % Jet Penetration Depth at end of time interval 1 (m)
% Time Interval 2 (t>to)
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------
% Determine inlet mixing region temperature, TI, for t>to
% ----------------------------------------------------------------------
CTIn = TIto*hto - (mf*Th*to)/(rho*Ac);
for i = 1:tEND
TIn(i,1) = ((mf*Th*t(i,1))/(rho*Ac) + CTIn)/(u*(t(i,1) - to) + hto);
end
% Determine Overall TI
% ----------------------------------------------------------------------
for i = 1:tEND
if i <= n
TI(i,1) = T2(i,1);
else
TI(i,1) = TIn(i,1);
end
end
% Thermocline Location for t>to
% ----------------------------------------------------------------------
for i = 1:tEND
if t(i,1) <= to
C(i,1) = up*t(i,1) + Ci;
else
C(i,1) = C(i-1,1) + u*delt;
end
end
% Final Temperature Profile Calculation - Thermocline Thickness Method
% ----------------------------------------------------------------------
% ----------------------------------------------------------------------
TaveS = 0.5*(TI + Tc);
time(2,1) = delt;
Tth(1,1:N) = Tc;
for i = 2:tEND-1
for j = 1:N
if x(j,1) > C(i,1) % Below Thermocline
Tth(i,j) = TaveS(i,1) + (Tc-TaveS(i,1))*erf((x(j,1)-C(i,1))/sqrt(4*alfa*time(i,1)));
elseif x(j,1) < C(i,1) % Above Thermocline
Tth(i,j) = (Tc + TI(i,1)) - (TaveS(i,1) +(Tc-TaveS(i,1))*erf((C(i,1)-x(j,1))/sqrt(4*alfa*time(i,1))));
end
end
% Determine Thermocline Thickness at the end of the time step
% ------------------------------------------------------------------
Wth(i,1) = 2*sqrt(4*alfa*time(i,1))*erfinv(1 + 0.001*(Tc/(Tc - TaveS(i,1))));
% Determine Time for profile with temperatures of next time step to
% obtain the same thermocline thickness
% ------------------------------------------------------------------
tnew(i,1) = (1/(4*alfa))*(Wth(i,1)/(2*erfinv(1 + 0.001*(Tc/(Tc-TaveS(i+1,1))))))^2;
% Calculate next 'fictitious' simulation time value
% ------------------------------------------------------------------
time(i+1,1) = tnew(i,1) + delt;
end
plot(x(1:N),Tth(3600,1:N))
title('temperature profile along the tank: Inlet Jet Mixing')
xlabel('Tank Depth(m)')
ylabel('Temperature(°C)')
0 个评论
回答(1 个)
Darren Wethington
2020-7-24
I'd put a brakpoint at line 173 or run it with the option "Pause on Errors" selected. What are the values of Tc and TaveS(i+1,1) when the error occurs? Are they real? Are they equal? (if they're equal, this would cause the value 1/(Tc-TaveS(i+1,1)) = Nan)
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!