This notebook covers the following aspects: * Define the concept of static elasticity * Solution using finite element scheme * Solution using the finite difference scheme and its comparion with FEM solution
Static elasticity can be considered a particular case derived from the elastic wave equation when the displacement does not depend on time, i.e \(\partial_t^2 u(x,t) = 0\). Under this assumption and departing from the 1D elastic wave equation, the differential equation turns into the Poisson equation
\[\begin{equation} -\mu \partial_x^2 u = f, \end{equation}\]
where \(\mu\) is the shear modulus for a homogeneous media, \(u\) is the displacement field, and \(f\) is the external force. The solution for this problems is found after bringing this equation into its weak form, applying the free boundary condition, and using the Galerkin principle with a suitable basis. Then, the displacement defined in a discrete set of points \(x_i\) is given as the solution of a system of N equations, with
\[\begin{equation} \mathbf{u} = (\mathbf{K}^{T})^{-1} \mathbf{f} \end{equation}\]
where \(\mathbf{K}\) is the stiffness matrix. For an elastic physical system with constant shear modulus \(\mu\) and uniform element size \(h\), it is given as
\[\begin{equation} K_{ij} = \frac{\mu}{h} \begin{pmatrix} 1 & -1 & & & \\ -1 & 2 & -1 & & \\ & & \ddots & & \\ & & -1 & 2 & -1 \\ & & & -1 & 1 \end{pmatrix} \end{equation}\]
The purpose of this notebook is to illustrate how the problem of static elasticity is solved with the finite- element method. We also compare the solution using finite - differences, the so-called relaxation method given by (derivation in the book):
\[\begin{equation} u_{i}^{k+1} = \dfrac{u_{i+1}^{k} + u_{i-1}^{k}}{2} + \dfrac{h^2}{2 \mu}f_i \end{equation}\]
# Import all necessary libraries, this is a configuration step for the exercise.
# Please run it before the simulation code!
import numpy as np
import matplotlib
# Show Plot in The Notebook
matplotlib.use("nbagg")
import matplotlib.pyplot as plt# ---------------------------------------------------------------
# Initialization of setup
# ---------------------------------------------------------------
nx = 20 # Number of element boundary points
u = np.zeros(nx) # Solution vector
f = np.zeros(nx) # Source vector
mu = 1 # Constant shear modulus
# Element boundary points
x = np.linspace(0, 1, nx) # x in [0,1]
h = x[2] - x[1] # Constant element size
# ---------------------------------------------------------------
# Assemble stiffness matrix K_ij
# ---------------------------------------------------------------
K = np.zeros((nx, nx))
for i in range(1, nx-1):
for j in range(1, nx-1):
if i == j:
K[i, j] = 2*mu/h
elif i == j + 1:
K[i, j] = -mu/h
elif i + 1 == j:
K[i, j] = -mu/h
else:
K[i, j] = 0
# ---------------------------------------------------------------
# Souce term is a spike at i = 3*nx/4
f[int(3*nx/4)] = 1
# Boundary condition at x = 0
u[0] = 0.15 ; f[1] = u[0]/h
# Boundary condition at x = 1
u[nx-1] = 0.05 ; f[nx-2] = u[nx-1]/h
# ---------------------------------------------------------------
# Finite element solution
# ---------------------------------------------------------------
u[1:nx-1] = np.linalg.inv(K[1:nx-1, 1:nx-1]) @ np.transpose(f[1:nx-1])
# ---------------------------------------------------------------
# Plotting section
# ---------------------------------------------------------------
plt.figure(figsize=(8,6))
xfe = u
plt.plot(x, xfe, color='r', lw=2, label='Finite elements')
plt.title('Static Elasticity', size=16)
plt.ylabel('Displacement $u(x)$', size=16)
plt.xlabel('Position $x$', size=16)
plt.axis([0, 1, 0.04, .28])
plt.legend()
plt.grid(True)
plt.show()
# Poisson's equation with relaxation method
# ---------------------------------------------------------------
nt = 500 # Number of iterations
iplot = 2 # Snapshot frequency
# non-zero boundary conditions
u = np.zeros(nx) # set u to zero
du = np.zeros(nx) # du/dx
f = np.zeros(nx) # forcing
f[int(3*nx/4)] = 1./h
xfd = np.arange(0, nx)*h
# ---------------------------------------------------------------
# Initialize animated plot
# ---------------------------------------------------------------
plt.figure(figsize=(8,6))
line1 = plt.plot(x, xfe, color='r', lw=2, label='FE')
line2 = plt.plot(xfd, u, color='k', ls='-.', label='FD relaxation')
plt.title('Static Elasticity with relaxation method', size=16)
plt.ylabel('Displacement, $u$', size=12)
plt.xlabel('Position, $x$', size=12)
plt.legend(loc=4)
plt.grid(True)
plt.ion() # set interective mode
plt.show()
# ---------------------------------------------------------------
for it in range(nt):
# Calculate the average of u (omit boundaries)
for i in range(1, nx-1):
du[i] =u [i+1] + u[i-1]
u = 0.5*( f*h**2/mu + du )
u[0] = 0.15 # Boundary condition at x=0
u[nx-1] = 0.05 # Boundary condition at x=1
fd = u
# --------------------------------------
# Animation plot. Display both solutions
if not it % iplot:
for l in line2:
l.remove()
del l
line1 = plt.plot(x, xfe, color='r', lw=2)
line2 = plt.plot(xfd, fd, color='k', ls='-.')
plt.gcf().canvas.draw() 