Consider two parallel Cylinders with Diameters of $R_1$ and $R_2$:

Contact between Two Cylinders with Parallel Axes

The contact width can be calculated from 2D Hertz formula:

$$ a=2 \sqrt{\frac{PR}{\pi E_c}} \tag{1}$$

Where $P=\frac{F}{L}$ is the force per unit length, and


and $$\frac{1}{E_c}=\frac{1-\nu_1^2}{E_1}+\frac{1-\nu_2^2}{E_2}\tag{3}$$

From here I know the total displacement can be written as:

$$ \delta=\frac{P}{\pi E_c}\left( \ln\left(\frac{8R_1}{a} \right)+\ln\left(\frac{8R_2}{a} \right)-1 \right) \tag{4}$$

I need to calculate $P$ as a function of $\delta$ so :

R := (R1*R2)/(R1 + R2)
a := 2*Sqrt[P*R/(pi*Ec)]
Solve[delta == P*(Log[8*R1/a] + Log[8*R2/a] - 1)/(pi*Ec), P]

Which gives me "a" solution and a warning:

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

enter image description here

To have a general understanding of the solution I plotted a simplified version:

Plot[-x/ProductLog[-x], {x, 0, 1/E}]

enter image description here

Which does not make sense because at $\delta=0$ the force must be zero, and it must increase by displacement increasing.

Trying to solve the equation with Reduce, as sugested in the waring:

Reduce[delta == P*(Log[8*R1/a] + Log[8*R2/a] - 1)/(pi*Ec), P]

also does not yield any results after a long time of calculation. I would appreciate if you could help me know if I'm making any mistakes and/or how to calculate $P$ versus $\delta$ for parallel cylinders.

  • $\begingroup$ If you add Reals as domain specification into Reduce, you are quickly able to see that Mathematica admits it cannot solve for P using Reduce. $\endgroup$ – Subho Apr 4 '18 at 15:49
  • $\begingroup$ I actually did that. That's why I posted the question. $\endgroup$ – Foad Apr 4 '18 at 16:09
  • $\begingroup$ There is typo in pi. E is in-build command so you have to replace with another symbole $\endgroup$ – Gopal Verma Apr 4 '18 at 16:17
  • $\begingroup$ @GopalVerma ah! silly mistake. I'm gonna fix that and come back here to report. $\endgroup$ – Foad Apr 4 '18 at 16:18
  • 2
    $\begingroup$ Pi should be capitalized. $\endgroup$ – george2079 Apr 4 '18 at 17:45

this takes a bit of extra work.. first factor the lead P into a new symbol ( deltabyp ) and solve:

R = (R1*R2)/(R1 + R2);
a = 2*Sqrt[P*R/(Pi*Ec)];
sol = Solve[deltabyp == (Log[8*R1/a] + Log[8*R2/a] - 1)/(Pi*Ec), 

{P -> 16 E^(-1 - deltabyp Ec Pi) Ec Pi (R1 + R2)}

now sub back in for deltabyp and Reduce yields a pretty reasonable form:

 Simplify[Reduce[(P /. sol /. deltabyp -> delta/P  ) == P, P, Reals], 
   Assumptions -> {R1 > 0, R2 > 0, Ec > 0, delta > 0}]

enter image description here

plot the two solutions:

With[{R1 = 1, R2 = 1, Ec = 1}, 
  -((delta Ec \[Pi])/ProductLog[-((delta E)/(16 (R1 + R2)))]),
  -((delta Ec \[Pi])/ProductLog[-1, -((delta E)/(16 (R1 + R2)))])},
  {delta, 0, 16 (R1 + R2)/E^2}]]

enter image description here

so you see the desired result is the second one,

enter image description here

note that ProductLog is mathematica's name for the Lambert-W function.


$$\implies P=-\frac{\pi \delta E_c}{W_{-1}\left(- \frac{e \delta}{16\left( R_1 +R_2 \right)} \right)}$$

They curiously avoid mentioning that in the docs for ProductLog , however there is an undocumented LambertW that simplifies to ProductLog

  • 1
    $\begingroup$ solve returns a list of solutions, in this case there is only one solution, so the [[1]] effectively gets rid of an extra set of { }. $\endgroup$ – george2079 Apr 4 '18 at 21:06
  • 1
    $\begingroup$ mathworld.wolfram.com/LambertW-Function.html $\endgroup$ – george2079 Apr 5 '18 at 11:50
  • 1
    $\begingroup$ my understanding ProductLog[-1, z] is LambertW[-1, z] which is W subscript -1 on the wiki page. $\endgroup$ – george2079 Apr 5 '18 at 14:07
  • 1
    $\begingroup$ @Foad You can think of $W$ like you do $\sqrt{\cdot}$. Both give you a solution to a certain equation; $W(x)$ satisfies $ x = W(x) \exp{W(x)}$, and $\sqrt{x}$ satisfies $x = (\sqrt{x})^2$. However, $-\sqrt{x}$ also satisfies the second equation. It turns out that, similarly, there are multiple solutions to $y=xe^x$, so if you want $W$ to return an unambiguous number (i.e., be a bona fide function), you need to specify which of the multiple solutions you would like. $k$ is an index that lets you do this. $\endgroup$ – sourisse Apr 5 '18 at 15:36
  • 1
    $\begingroup$ @Foad It doesn't "mean" anything, it's just a convention. You could just as well have called it the $\star$th solution, as long as everybody agreed on which one that was. You can blame Knuth et. al for popularizing this notation: link.springer.com/article/10.1007%2FBF02124750 $\endgroup$ – sourisse Apr 5 '18 at 15:47

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.