Griffith

In [1]:

import sympy
import numpy as np
import matplotlib.pyplot as plt

sympy.init_printing()

import sympy import numpy as np import matplotlib.pyplot as plt sympy.init_printing()

According to Linear elastic mechanics, ahead of a mode-I crack tip, the stress distribution is given by: $$ \sigma_x = \frac{K_I}{\sqrt{2\pi r}}\cos\left(\frac{\theta}{2}\right)\left[1 - \sin\left(\frac{\theta}{2}\right)\sin\left(\frac{3\theta}{2}\right)\right] $$ $$ \sigma_y = \frac{K_I}{\sqrt{2\pi r}}\cos\left(\frac{\theta}{2}\right)\left[1 + \sin\left(\frac{\theta}{2}\right)\sin\left(\frac{3\theta}{2}\right)\right] $$ $$\tau_{xy} = \frac{K_I}{\sqrt{2\pi r}}\cos\left(\frac{\theta}{2}\right)\sin\left(\frac{\theta}{2}\right)\cos\left(\frac{3\theta}{2}\right)$$ $$\tau_{xz} = \tau_{yz} = 0$$ For plane stress condition: $$\sigma_z = 0,$$ for plane strain condition: $$\sigma_z = \nu(\sigma_x + \sigma_y).$$

In [2]:

Y, sigma, a, r, theta, sigma_s = sympy.symbols("Y, sigma, a, r, theta, sigma_s", positive = True)
#
KI = Y*sigma*sympy.sqrt(sympy.pi*a)

# Plane stress condition
sigma_x = KI/sympy.sqrt(2*sympy.pi*r)*sympy.cos(theta/2) * ( 1 - sympy.sin(theta/2)*sympy.sin(3.*theta/2))
sigma_y = KI/sympy.sqrt(2*sympy.pi*r)*sympy.cos(theta/2) * ( 1 + sympy.sin(theta/2)*sympy.sin(3.*theta/2))
sigma_z = 0.
tau_xy  = KI/sympy.sqrt(2*sympy.pi*r)*sympy.cos(theta/2) * sympy.sin(theta/2)*sympy.cos(3.*theta/2)
# tau_xz = tau_yz = 0, skipped in the equations below

# Von Mises stress
sigma_v = sympy.sqrt(0.5*((sigma_x - sigma_y)**2 + (sigma_x - sigma_z)**2 + (sigma_y-sigma_z)**2 + 6.*tau_xy**2))

Y, sigma, a, r, theta, sigma_s = sympy.symbols("Y, sigma, a, r, theta, sigma_s", positive = True) # KI = Y*sigma*sympy.sqrt(sympy.pi*a) # Plane stress condition sigma_x = KI/sympy.sqrt(2*sympy.pi*r)*sympy.cos(theta/2) * ( 1 - sympy.sin(theta/2)*sympy.sin(3.*theta/2)) sigma_y = KI/sympy.sqrt(2*sympy.pi*r)*sympy.cos(theta/2) * ( 1 + sympy.sin(theta/2)*sympy.sin(3.*theta/2)) sigma_z = 0. tau_xy = KI/sympy.sqrt(2*sympy.pi*r)*sympy.cos(theta/2) * sympy.sin(theta/2)*sympy.cos(3.*theta/2) # tau_xz = tau_yz = 0, skipped in the equations below # Von Mises stress sigma_v = sympy.sqrt(0.5*((sigma_x - sigma_y)**2 + (sigma_x - sigma_z)**2 + (sigma_y-sigma_z)**2 + 6.*tau_xy**2))

Stress distribution along x direction ($\theta = 0$), and $y=x$ direction ($\theta = \pi/4$):

In [3]:

r_list = np.linspace(0.01, 5, 1000)
sx_list = list()
sy_list = list()
tau_list = list()
for ri in r_list:
    sx_list.append(sigma_x.subs({Y:1, a:1, sigma:1, r:ri, theta:0}))
    sy_list.append(sigma_y.subs({Y:1, a:1, sigma:1, r:ri, theta:0}))
    tau_list.append(tau_xy.subs({Y:1, a:1, sigma:1, r:ri, theta:0}))
plt.plot(r_list, sx_list, r_list, sy_list, r_list, tau_list)

r_list = np.linspace(0.01, 5, 1000) sx_list = list() sy_list = list() tau_list = list() for ri in r_list: sx_list.append(sigma_x.subs({Y:1, a:1, sigma:1, r:ri, theta:0})) sy_list.append(sigma_y.subs({Y:1, a:1, sigma:1, r:ri, theta:0})) tau_list.append(tau_xy.subs({Y:1, a:1, sigma:1, r:ri, theta:0})) plt.plot(r_list, sx_list, r_list, sy_list, r_list, tau_list)

Out[3]:

[<matplotlib.lines.Line2D at 0x243bead4580>,
 <matplotlib.lines.Line2D at 0x243bead4520>,
 <matplotlib.lines.Line2D at 0x243bead4640>]

In [4]:

