fzero solver within a for loop stops working at a certain step

Hello!

This is my first time posting a question. I am usually good at scouring the web and finding helpful threads for my MATLAB issues, but this time I can't find something that resemble what I am facing.

To give you some context first, I have a blowdown system simulation, and at the moment I have an explicit procedure to compute the state of the gas flowing through my system in a certain time interval. I use a for loop, where at each step I update the state of the gas in the tank and then go through a series of computations to obtain the state (pressure, temperature, density, Mach number) of the gas in strategic points along the piping system (nodes).

My problem shows up when I try to calculate the Mach number at a certain node. The expression I am using is the following:

M7(k) = fzero(@(x2) m_dot(k)-A_7*x2*p7_stag(k)*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d*x2^2)^c, 0.1);

I basically try to solve the mass flow rate expression (m_dot(k), which changes at each step) for the Mach number M7.

A_7, gamma_he, R_he, d and c are known constants, while m_dot(k), p7_stag(k), T_stag(k) are values that change at every step k of the for loop but are all computed before the incriminated line.

The weird thing is that the solver works perfectly up until a certain value of k, where it suddenly gives me the following error:

Error using fzero (line 308). Initial function value must be finite and real.

I'll add a plot of M7 just to show you what happens:

The peak you see at x=65 is correct since I have an increase in the mass flow rate. Afterwards, though, it should continue to decrease gradually towards zero, instead it drops down at x=195.

Any idea what may be the cause of it and how I can solve it?

Thanks in advance!

[Davide Masiello]

@Valentina Lo Gatto It would be useful if you could share the entire code.

[Valentina Lo Gatto]

close all
syms x1 x2 x3 u x y z p q r w

%% Helium properties
R_he = 2077.1; % [J/kg K], Helium gas constant
he_mol_weight = 4.003*10^(-3); % [kg/mol], helium molecular weight
gamma_he = 1.6667; % helium's specific heats ratio
Cp = 5190; % [J/(kg K)] helium's specific heat at constant pressure
Cv = Cp - R_he; % helium's specific heat at constant volume
alpha_he = 1.36174*10^(-6); % [K/Pa], helium coefficient for real gas law 

%% Tank properties
D_tank = 0.0254*10.6; % 0.269; % [m], inner diameter of the tank
L_tank = 0.0254*42.2; % 1.072; % [m], length of the tank
A_tank = 4*(pi*D_tank^2 + D_tank*L_tank); % [m^2], surface area of 4 tanks
V_tank = 4*((1/39.37)^3)*2756; % 0.04516; % [m^3], total volume for 4 tanks
p1_0 = 3844*6895; % [Pa], tank initial pressure
T1_0 = 14.6 + 273.1; % [K], tank initial temperature

%% rho calculations from real gas law:
z_real_gas_law_0 = 1+alpha_he*(p1_0/T1_0);
rho1_0 = p1_0/(z_real_gas_law_0*R_he*T1_0); % [kg/m^3], initial density inside the tank real gas

m_tank_0 = rho1_0*V_tank; % [kg], mass of helium in the tank before the start of inflation

%% Pipes dimensions
D2 = (0.250 - 2*0.035)*0.0254; 
D5 = (0.500 - 2*0.049)*0.0254; 

L_23 = 1; % [m], length of the pipe from node 2 to node 3
L_56 = 4; % [m], length of the pipe from node 5 to node 6

A_23 = (pi/4)*D2^2; % [m^2], cross-sectional area of the pipe from node 2 to node 3
A_56 = (pi/4)*D5^2; % [m^2], cross-sectional area of the pipe from node 5 to node 6

k_minor_23 = 1*0.1 + 1*0.3 + 4* 0.5; % minor losses due to tube bending, opening of valves, tee

%% Orifice plates dimensions and times of opening
Cd = 0.84; % discharge coefficient

d1 = 0.08*0.0254; % [m], first orifice plate diameter.  
d2 = 0.1*0.0254; % [m], second orifice plate diameter.
d3 = 0.0*0.0254; % [m], third orifice plate diameter. We are not using it in this instance.

t0 = 1; % [s], time before the first orifice is fully open
t1 = 65; % [s], between t0 and t1 only the first orifice valve is open
t2 = 66; % [s], between t1 and t2 the first orifice is fully open and the second is opening
t3 = 246; % [s], between t2 and t3 both the first and second orifices are fully open
t4 = 247; % [s], between t3 and t4 the first and second orifices are fully open while the third orifice is opening. After t4 all three are fully open.

