[SOLVED] CS elasticity1d

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. # Youngsmodulus
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.]