Four equations and four unknown variables in Mathematica

Question

I'm trying to solve a problem in Mathematica where I have 4 unknown variables and 4 equations, but it doesn't quite work for me. I would appreciate if someone could tell me what I'm doing wrong.

Clear["Global`*"]
vdot = 20;
Temp = 10;
ρ = QuantityMagnitude[
   ThermodynamicData["Water", 
    "Density", {"Temperature" -> Quantity[Temp, "DegreesCelsius"]}]];
μ = QuantityMagnitude[
   ThermodynamicData["Water", 
    "Viscosity", {"Temperature" -> Quantity[Temp, "DegreesCelsius"]}]];
ϵ = 2*10^-6;
g = 9.81;
Le = 200;
eq1 = Reyn == (ρ v Diam)/μ;
eq2 = 1/Sqrt[
   f] == -2.0 Log10[(ϵ/Diam)/3.7 + 2.51/(Reyn Sqrt[f] )];
eq3 = 17 == f Le/Diam v^2/(2 g);
eq4 = v == vdot/((Diam/2)^2 π);
Solve[{eq1, eq2, eq3, eq4}, {Reyn, v, f, Diam}]

Answer 1

This works

sol = FindRoot[{eq1, eq2, eq3, eq4}, {{Reyn, 10^7}, {v, 1}, {f, 1}, {Diam, 10}}]

{Reyn -> 1.55996*10^7, v -> 16.3066, f -> 0.00783754, Diam -> 1.24965}

{eq1, eq2, eq3, eq4} /. sol

{True, True, True, True}

Answer 2

You can try to act as follows. First, let us leave your parameters of the equation to stay symbols and solve the first, third and fourth equations:

eq1 = Reyn == (\[Rho] v Diam)/\[Mu]
eq2 = 1/Sqrt[f] == -2.0 Log10[(\[Epsilon]/Diam)/3.7 + 2.51/(Reyn Sqrt[f])]
eq3 = 17 == f Le/Diam v^2/(2 g)
eq4 = v == vdot/((Diam/2)^2 \[Pi])

sl = Solve[{eq1, eq3, eq4}, {Reyn, v, Diam}]

I guess you only need the solution in real values. The third of the solutions seems to be like this:

sl[[3]]

(* {Reyn -> (2 2^(2/5) 17^(1/5) g^(1/5) vdot^(3/5) ρ)/(
  f^(1/5) Le^(1/5) π^(3/5) μ), 
 v -> (2^(4/5) 17^(2/5) g^(2/5) vdot^(1/5))/(
  f^(2/5) Le^(2/5) π^(1/5)), 
 Diam -> (2^(3/5) f^(1/5) Le^(1/5) vdot^(2/5))/(
  17^(1/5) g^(1/5) π^(2/5))}  *)

Let us substitute it into eq2, and after that let us substitute your numbers instead of the symbols. For the latter aim let us first prepare a rule:

rule = {vdot -> 20, Temp -> 10, ρ ->QuantityMagnitude[
    ThermodynamicData["Water","Density", {"Temperature" ->Quantity[Temp, "DegreesCelsius"]}]], μ ->QuantityMagnitude[ThermodynamicData["Water", 
     "Viscosity", {"Temperature" ->Quantity[Temp, "DegreesCelsius"]}]],ϵ -> 2*10^-6,g -> 9.81, Le -> 200}

Now let us substitute:

eq2New = eq2 /. sl[[3]] /. rule

(*  1/Sqrt[f] == -0.868589 Log[
   4.24351*10^-7/f^(3/10) + 1.64012*10^-7/f^(1/5)]  *)

Now it becomes clear, why Mma does not return you any solution. Such transcendental equations are typically not solved with the function Solve. It can, however, be solved numerically. Let us first draw the parts of the equation:

[![Plot\[{-0.869 Log\[4.2435*^-7/f^(3/10) + 1.64`*^-7/f^(1/5)\], 1/Sqrt\[
  f\]}, {f, 0, 0.02}\]][1]][1]

we see that there is the cross-section somewhere below 0.01. Now we can safely solve it:

FindRoot[eq2New, {f, 0.05}]

(*  {f -> 0.00783754}  *)

Done. Have fun!

Answer 3

Solve is an exact solver; use exact values:

Clear["Global`*"]
vdot = 20;
Temp = 10;

ρ = QuantityMagnitude[
    ThermodynamicData["Water", 
     "Density", {"Temperature" -> Quantity[Temp, "DegreesCelsius"]}]] // 
   Rationalize[#, 0] &;

μ = QuantityMagnitude[
    ThermodynamicData["Water", 
     "Viscosity", {"Temperature" -> Quantity[Temp, "DegreesCelsius"]}]] // 
   Rationalize[#, 0] &;

ϵ = 2*10^-6;
g = 981/100;
Le = 200;
eq1 = Reyn == (ρ v Diam)/μ;
eq2 = 1/Sqrt[
     f] == -2 Log10[(10 ϵ/Diam)/37 + 251/(100 Reyn Sqrt[f])];
eq3 = 17 == f Le/Diam v^2/(2 g);
eq4 = v == vdot/((Diam/2)^2 π);

Include the constraints that all variables are positive

sol = Solve[{eq1, eq2, eq3, eq4, Reyn > 0, v > 0, f > 0, Diam > 0}, 
 {Reyn, v, f, Diam}];

Verifying the solution

{eq1, eq2, eq3, eq4, Reyn > 0, v > 0, f > 0, Diam > 0} /. sol // 
 FullSimplify

(* {{True, True, True, True, True, True, True, True}} *)

The approximate numeric values are

sol[[1]] // N

(* {Reyn -> 1.55996*10^7, v -> 16.3066, f -> 0.00783754, Diam -> 1.24965} *)