t_max_infl = 600; % [s], max time for the inflation tests (10 minutes)

%% Diffuser plates
d_diff = 0.04; % [m], nominal diffuser diameter
A_7 = (pi/4)*d_diff^2; % [m^2], area of the diffuser chamber
n_diff_holes = 16; % number of holes on the diffuser's exit surface
d_diff_holes = 0.0254*0.25; % 6/1000; % [m], holes diameter --> 6 mm
A_8 = n_diff_holes*(pi/4)*d_diff_holes^2; % [m^2], diffuser's exit holes area

%% Computation
p10 = 67000; % [Pa], pressure at node 10

a = -1/(gamma_he-1);
b = gamma_he/(gamma_he-1);
c = -0.5*(gamma_he+1)/(gamma_he-1); 
d = (gamma_he-1)/2;

delta_t = 0.2;
l_array = t_max_infl/delta_t;

t = zeros (1,l_array);
m_dot = zeros (1,l_array);
T_dot = zeros (1,l_array);
m_tank = zeros(1,l_array);
m_inj = zeros(1,l_array);
m_perc = zeros(1,l_array);
T_tank = zeros(1,l_array);
k_he = zeros(1,l_array);
n_he = zeros(1,l_array);
p_tank = zeros(1,l_array);
rho_tank = zeros(1,l_array);
A_4 = zeros(1,l_array);
T_stag = zeros(1,l_array);
M2 = zeros(1,l_array);
p2_stag = zeros(1,l_array);
p2 = zeros(1,l_array);
T2 = zeros(1,l_array);
rho_2 = zeros(1,l_array);
Re_23 = zeros(1,l_array);
M3 = zeros(1,l_array);
p3_stag = zeros(1,l_array);
p3 = zeros(1,l_array);
T3 = zeros(1,l_array);
rho_3= zeros(1,l_array);
M4 = zeros(1,l_array);
p4_stag = zeros(1,l_array);
M5 = zeros(1,l_array);
p5_stag = zeros(1,l_array);
p5 = zeros(1,l_array);
T5 = zeros(1,l_array);
rho_5 = zeros(1,l_array);
M6 = zeros(1,l_array);
p6_stag = zeros(1,l_array);
p6 = zeros(1,l_array);
T6 = zeros(1,l_array);
rho_6= zeros(1,l_array);
M7 = zeros(1,l_array);
p_threshold = zeros(1,l_array);
p7_stag = zeros(1,l_array);
p7 = zeros(1,l_array);
T7 = zeros(1,l_array);
rho_7= zeros(1,l_array);
M8 = zeros(1,l_array);
p8_stag = zeros(1,l_array);
p8 = zeros(1,l_array);
T8 = zeros(1,l_array);
rho_8 = zeros(1,l_array);
f_factor56 = zeros(1,l_array);
L_star = zeros(1,l_array);
grad_T = zeros(1,l_array);
Q_dot = zeros(1,l_array);

