How to resolve error using erfinv: Input must be real and full?

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)')

Answer 1

[Darren Wethington]

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)