Solving for flow speed?

I'd like to use Mathematica to verify the results of a Matlab code that iteratively solves for a flow velocity, u, given some fluid properties and a pressure drop. Ideally, I would like solve the system explicity for u.

I've set up the equations like so:

ReNum = (rho u L) / mu;


In the expression above, rho, L, and mu are constants, u is the variable of interest. The next expression is a description for a "friction factor", f.

ffEqn = 1/Sqrt[f] == -2.0 Log10[(eps / d)/3.7 + 2.51/(ReNum Sqrt[f])]


Here, d, and eps are constants, and the factor f is influenced by the flow velocity and properties (since ReNum appears in the description). The equation below relates the pressure drop to the flow velocity, adjusted by the factor, f. Here, L, d, and rho are constants related to the system and fluid properties.

dwEqn = deltaP == f (L/d (rho u^2)/2)


I've been naively trying to apply Solve and Reduce, but without success. Solving ffEqn for f, and subbing that in to dwEqn, and applying solve results in Solve::ifun and Solve::inex warnings and Solve returns unevaluated:

Solve[deltaP == (1.52648*10^21 L u^2 rho)/(
d ((3.53972*10^19 u^2 eps^2 rho^2)/mu^2 - (
5.7107*10^20 u eps rho ProductLog[(
0.458682 d E^((0.123968 u eps rho)/mu)
u rho)/mu])/mu +
2.30329*10^21 ProductLog[(
0.458682 d E^((0.123968 u eps rho)/mu)
u rho)/mu]^2)), u]


Is there some reason why this approach shouldn't work?

• "Is there some reason why this approach shouldn't work?" - Because it's a transcendental equation? – xzczd Mar 20 '15 at 2:42

You can't solve that monster algebraically. One could create a function that takes all the parameters and numerically finds the velocity. u0is a starting point

solveU[u0_, {deltaP_, L_, rho_, d_, eps_, mu_}] :=
u /. FindRoot[deltaP == (1.52648*10^21 L u^2 rho)/(
d ((3.53972*10^19 u^2 eps^2 rho^2)/mu^2 - (
5.7107*10^20 u eps rho ProductLog[(
0.458682 d E^((0.123968 u eps rho)/mu)
u rho)/mu])/mu +
2.30329*10^21 ProductLog[(
0.458682 d E^((0.123968 u eps rho)/mu)
u rho)/mu]^2)),
{u, u0}]


For example:

solveU[1, {100, 10, 1000, 0.4, 10^-5, 8.9*10^-4}]
(*Output*)
0.744147


You can even make a Plot

Plot[
solveU[1, {400, 10, 1000, d, 10^-5, 8.9*10^-4}],
{d, 0.001, 1},
AxesLabel -> {"Diameter (m)", "Velocity (m/s)"}, GridLines -> Automatic]