Hertzian contact mechanics: symbolically calculating the required force for a specific displacement

Consider two parallel Cylinders with Diameters of $R_1$ and $R_2$: 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

$$\frac{1}{R}=\frac{1}{R_1}+\frac{1}{R_2}\tag{2}$$

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. To have a general understanding of the solution I plotted a simplified version:

Plot[-x/ProductLog[-x], {x, 0, 1/E}] 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.

• 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. – Subho Apr 4 '18 at 15:49
• I actually did that. That's why I posted the question. – Foad Apr 4 '18 at 16:09
• There is typo in pi. E is in-build command so you have to replace with another symbole – Gopal Verma Apr 4 '18 at 16:17
• @GopalVerma ah! silly mistake. I'm gonna fix that and come back here to report. – Foad Apr 4 '18 at 16:18
• Pi should be capitalized. – 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][]

{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}] plot the two solutions:

With[{R1 = 1, R2 = 1, Ec = 1},
Plot[{
-((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}]] so you see the desired result is the second one, 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)}$$
• @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. – sourisse Apr 5 '18 at 15:36
• @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 – sourisse Apr 5 '18 at 15:47