I'm going to ask my question, explain my problem, show how I thought to solve the problem and where I got stuck, and then re-state my question.

Question: how might I solve this equation:

$$\tag{1} Re_{p}(S)=\frac{S}{1+0.15Re_{p}^{0.687}}$$ for $Re_{p}(S)$?

Explanation of my problem: I am trying to follow a paper, it goes like this:

The settling velocity for a particle falling in a viscous fluid under gravity is $$\tag{2} v_g=\sqrt{\frac{4}{3} \frac{dg}{C_D} \left(\frac{\rho_p-\rho}{\rho}\right)}$$

where $C_D$ is the drag coefficient, a function of the particle Reynolds number

$$\tag{3} Re_p=\frac{v_g d}{\nu}$$

The correlation for $C_D$ is given as

$$\tag{4} C_D=\frac{24}{Re_p}\left(1+0.15Re_{p}^{0.687}\right)$$

The three equations (2)-(4) can be reduced to a single nonlinear equation with respect to the Reynolds number

$$\tag{5} Re_{p}(S)=\frac{S}{1+0.15Re_{p}^{0.687}}$$


$$\tag{6} S=\frac{d^3 g}{18 \nu^2}\left(\frac{\rho_p-\rho}{\rho}\right)$$

The authors of the paper then say the solution of Eqn (5) is shown in the Figure below (as the solid line), which allows one to find the settling velocity

$$\tag{7} v_g=\frac{\nu Re_p(S)}{d}$$

as a function of the settling number $S$.

enter image description here

How I thought to solve the problem and where I got stuck: I thought I might be able to solve for $Re_p(S)$ by rearranging Eqn (1) as so

$$\tag{8} S=Re_p(1+0.15Re_p^{0.687})$$



and then procede by performing the "completing the square" method, e.g., dividing through by $0.15$:

$$\tag{10} Re_p^{1.687}+\frac{Re_p}{0.15}=\frac{S}{0.15}$$

and then completing the square by adding $\left(\frac{Re_p}{0.3}\right)^2$ to both sides:

$$\tag{11} Re_p^{1.687}+\frac{Re_p}{0.15}+\left(\frac{Re_p}{0.3}\right)^2=\frac{S}{0.15}+\left(\frac{Re_p}{0.3}\right)^2$$

However, since my first term on the LHS is not to the power of 2 I can't proceed to simplify the LHS.

Back to the question: How did the authors of the paper solve Eqn(1) to produce the solid line in the figure given above? (Solved for $Re_p(S)$)

(Also, please let me know if my question title is appropriate for the question I am asking)

  • $\begingroup$ Are you sure they didn't solve it numerically? $\endgroup$ – stochastic Jan 5 '18 at 21:36
  • 1
    $\begingroup$ There is no "completing the square" or related solution since the exponent is a numeric value. The solution will be done numerically. $\endgroup$ – Andreas Jan 5 '18 at 21:38
  • $\begingroup$ @Andreas what sort of numerical method should I look up to begin to understand how I might solve Eqn (1)? $\endgroup$ – Armadillo Jan 5 '18 at 22:02

This is similar to solving equations like $\tan(x) = \log(x)$. There is no way to "isolate" for $x$. In order to solve for the Reynold's numer in your equation, an iterative (i.e. numerical) approach must be taken. I assume you are in engineering, so eventually you will take a Numerical Methods course where you will learn all about solving non-linear equations.

For now, you can take the following simple iterative approach. Recall your equation $$ Re_p = \frac{S}{1 + 0.15Re_p^{0.687}}$$ Suppose you want to solve this for $S = 10^0 = 1$. Then, $$ Re_p = \frac{1}{1 + 0.15Re_p^{0.687}}$$ Make an initial guess for $Re_p$, call it $Re_p^0 = 10$. Use the above equation to iteratively get better guesses for $Re_p$. That is, $$ Re_p^{i + 1} = \frac{1}{1 + 0.15{Re_p^{i}}^{0.687}}$$ The following table summarizes the results for the first few iterations. $$\begin{array}{c | c} i & Re_p^i \\ \hline 0 & 10\\ 1 & 0.5772098\\ 2 & 0.9068355\\ 3 & 0.8770160\\ 4 & 0.8794782\\ 5 & 0.8792734\\ 6 & 0.8792904\\ 7 & 0.8792890\\ 8 & 0.8792891\\ 9 & 0.8792891\\ \end{array}$$ Of course, you'd want to automate this by writing a simple computer program.

