CoCalc Shared Files2015-10-23-140959.html
Author: Jason Moore
Views : 12
Description: Jupyter html version of 2015-10-23-140959.ipynb
2015-10-23-140959

Plot the ratio of the principal normal and shear stresses to the maximum contact area pressure as a function of the depth below the surface along the z axis for two spheres.

In [51]:
from sympy import symbols, Eq, lambdify, init_printing, pi, Abs, atan, S

In [52]:
from IPython.html.widgets import interactive
from IPython.display import display

In [53]:
import matplotlib.pyplot as plt
import numpy as np

In [54]:
%matplotlib nbagg

In [55]:
init_printing()


First define all of the known and unknown symbols needed for the calculations for sphere to sphere contact.

In [56]:
a, F, v, v1, v2, E1, E2, d1, d2 = symbols('a, F, nu, nu1, nu2, E1, E2, d1, d2', real=True)
p_max, s1, s3, t13, z = symbols('p_max, sigma1, sigma3, tau_13, z', real=True)


Now define the equations for the contact area radius, $a$, the maximum pressure at the contact area, $p_{max}$, and the principal stresses along the $z$ axis.

In [57]:
radius_eq = Eq(a, (3 * F / 8 * ((1 - v1**2) / E1 + (1 - v2**2)) / (1 / d1 + 1 / d2))**3)

Out[57]:
a == 27*F**3*(-nu2**2 + 1 + (-nu1**2 + 1)/E1)**3/(512*(1/d2 + 1/d1)**3)
In [58]:
pressure_eq = Eq(p_max, 3 * F / 2 / pi / a**2)
pressure_eq

Out[58]:
p_max == 3*F/(2*pi*a**2)
In [59]:
s1_eq = Eq(s1, -p_max * ((1 - Abs(z / a) * atan(1 / Abs(z / a))) * (1 + v) - S(1) / 2 / (1 + z**2 / a**2)))
s1_eq

Out[59]:
sigma1 == -p_max*((nu + 1)*(-Abs(z/a)*atan(1/Abs(z/a)) + 1) - 1/(2*(1 + z**2/a**2)))
In [60]:
s3_eq = Eq(s3, -p_max / (1 + z**2 / a**2))
s3_eq

Out[60]:
sigma3 == -p_max/(1 + z**2/a**2)
In [61]:
t_max_eq = Eq(t13, (s1 - s3) / 2)
t_max_eq

Out[61]:
tau_13 == sigma1/2 - sigma3/2
In [74]:
t_max_eq.rhs.subs({s1: s1_eq.rhs, s3: s3_eq.rhs}).diff(z).simplify()

Out[74]:
-p_max*(3*a**3*z + (a**2 + z**2)*(nu + 1)*(a**2*Abs(z/a) - (a**2 + z**2)*atan(1/Abs(z/a)))*sign(z/a))/(2*a*(a**2 + z**2)**2)
In [63]:
(s1_eq.rhs / pressure_eq.rhs).simplify()

Out[63]:
pi*a**2*p_max*(a**2 + 2*(a**2 + z**2)*(nu + 1)*(Abs(z/a)*atan(1/Abs(z/a)) - 1))/(3*F*(a**2 + z**2))

Create some functions to evaluate the expressions that we derived above.

In [64]:
f_a = lambdify((F, v1, v2, E1, E2, d1, d2), radius_eq.rhs, modules='numpy')

In [65]:
f_p_max = lambdify((F, a), pressure_eq.rhs, modules='numpy')

In [66]:
f_s1 = lambdify((p_max, z, a, v), s1_eq.rhs, modules='numpy')

In [67]:
f_s3 = lambdify((p_max, z, a), s3_eq.rhs, modules='numpy')

In [68]:
f_t_max = lambdify((s1, s3), t_max_eq.rhs, modules='numpy')


We need a graph of $\frac{\sigma_1, \sigma_3, \tau_{1/2}}{p_{max}}$ vs $z$ so define a function that produces the plot.

In [69]:
def plot_stresses(max_depth=3.0, v1=0.326, v2=0.324, E1=17.2E6, E2=15.4E6, d1=1.0, d2=2.0):
"""Plots the principal stresses in each sphere versus the depth along the z axis.

Parameters
===========
max_depth : float
The depth is defined as z = max_depth * a.
v1, v2 : float
Poisson's ratio for sphere 1 and 2.
E1, E2 : float
Modulus of elasticity for sphere 1 and 2.
d1, d2 : float
Diameters of sphere 1 and 2.

"""

fig, axes = plt.subplots(1, 2, sharey=True)

F = 1.0

a = f_a(F, v1, v2, E1, E2, d1, d2)
z = np.linspace(1E-14, max_depth * a, num=100)
p_max = f_p_max(F, a)

s1_1 = f_s1(p_max, z, a, v1)
s3_1 = f_s3(p_max, z, a)
t_13_1 = f_t_max(s1_1, s3_1)

s1_2 = f_s1(p_max, z, a, v2)
s3_2 = s3_1
t_13_2 = f_t_max(s1_2, s3_2)

z = z / a

axes[0].plot(z, np.abs(s1_1 / p_max),
z, np.abs(s3_1 / p_max),
z, np.abs(t_13_1/ p_max))
axes[1].plot(z, np.abs(s1_2 / p_max),
z, np.abs(s3_2 / p_max),
z, np.abs(t_13_2/ p_max))

axes[0].legend([r'$\sigma_1$', r'$\sigma_2$', r'$\tau_{max}$'])
axes[1].legend([r'$\sigma_1$', r'$\sigma_2$', r'$\tau_{max}$'])

axes[0].set_xlabel('$z/a$')
axes[1].set_xlabel('$z/a$')

axes[0].set_ylabel(r'$\frac{\sigma_1, \sigma_3, \tau_{1/2}}{p_{max}}$')

axes[0].set_title('Sphere 1')
axes[1].set_title('Sphere 2')

plt.tight_layout()

plt.show()

In [70]:
plot_stresses()

In [71]:
w = interactive(plot_stresses, v1=(0.2, 0.5), v2=(0.2, 0.5), E1=(0.0, 50.0E6),
E2=(0.0, 50.0E6), d1=(-10.0, 10.0), d2=(-10.0, 10.0))

In [72]:
display(w)

None
In [ ]: