Overall deflection
Considering that there is a small region of contact, and we can use the Hertzian model it seems that there is an analytical solution 1 (although I would not call this crushing)
$$2\delta = \frac{P}{L} (V_1 + V_2)\left[1 + \log\left\{\frac{2L^3}{(V_1 +V_2) P d}\right\}\right]$$
where $V_i = (1 - \nu_i^2)/(\pi E)$. If we assume that the planes are infinitely rigid compared to the cylinder we obtain
$$2\delta = \frac{P V_1}{L} \left[1 + \log\left\{\frac{2L^3}{V_1 P d}\right\}\right]$$
or
$$2 \delta = \frac{P}{\pi E_1 L} \left(1 - \nu_1^2\right) \log{\left[\frac{2 \pi E_1 L^3}{d P \left(1 - \nu_1^2\right)} \right]}$$
This equation can be inverted to obtain
$$P = \frac{2 \pi E_1 L^3}{d \left(1 -\nu_1^2\right)} e^{\operatorname{LambertW}{\left (- \frac{d \delta}{L^2} \right )}}$$
Stress at the interior
We can model the cylinder as a 2D problem: a disk with radial forces in the poles. The stress function for a disk of diameter $d$ with center in the origin, and radial inward and opposite forces $P$ placed at $(0, d/2)$ and $(0, -d/2)$ is given by
$$\phi = x\arctan\left[\frac{x}{d/2 - y}\right] + x\arctan\left[\frac{x}{d/2 + y}\right] + \frac{P}{\pi d}(x^2 + y^2)$$
We know that the stresses are given by
\begin{align}
\sigma_{xx} = \frac{\partial^2 \phi}{\partial x^2}\\
\sigma_{yy} = \frac{\partial^2 \phi}{\partial y^2}\\
\sigma_{xy} = -\frac{\partial^2 \phi}{\partial x \partial y}
\end{align}
that gives
$$\sigma_{xx} = 2 \left[\frac{P}{\pi d} - \frac{32 x^{4}}{\left(d + 2 y\right)^{5} \left(\frac{4 x^{2}}{\left(d + 2 y\right)^{2}} + 1\right)^{2}} - \frac{32 x^{4}}{\left(d - 2 y\right)^{5} \left(\frac{4 x^{2}}{\left(d - 2 y\right)^{2}} + 1\right)^{2}} + \frac{8 x^{2}}{\left(d + 2 y\right)^{3} \left(\frac{4 x^{2}}{\left(d + 2 y\right)^{2}} + 1\right)} + \frac{8 x^{2}}{\left(d - 2 y\right)^{3} \left(\frac{4 x^{2}}{\left(d - 2 y\right)^{2}} + 1\right)}\right]$$
$$\sigma_{yy} = 2 \left[\frac{P}{\pi d} - \frac{8 x^{2}}{\left(d + 2 y\right)^{3} \left(\frac{4 x^{2}}{\left(d + 2 y\right)^{2}} + 1\right)^{2}} - \frac{8 x^{2}}{\left(d - 2 y\right)^{3} \left(\frac{4 x^{2}}{\left(d - 2 y\right)^{2}} + 1\right)^{2}} + \frac{2}{\left(d + 2 y\right) \left(\frac{4 x^{2}}{\left(d + 2 y\right)^{2}} + 1\right)} + \frac{2}{\left(d - 2 y\right) \left(\frac{4 x^{2}}{\left(d - 2 y\right)^{2}} + 1\right)}\right]$$
$$\sigma_{xy} = - 8 x \left[\frac{4 x^{2}}{\left(d + 2 y\right)^{4} \left(\frac{4 x^{2}}{\left(d + 2 y\right)^{2}} + 1\right)^{2}} - \frac{4 x^{2}}{\left(d - 2 y\right)^{4} \left(\frac{4 x^{2}}{\left(d - 2 y\right)^{2}} + 1\right)^{2}} - \frac{1}{\left(d + 2 y\right)^{2} \left(\frac{4 x^{2}}{\left(d + 2 y\right)^{2}} + 1\right)} + \frac{1}{\left(d - 2 y\right)^{2} \left(\frac{4 x^{2}}{\left(d - 2 y\right)^{2}} + 1\right)}\right]$$
and for strains
\begin{align}
\epsilon_{xx} &= \frac{1}{E}(\sigma_{xx} - \nu \sigma_{yy})\\
\epsilon_{yy} &= \frac{1}{E}(\sigma_{yy} - \nu \sigma_{xx})\\
\epsilon_{xy} &= \frac{\sigma_{xy}}{G} \, .
\end{align}
For displacements, there are two options that come to my mind.
- Rewrite the stress function in polar coordinates, and then use the Mitchell solution for displacements. The stress function should look something like
$$\phi(r,\theta) = r\theta \sin\theta + \frac{2P}{\pi d}r^2$$
- Integrate the strains
\begin{align}
u_x = \int\epsilon_{xx} dx + f(y)\\
u_y = \int\epsilon_{yy} dy + g(x)
\end{align}
with $2\epsilon_{xy} = \partial u_x/\partial x + \partial u_y/\partial y$, differentiate this equation w.r.t $y$ and $x$ and solve for $f$ and $g$.
References
- Puttock, M. J., & Thwaite, E. G. (1969). Elastic compression of spheres and cylinders at point and line contact. Melbourne, Australia: Commonwealth Scientific and Industrial Research Organization.
