elasticity1d
FEA¶
Worked example – linear elastic bar in 1D
Import some libraries needed for numerical computing and plotting
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
Function to compute the element stiffness matrix $\mathbf{K}^e =\dfrac{A^eE^e}{l^e}
\begin{bmatrix}
1 & -1 \\
-1 & 1
\end{bmatrix}$
In [ ]:
def get_Ke(Ae, Ee, le):
return (Ae * Ee / le) * np.array([[1., -1.], [-1., 1.]])
Function to compute the element body force vector $\mathbf{f}^e_\Omega = \dfrac{l^e}{6}
\begin{bmatrix}
2 & 1 \\
1 & 2
\end{bmatrix}
\begin{bmatrix}
b_1^e \\ b_2^e
\end{bmatrix}$
In [ ]:
def get_fe_omega(le, b1, b2):
return (le / 6.) * np.array([[2., 1.], [1., 2.]]).dot(np.array([[b1], [b2]]))
Function to return the global dofs for the element
In [ ]:
def get_dof_index(e):
return [e, e + 1]
Compute the strain from the displacement, $\dfrac{\mathsf{d}u^e(x)}{\mathsf{d} x} = \mathbf{B}^e \mathbf{d}^e$
In [ ]:
def get_strain_e(le, de):
return (1. / le) * np.array([-1., 1.]).dot(de)
Define the material properties and the loading
In [ ]:
E = 8. # Young’s modulus
l = 4. # length of bar
A = 2. # cross sectional area
b = 3. # body force
t_bar = 1. # applied traction
Determine the coordinates of the nodes $x$ based on the number of elements. Assume uniform spacing of the nodes.
In [ ]:
n_el = 5 # number of elements
n_np = n_el + 1 # number of nodal points
x = np.linspace(0, l, n_np) # nodal coordinates
le = l / n_el # length of element
Initialise the global stiffness matrix $\mathbf{K}$ and the force vector $\mathbf{f}$
In [ ]:
K = np.zeros((n_np, n_np))
f = np.zeros((n_np, 1))
Loop over all the elements and assemble $\mathbf{K}$ and $\mathbf{f}$ from the element contributions
In [ ]:
for ee in range(n_el):
dof_index = get_dof_index(ee)
K[np.ix_(dof_index, dof_index)] += get_Ke(A, E, le)
f[np.ix_(dof_index)] += get_fe_omega(le, b, b )
Add the applied traction
In [ ]:
f[-1] += t_bar * A
Apply the boundary conditions and solve the problem
In [ ]:
free_dof = np.arange(1,n_np)
K_free = K[np.ix_(free_dof, free_dof)]
f_free = f[np.ix_(free_dof)]
# solve the linear system
d_free = np.linalg.solve(K_free,f_free)
d = np.zeros((n_np, 1))
d[1:] = d_free
Postprocessing.
Compute the reaction force.
Display the structure of $\mathbf{K}$.
Plot the stress distribution.
The exact solution is given by $u^{\text{ex}} = -\dfrac{3}{32}x^2 + \dfrac{7}{8}x$ and
$\sigma^\text{ex} = -\dfrac{3}{2}x + 7 = \bar{t} + \dfrac{b(l-x)}{A}$.
In [ ]:
r = K[0,:].dot(d) – f[0]
print(‘Reaction force ‘, r)
x_ex = np.linspace(0, l, 100*n_el) # nodal coordinates
u_ex = -(3/32)*np.power(x_ex,2) + (7/8)*x_ex
stress_ex = (-3/2)*x_ex + 7
line1, = plt.plot(x, d, label = ‘$u^h$’, marker = ‘o’)
line2, = plt.plot(x_ex, u_ex, label = ‘$u^\mathrm{ex}$’)
plt.xlabel(‘$x$’)
plt.ylabel(‘$u$’)
plt.legend(handles=[line1, line2])
plt.grid(True)
plt.show()
plt.spy(K)
plt.show()
x_post = np.zeros((2 * n_el))
x_post[0::2] = x[0:-1]
x_post[1::2] = x[1:]
stress_post = np.zeros((2 * n_el))
# compute the strain in each element
for ee in range(n_el):
dof_index = get_dof_index(ee)
index = [2 * ee, 2 * ee + 1]
de = d[dof_index]
stress_e = E * get_strain_e(le, de)
stress_post[index] = stress_e
line1, = plt.plot(x_post, stress_post, label = ‘$\sigma^h$’, marker = ‘o’)
line2, = plt.plot(x_ex, stress_ex, label = ‘$\sigma^\mathrm{ex}$’)
plt.legend(handles=[line1, line2])
plt.xlabel(‘$x$’)
plt.ylabel(‘$\sigma$’)
plt.grid(True)
plt.show()
Reaction force [-14.]