for k = 1:(l_array+1)

    if k == 1 
        t(k) = 0;
        m_dot(k) = 0;
        T_dot(k) = 0;
        m_tank(k) = m_tank_0;
        T_tank(k) = T1_0;

    else 
        t(k)=t(k-1)+delta_t;
        m_tank(k) = m_tank(k-1) - m_dot(k-1)*delta_t;
        T_tank(k) = T_tank(k-1) + T_dot(k-1)*delta_t;
    end

    m_inj(k) = m_tank_0 - m_tank(k);
    m_perc(k) = m_inj(k)/m_tank_0;

    rho_tank(k)= m_tank(k)/V_tank;

   
    % pressure in the tank calculated using the real gas state equation
    p_tank(k) = (rho_tank(k)*R_he*T_tank(k))/(1-alpha_he*rho_tank(k)*R_he);
    
    % Stagnation Temperature
    T_stag(k) = T_tank(k);

    % Node 2, stagnation pressure
    p2_stag(k) = p_tank(k);

    % Node 4, orifices area calculation, A_4
    if (t(k)<t0)
        A_4(k) = 0.1*(pi/4)*d1^2 + 0.9*(pi/4)*d1^2*(t(k)/t0);
    elseif (t(k)<t1)
        A_4(k) = (pi/4)*d1^2;
    elseif (t(k)<t2)
        A_4(k) = (pi/4)*d1^2 + (pi/4)*d2^2*(t(k)-t1);
    elseif (t(k)<t3)
        A_4(k) = (pi/4)*d1^2 + (pi/4)*d2^2;
    elseif (t(k)<t4)
        A_4(k) = (pi/4)*d1^2 + (pi/4)*d2^2 + (pi/4)*d3^2*(t(k)-t3);
    else
        A_4(k) = (pi/4)*d1^2 + (pi/4)*d2^2 + (pi/4)*d3^2;
    end
    
    % NOTE: M4 = 1 --> choked condition at the orifice(s)

    % I assume, then, that m_dot_3=m_dot_4 plus the condition p3_stag=p4_stag so I can compute M3

    f1 = @(x1) x1*(1+d*x1^2)^c-(Cd*A_4(k)*(1+d)^c)/A_23;
    M3(k) = fzero (f1, 1);

    T3(k) = T_stag(k)/(1 + d*M3(k)^2);

    % Next step is to compute m_dot, p3_stag and M2 from a set of 3 equations.
    eq1 = x-A_23*y*p2_stag(k)*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d*y^2)^c == 0; % m_dot_2 expression
    eq2 = x-Cd*A_4(k)*z*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d)^c == 0; % m_dot_4 expression
    eq3 = (p2_stag(k)*(1+d*y^2)^(-b)-z*(1+d*M3(k)^2)^(-b))*(p2_stag(k)*(1+d*y^2)^a+z*(1+d*M3(k)^2)^a)-gamma_he*(((gamma_he+1)/2)^(2*c))*((Cd*A_4(k)/A_23)^2)*(z^2)*(k_minor_23+(L_23/D2)*(1.8*log10((4/(pi*D2*6.9*(6.731*10^-6 + (4.569*10^(-14))*p2_stag(k)*(1+d*y^2)^(-b) + (4.330*10^-8)*T_stag(k)/(1+d*y^2))))*(Cd*A_4(k)*z*sqrt(gamma_he/(R_he*T_stag(k)))*((0.5*(gamma_he+1))^c))))^(-2)) == 0;

    [m_dot(k),M2(k),p3_stag(k)] = vpasolve([eq1, eq2, eq3], [x, y, z], [0 Inf; -1 1; 0 Inf]);
    
    % ) Node 2, other conditions
    p2(k) = p2_stag(k)*(1 + d*M2(k)^2)^(-b);
    T2(k) = T_stag(k)/(1 + d*M2(k)^2);
    rho_2(k) = p2(k)/(R_he*T2(k));
    % ) Node 3 other conditions
    p3(k) = p3_stag(k)*(1 + d*M3(k)^2)^(-b);
    rho_3(k) = p3(k)/(R_he*T3(k));

    % Next step is to compute p5_stag, M5 and p6_stag from another set of 3 equations. !!!!! NOTE: we are assuming M6=1 
    eq4 = A_56*p*q*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d*q^2)^c-m_dot(k) == 0; % m_dot_5 expression
    eq5 = A_56*r*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d)^c-m_dot(k) == 0;
    % eq6 = (1.8*log10((4*m_dot(k))/(pi*D5*(6.731*10^-6+(4.569*10^(-14))*p*(1+d*q^2)^(-b)+(4.330*10^-8)*T_stag(k)/(1+d*q^2))*6.9)))^(-2)*(L_56/D5)*(R_he*T_stag(k))*(m_dot(k)/A_56)^2-(p*(1+d*q^2)^(-b)-r*(1+d)^(-b))*(p*(1+d*q^2)^a+r*(1+d)^a) == 0;
    eq6 = ((1.8*log10((4*m_dot(k))/(pi*D5*(6.731*10^-6+(4.569*10^(-14))*p*(1+d*q^2)^(-b)+(4.330*10^-8)*T_stag(k)/(1+d*q^2))*6.9)))^(-2)*L_56)/D5-(1-q^2)/(gamma_he*q^2)+((gamma_he+1)/(2*gamma_he))*log(((gamma_he+1)*q^2)/(2+(gamma_he-1)*q^2)) == 0;
    
    [p5_stag(k),M5(k),p6_stag(k)] = vpasolve([eq4, eq5, eq6], [p, q, r], [0 Inf; 0 0.5; 0 Inf]);
    
    % ) Node 5 other conditions
    p5(k) = p5_stag(k)*(1 + d*M5(k)^2)^(-b);
    T5(k) = T_stag(k)/(1 + d*M5(k)^2);
    rho_5(k) = p5(k)/(R_he*T5(k));

    % ) Node 6 other conditions (Remember we are assuming M6=1)
    p6(k) = p6_stag(k)*(1 + d)^(-b);
    T6(k) = T_stag(k)/(1 + d);
    rho_6(k) = p6(k)/(R_he*T6(k));

    % ) Node 7 and 8. Note that we assume no loss happens between the diffuser body and the exit plate, so: p7_stag(k) = p8_stag(k)
    
    % ) Node 8
    p_threshold(k) = m_dot(k)/(A_8*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d)^c);

    if p_threshold(k)>2.049*p10
        M8(k) = 1;
        p8_stag(k) = p_threshold(k);
        p8(k) = p8_stag(k)*(1+d)^(-b);        
    else
        p8(k) = p10;

        f4 = @(x3) p8(k)*(1+d*x3^2)^b-m_dot(k)/(A_8*x3*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d*x3^2)^c); 
        M8(k) = fzero (f4, 0.8);

        p8_stag(k) = p8(k)*(1+d*M8(k)^2)^b;      
    end

    T8(k) = T_stag(k)/(1 + d*M8(k)^2);
    rho_8(k) = p8(k)/(R_he*T8(k));
    
    % ) Node 7
    p7_stag(k) = p8_stag(k);

    M7(k) = fzero(@(x2) m_dot(k)-A_7*x2*p7_stag(k)*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d*x2^2)^c, 0.1);
   
    p7(k) = p7_stag(k)*(1 + d*M7(k)^2)^(-b);
    T7(k) = T_stag(k)/(1 + d*M7(k)^2);
    rho_7(k) = p7(k)/(R_he*T7(k));

    % Other updates before next step in the loop:
    k_he(k) = 0.04876 + p_tank(k)*6.318*10^(-14) + T_tank(k)*4.330*10^(-8);
    Q_dot(k) = k_he(k)*A_tank*(T_0_atm - T_tank(k))/(D_tank/4); % [W/m^2] conductive heat flow over time from environment to tank (differential form of Fourier's law of thermal conduction) 

    T_dot(k) = -(m_dot(k)*(gamma_he-1)*T_tank(k))/m_tank(k) + Q_dot(k)/(m_tank(k)*Cv);

end

It is quite a long code, so I extracted the relevant sections. Also apologies for the (likely many) inefficiencies you may find along the way.

I should mention that the code works fine if I comment out the equations related to Node 7 calculations (to be fair, I do have a second issue with the solution of the system (eq4, eq5, eq6) when I use the commented eq.6 instead of the uncommented one, but it is better not to mix the topics).

[Riccardo Scorretti]

When I try to run your script, I get an error at line 130:

T_tank(k) = T_tank_0;

Indeed the variable T_tank_0 is undefined at this point. Please fixe the code, so that it can be executed.

[Valentina Lo Gatto]

Sorry, you are right. As I said, I have different sections in this code and T_tank_0 was defined somewhere else but it is the same as T1_0 (defined in the code above). I have fixed it now.

T_tank(k) = T1_0;

Answer 1

[Matt J]

I am not sure why you are using fzero. It looks like the root can be solved explicitly for x2.

In any case, you can more easily see what is happening if you evaluate the function

fun = @(x2) m_dot(k)-A_7*x2*p7_stag(k)*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d*x2^2)^c

at the offending k and initial x2 and see what value it returns. If d is negative and c is non-integer, it seems likely that (1+d*x2^2)^c will be complex valued for sufficiently large x2. The sqrt could also be returning complex values, depending on the specific values of the arguments.

[Valentina Lo Gatto]

A_7, gamma_he, R_he, d and c are known constants, specifically

c = -0.5*(gamma_he+1)/(gamma_he-1) = -1.9999
d = (gamma_he-1)/2 = 0.333

At the point of failure (step k = 976 in the loop), T_stag = 77.6274 so

sqrt(gamma_he/(R_he*T_stag(k))) = 0.0032

However, upon your input, I went to check each component and it seems that p7_stag is a NaN value and the failure actually happens earlier on in the calculation of M8:

f4 = @(x3) p8(k)*(1+d*x3^2)^b-m_dot(k)/(A_8*x3*sqrt(gamma_he/(R_he*T_stag(k)))*(1+d*x3^2)^c); 
M8(k) = fzero (f4, 0.8);

I'll look into this to see if I can sort whatever is wrong with it.

Thanks for your contribution!