I am getting unrealistic values using vpasolve

clc, clear
d = 5/12; % Pipe diameter, ft
L = 1; % Pipe length, ft
h = 2/12; % Height difference of mercury, ft
rho_w = 1.937; % density of water, slug/ft^3
rho_m = 26.3; % density of mercury, slug/ft^3
mu = 2.09*10^-5; % slug/ft*second
g = 32.2; % gravity, ft/s^2

syms Vcl hf f V Re_d
x = vpasolve([(rho_w*Vcl^2)/2 == ((rho_m-rho_w)*g*h)+(rho_w*g*hf), ...
              hf == f*(L/d)*((V^2)/(2*g)), V/Vcl == (1+1.3*sqrt(f))^-1, ...
              1/sqrt(f) == -2*log10(2.51/(Re_d*sqrt(f))), ...
              Re_d == (rho_w*V*d)/mu], ...
                                        [Vcl, hf, f, V, Re_d])

I am getting values with an i in them, and I am unsure how to fix this. I have confirmed that the equations in the solver are correct.

This is the answer I am getting:

x =

struct with fields:

Vcl: - 11.618976698556484163757444664859 - 0.0000000000000000031159299773113211212032247085422i

hf: 0.00000000010691099774793287072143003602414 + 0.0000000000000000011243452733140998420596898170279i

f: 4491158430069173668821.2496305322 + 12026888742356199124.612735277612i

V: 0.00000000000000054448854808565276662309009580122 - 0.00000000000000000072904648955139354037146343208172i

Re_d: 0.000000000021026202504822758976257513492856 - 0.00000000000002815317085847387421157383877638i

Answer 1

[Paul]

I had to rearrange the second and third equations to muliply through by Vcl and sqrt(f) respectively and restrict the range of the solution space in the call to vpasolve()

d = 5/12; % Pipe diameter, ft
L = 1; % Pipe length, ft
h = 2/12; % Height difference of mercury, ft
rho_w = 1.937; % density of water, slug/ft^3
rho_m = 26.3; % density of mercury, slug/ft^3
mu = 2.09*10^-5; % slug/ft*second
g = 32.2; % gravity, ft/s^2
syms Vcl hf f V Re_d
x = vpasolve([(rho_w*Vcl^2)/2 == ((rho_m-rho_w)*g*h)+(rho_w*g*hf), ...
              hf              == f*(L/d)*((V^2)/(2*g)), ...
              V               == Vcl*(1+1.3*sqrt(f))^-1, ...
              1               == -2*log10(2.51/(Re_d*sqrt(f)))*sqrt(f), ...
              Re_d            == (rho_w*V*d)/mu], ...
              [Vcl, hf, f, V, Re_d], ...
              [ -inf inf; ...
               -inf inf; ...
                  0 inf; ...
                  -inf inf
                  -inf inf])

x = struct with fields:
     Vcl: 11.766003641583686200483461147558
      hf: 0.053388543152905222463767286417839
       f: 0.013742615972609076804165352881268
       V: 10.210021170256005402469178079995
    Re_d: 394274.54160258943110540850414583

Answer 2

[Walter Roberson]

As a third parameter to vpasolve pass in a 5 x 2 array of constraints, lower bounds and upper bounds. In particular f should be used constrained to lower bound 0.