Point source on top of a halfspace¶
To illustrate the use of the package we are going to play with the solutions for a concentrated force located on top of a halfspace. The origin, \(\mathbf{x} = (0,0,0)\), is placed on the free surface and positive \(z\) is inside the medium. This problem is of interest when modeling the deformation/stress around a localized load, e.g., the load caused by the weigth of a building on top of a soil.
The derivations for the strain and stress tensors are not too difficult, but it can get cumbersome really fast because of the lengthy calculations. Using the package we can simplify the whole process.
import numpy as np
from sympy import init_printing, symbols, lambdify, S, simplify
from sympy import pi, Matrix, sqrt, oo
from continuum_mechanics.solids import sym_grad, strain_stress
import matplotlib.pyplot as plt
from matplotlib import colors
The following snippet allows to format the graphs.
repo = "https://raw.githubusercontent.com/nicoguaro/matplotlib_styles/master"
style = repo + "/styles/minimalist.mplstyle"
plt.style.use(style)
x, y, z, r, E, nu, Fx, Fy, Fz = symbols('x y z r E nu F_x F_y F_z')
The components of the displacement vector are given by [LANDAU]
\[\begin{split}\begin{align}
&u_x = \frac{(1 + \nu)}{2 \pi E} \left\{\left[\frac{xz}{r^3} - \frac{(1 - 2\nu)x}{r(r + z)}\right]F_z +
\frac{2(1 - \nu)r +z}{r(r + z)}F_x +\frac{[2r(\nu r + z) + z^2]x}{r^3(r + z)^2}(xF_x + y F_y)\right\}\, ,\\
&u_y = \frac{(1 + \nu)}{2 \pi E} \left\{\left[\frac{yz}{r^3} - \frac{(1 - 2\nu)y}{r(r + z)}\right]F_z +
\frac{2(1 - \nu)r +z}{r(r + z)}F_y +\frac{[2r(\nu r + z) + z^2]y}{r^3(r + z)^2}(xF_x + y F_y)\right\}\, ,\\
&u_z = \frac{(1 + \nu)}{2 \pi E} \left\{\left[\frac{2(1 - \nu)}{r} - \frac{z^2}{r^3}\right]F_z +\left[\frac{1 - 2\nu}{r(r + z)} + \frac{z}{r^3}\right](xF_x + y F_y)\right\}\, ,
\end{align}\end{split}\]
with \(r = \sqrt{x^2 + y^2 + z^2}\).
ux = (1+nu)/(2*pi*E)*((x*z/r**3 - (1-2*nu)*x/(r*(r + z)))*Fz +
(2*(1 - nu)*r + z)/(r*(r + z))*Fx +
((2*r*(nu*r + z) + z**2)*x)/(r**3*(r + z)**2)*(x*Fx + y*Fy))
ux
\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(\frac{F_{x}}{r \left(r + z\right)} \left(r \left(- 2 \nu + 2\right) + z\right) + F_{z} \left(- \frac{x \left(- 2 \nu + 1\right)}{r \left(r + z\right)} + \frac{x z}{r^{3}}\right) + \frac{x}{r^{3} \left(r + z\right)^{2}} \left(F_{x} x + F_{y} y\right) \left(2 r \left(\nu r + z\right) + z^{2}\right)\right)\]
uy = (1+nu)/(2*pi*E)*((y*z/r**3 - (1-2*nu)*y/(r*(r + z)))*Fz +
(2*(1 - nu)*r + z)/(r*(r + z))*Fy +
((2*r*(nu*r + z) + z**2)*y)/(r**3*(r + z)**2)*(x*Fx + y*Fy))
uy
\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(\frac{F_{y}}{r \left(r + z\right)} \left(r \left(- 2 \nu + 2\right) + z\right) + F_{z} \left(- \frac{y \left(- 2 \nu + 1\right)}{r \left(r + z\right)} + \frac{y z}{r^{3}}\right) + \frac{y}{r^{3} \left(r + z\right)^{2}} \left(F_{x} x + F_{y} y\right) \left(2 r \left(\nu r + z\right) + z^{2}\right)\right)\]
uz = (1+nu)/(2*pi*E)*((2*(1 - nu)/r + z**2/r**3)*Fz +
((1 - 2*nu)/(r*(r + z)) + z/r**3)*(x*Fx + y*Fy))
uz
\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(F_{z} \left(\frac{1}{r} \left(- 2 \nu + 2\right) + \frac{z^{2}}{r^{3}}\right) + \left(F_{x} x + F_{y} y\right) \left(\frac{- 2 \nu + 1}{r \left(r + z\right)} + \frac{z}{r^{3}}\right)\right)\]
Withouth loss of generality we can assume that \(F_y=0\), this is equivalent a rotate the axes until the force is in the plane \(y=0\).
ux = ux.subs(Fy, 0)
ux
\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(\frac{Fx}{r \left(r + z\right)} \left(r \left(- 2 \nu + 2\right) + z\right) + \frac{Fx x^{2}}{r^{3} \left(r + z\right)^{2}} \left(2 r \left(\nu r + z\right) + z^{2}\right) + Fz \left(- \frac{x \left(- 2 \nu + 1\right)}{r \left(r + z\right)} + \frac{x z}{r^{3}}\right)\right)\]
uy = ux.subs(Fy, 0)
uy
\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(\frac{Fx}{r \left(r + z\right)} \left(r \left(- 2 \nu + 2\right) + z\right) + \frac{Fx x^{2}}{r^{3} \left(r + z\right)^{2}} \left(2 r \left(\nu r + z\right) + z^{2}\right) + Fz \left(- \frac{x \left(- 2 \nu + 1\right)}{r \left(r + z\right)} + \frac{x z}{r^{3}}\right)\right)\]
uz = uz.subs(Fy, 0)
uz
\[\frac{1}{2 \pi E} \left(\nu + 1\right) \left(Fx x \left(\frac{- 2 \nu + 1}{r \left(r + z\right)} + \frac{z}{r^{3}}\right) + Fz \left(\frac{1}{r} \left(- 2 \nu + 2\right) + \frac{z^{2}}{r^{3}}\right)\right)\]
The displacement vector is then
u = Matrix([ux, uy, uz]).subs(r, sqrt(x**2 + y**2 + z**2))
We can check that the displacement vanish when \(x,y,z \rightarrow \infty\)
u.limit(x, oo)
\[\begin{split}\left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]
u.limit(y, oo)
\[\begin{split}\left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]
u.limit(z, oo)
\[\begin{split}\left[\begin{matrix}0\\0\\0\end{matrix}\right]\end{split}\]
We can compute the strain and stress tensors using the symmetric
gradient
(vector.sym_grad())
and strain-to-stress
(solids.strain_stress())
functions.
lamda = E*nu/((1 + nu)*(1 - 2*nu))
mu = E/(2*(1 - nu))
strain = sym_grad(u)
stress = strain_stress(strain, [lamda, mu])
The expressions for strains and stresses are lengthy and difficult to work with. Nevertheless, we can work with them. For example, we can evaluate the stress tensor at a point \(\mathbf{x} = (1, 0, 1)\) for a vertical load and a Poisson coefficient \(\nu = 1/4\).
simplify(stress.subs({x: 1, y: 0, z:1, nu: S(1)/4, Fx: 0}))
\[\begin{split}\left[\begin{matrix}- \frac{F_{z} \left(73 \sqrt{2} + 108\right)}{48 \pi \left(7 + 5 \sqrt{2}\right)} & - \frac{5 F_{z} \left(2 \sqrt{2} + 3\right)}{24 \pi \left(7 + 5 \sqrt{2}\right)} & - \frac{5 F_{z} \left(4 + 3 \sqrt{2}\right)}{16 \pi \left(2 \sqrt{2} + 3\right)}\\- \frac{5 F_{z} \left(2 \sqrt{2} + 3\right)}{24 \pi \left(7 + 5 \sqrt{2}\right)} & - \frac{F_{z} \left(11 \sqrt{2} + 16\right)}{16 \pi \left(7 + 5 \sqrt{2}\right)} & 0\\- \frac{5 F_{z} \left(4 + 3 \sqrt{2}\right)}{16 \pi \left(2 \sqrt{2} + 3\right)} & 0 & - \frac{F_{z} \left(103 \sqrt{2} + 148\right)}{48 \pi \left(7 + 5 \sqrt{2}\right)}\end{matrix}\right]\end{split}\]
Visualization of the fields¶
Since it is difficult to handle these lengthy expressions we can visualize them. For that, we define a grid where to evaluate the expressions,
\[(x, z) \in [-2, 2]\times[0, 5]\, ,\]
in this case.
x_vec, z_vec = np.mgrid[-2:2:100j, 0:5:100j]
We can use lambdify() to turn the SymPy expressions to evaluatable functions.
def field_plot(expr, x_vec, y_vec, z_vec, E_val, nu_val, Fx_val, Fz_val, title=''):
"""Plot the field"""
# Lambdify the function
expr_fun = lambdify((x, y, z, E, nu, Fx, Fz), expr, "numpy")
expr_vec = expr_fun(x_vec, y_vec, z_vec, E_val, nu_val, Fx_val, Fz_val)
# Determine extrema
vmin = np.min(expr_vec)
vmax = np.max(expr_vec)
print("Minimum value in the domain: {:g}".format(vmin))
print("Maximum value in the domain: {:g}".format(vmax))
vmax = max(np.abs(vmax), np.abs(vmin))
# Plotting
fig = plt.gcf()
levels = np.logspace(-1, np.log10(vmax), 10)
levels = np.hstack((-levels[-1::-1], [0], levels))
cbar_ticks = ["{:.2g}".format(level) for level in levels]
cont = plt.contourf(x_vec, z_vec, expr_vec, levels=levels,
cmap="RdYlBu_r", norm=colors.SymLogNorm(0.1))
cbar = fig.colorbar(cont, ticks=levels[::2])
cbar.ax.set_yticklabels(cbar_ticks[::2])
plt.axis("image")
plt.gca().invert_yaxis()
plt.xlabel(r"$x$")
plt.ylabel(r"$z$")
plt.title(title)
return cont
Displacements¶
plt.figure()
field_plot(u.norm(), x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()
Minimum value in the domain: 0.0881197
Maximum value in the domain: 15.4645
plt.figure()
field_plot(u[0], x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()
Minimum value in the domain: -4.09665
Maximum value in the domain: 4.09665
plt.figure()
field_plot(u[2], x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()
Minimum value in the domain: 0.0869101
Maximum value in the domain: 14.3383
Stresses¶
We can plot the components of stress
for row in range(0, 3):
for col in range(row, 3):
plt.figure()
field_plot(stress[row,col], x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0,
title=r"$\sigma_{%i%i}$"%(row+1, col+1))
plt.show()
Minimum value in the domain: -41.4274
Maximum value in the domain: 406.682
Minimum value in the domain: -12.0021
Maximum value in the domain: 144.846
Minimum value in the domain: -95.9472
Maximum value in the domain: 95.9472
Minimum value in the domain: -59.0538
Maximum value in the domain: 116.991
Minimum value in the domain: -506.96
Maximum value in the domain: 506.96
Minimum value in the domain: -243.272
Maximum value in the domain: 116.991
Stress invariants¶
We can also plot the invariants of the stress tensor
I1 = S(1)/3 * stress.trace()
I2 = S(1)/2 * (stress.trace()**2 + (stress**2).trace())
I3 = stress.det()
Mises = sqrt(((stress[0,0] - stress[1,1])**2 + (stress[1,1] - stress[2,2])**2 +
(stress[2,2] - stress[0,0])**2 +
6*(stress[0,1]**2 + stress[1,2]**2 + stress[0,2]**2))/2)
plt.figure()
field_plot(I1, x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()
Minimum value in the domain: -107.797
Maximum value in the domain: 213.555
plt.figure()
field_plot(I2, x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()
Minimum value in the domain: 0.000977492
Maximum value in the domain: 579596
plt.figure()
field_plot(I3, x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()
Minimum value in the domain: -1.01409e+08
Maximum value in the domain: 419218
plt.figure()
field_plot(Mises, x_vec, 0, z_vec, 1.0, 0.3, 0.0, 1.0)
plt.show()
Minimum value in the domain: 0.0274784
Maximum value in the domain: 958.065