As a final note, I suspect that the authors might have gotten values of $S$ by simply substituting values for $Re_p$ in the above equation, and plotting the result for a range of $Re_p$ that are of interest.


$\color{Red}{\textrm{Mathematica code:}}$

R[S_] := Rep /. First@NSolve[S == Rep + 15/100*Rep^(1687/1000), Rep, Reals];

r1[S_] := S; (*  a dummy function *)

LogLogPlot[{R[S], r1[S]}, {S, 10^-1, 10^3}, 
PlotRange -> {Automatic, {10^-1, 10^3/2}}, AxesLabel -> {"Rep", "S"}, 
PlotLabel -> "Settlement of spherical particle in viscous fluid", 
PlotStyle -> {Blue, {Red, Dashed}}, PlotLegends -> {Automatic}]

enter image description here $\color{Red}{\textrm{Maple code:}}$

R := proc (S) fsolve(S = Rep+0.15*Rep^1.687, Rep, maxsols = 1) end proc:
plots:-loglogplot(['R(S)', S], S = .1 .. 10^3, view = [0.1 .. 1000, 0.1 .. 500]);

$\color{Red}{\textrm{Maxima code:}}$

R(S):= find_root (S=Rep + 0.15*Rep^1.678, Rep, 0.10,999)$;
plot2d ([R(S),S], [S, 0.11, 999],logx,logy)$

As already said, you need some numerical method to solve the nonlinear equation $$S=R+k\,R^{\,a}\qquad \text{with}\qquad k=0.15\qquad \text{and}\qquad a=1.687$$ and, as usual, if you want to save iterations, you need a "good" guess for starting.

Playing with inequalities, you know that the solution is somewhere between $$R_1=\frac S {1+k}\qquad \text{and}\qquad R_2=\left(\frac S {1+k}\right)^{1/a}$$

So, let us take as estimate the geometric mean of these bounds $$R_0=\sqrt{R_1 R_2}=\left(\frac S {1+k}\right)^{\frac{1+a}{2a}}$$

Now, consider that you look for the zero of the function $$f(R)=R+k\,R^{\,a}-S$$ $$f'(R)=1+a\, k\, R^{\,a-1}$$ Newton iterates will be given by $$R_{n+1}=R_n-\frac{R_n+k\, R_n^{\,a}-S}{1+a\, k\, R_n^{\,a-1}}$$

For illustration purposes, let us try with $S=10^3$; Newton method will provide the following iterates $$\left( \begin{array}{cc} n & x_n \\ 0 & 219.186 \\ 1 & 170.109 \\ 2 & 165.981 \\ 3 & 165.950 \end{array} \right)$$ which is the solution for six significant figures.

For sure, this was for computing a single value. If you want to build the curve, for a second point, use the previous result as estimate. Suppose that we want now $S=1100$; then, the iterates will be $$\left( \begin{array}{cc} n & x_n \\ 0 & 165.950 \\ 1 & 176.500 \\ 2 & 176.303 \end{array} \right)$$

For example, using $S=1$, Newton iterates would be $$\left( \begin{array}{cc} n & x_n \\ 0 & 0.894667 \\ 1 & 0.879285 \\ 2 & 0.879268 \end{array} \right)$$


You could do simpler (not to say much better) if you try to inverse the model. Generating $S$ from given $R$ and using least-square method to minimize the sum of squares of relative errors, with a very simplistic model, you could end with $$R\approx 1.35463\, S^{0.869369}-0.384677\, S$$ which would give pretty good guesses to be used for $R_0$.

For the worked examples, the values generated as estimates are $164.772$, $173.771$ and $0.969949$. Newton method will probably require a couple of iterations for convergence to high accuracy.

  • $\begingroup$ This answer is incomplete as it lacks a proof of convergence of the Newton iteration, which can be easily remedied. $\endgroup$ – Hans Jan 7 '18 at 4:35
  • $\begingroup$ @Hans. The first derivative is always positive; so only one root. $\endgroup$ – Claude Leibovici Jan 7 '18 at 4:46
  • $\begingroup$ The uniqueness of the root of the original equation does not guarantee the convergence of the Newton iteration. $\endgroup$ – Hans Jan 7 '18 at 5:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.