sx_list = list()
sy_list = list()
tau_list = list()
for ri in r_list:
    sx_list.append(sigma_x.subs({Y:1, a:1, sigma:1, r:ri, theta:np.pi/4}))
    sy_list.append(sigma_y.subs({Y:1, a:1, sigma:1, r:ri, theta:np.pi/4}))
    tau_list.append(tau_xy.subs({Y:1, a:1, sigma:1, r:ri, theta:np.pi/4}))
plt.plot(r_list, sx_list, r_list, sy_list, r_list, tau_list)

sx_list = list() sy_list = list() tau_list = list() for ri in r_list: sx_list.append(sigma_x.subs({Y:1, a:1, sigma:1, r:ri, theta:np.pi/4})) sy_list.append(sigma_y.subs({Y:1, a:1, sigma:1, r:ri, theta:np.pi/4})) tau_list.append(tau_xy.subs({Y:1, a:1, sigma:1, r:ri, theta:np.pi/4})) plt.plot(r_list, sx_list, r_list, sy_list, r_list, tau_list)

Out[4]:

[<matplotlib.lines.Line2D at 0x243c1da7df0>,
 <matplotlib.lines.Line2D at 0x243c1da7e50>,
 <matplotlib.lines.Line2D at 0x243c1da7f70>]

According to the von Mises yield criterion, when $$\sigma_v \geq \sigma_s$$ yielding occurs, where the von Mises stress for a general stress state is given by $$\sigma_v = \sqrt{\frac{1}{2}\left[(\sigma_x-\sigma_y)^2 + (\sigma_x - \sigma_z)^2 + (\sigma_y - \sigma_z)^2 + 6(\gamma_{xy}^2 + \gamma_{xz}^2 + \gamma_{yz}^2)\right]}.$$

In [5]:

# Let Sigma_{von_mises} = Sigma_s ()
rp = sympy.solve(sigma_v - sigma_s, r)[0]
rp.simplify()

# Let Sigma_{von_mises} = Sigma_s () rp = sympy.solve(sigma_v - sigma_s, r)[0] rp.simplify()

Out[5]:

$\displaystyle \frac{1.0 Y^{2} a \sigma^{2} \left(1.25 - 0.750000000000001 \cos{\left(\theta \right)}\right) \cos^{2}{\left(0.5 \theta \right)}}{\sigma_{s}^{2}}$

Along x direction （within the crack plane), $$r_p = \frac{1}{2\pi}\left(\frac{K_I}{\sigma_s}\right)^2.$$

In [6]:

# Assuming Y = 1, a = 1, sigma_1, and sigma_s = 1, and plot the rp as a function of theta
angles = np.linspace(0, 2.*np.pi, 360)
rp_1 = list()
for ang in angles:
    rp_1.append(rp.subs({Y:1, sigma:1, a:1, sigma_s:1, theta:ang}))

# Assuming Y = 1, a = 1, sigma_1, and sigma_s = 1, and plot the rp as a function of theta angles = np.linspace(0, 2.*np.pi, 360) rp_1 = list() for ang in angles: rp_1.append(rp.subs({Y:1, sigma:1, a:1, sigma_s:1, theta:ang}))

In [7]:

plt.axes(projection='polar')
plt.polar(angles, rp_1)
plt.show()

plt.axes(projection='polar') plt.polar(angles, rp_1) plt.show()

For plane strain conditions, assuming $\nu = 1/3$, we have:

In [8]:

nu = 1./3.
sigma_z = nu*(sigma_x + sigma_y)

sigma_v2 = sympy.sqrt(0.5*((sigma_x - sigma_y)**2 + (sigma_x - sigma_z)**2 + (sigma_y - sigma_z)**2 + 6.*tau_xy**2))

rp_plane_strain = sympy.solve(sigma_v2 - sigma_s, r)[0]
rp_plane_strain.simplify().subs({theta:0})

nu = 1./3. sigma_z = nu*(sigma_x + sigma_y) sigma_v2 = sympy.sqrt(0.5*((sigma_x - sigma_y)**2 + (sigma_x - sigma_z)**2 + (sigma_y - sigma_z)**2 + 6.*tau_xy**2)) rp_plane_strain = sympy.solve(sigma_v2 - sigma_s, r)[0] rp_plane_strain.simplify().subs({theta:0})

Out[8]:

$\displaystyle \frac{0.0555555555555556 Y^{2} a \sigma^{2}}{\sigma_{s}^{2}}$

In [9]:

rp_2 = list()
for ang in angles:
    rp_2.append(rp_plane_strain.subs({Y:1, sigma:1, a:1, sigma_s:1, theta:ang}))

rp_2 = list() for ang in angles: rp_2.append(rp_plane_strain.subs({Y:1, sigma:1, a:1, sigma_s:1, theta:ang}))

In [10]:

plt.axes(projection='polar')
plt.polar(angles, rp_1, angles, rp_2, '--')
plt.show()

plt.axes(projection='polar') plt.polar(angles, rp_1, angles, rp_2, '--') plt.show()