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}
egin{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}
egin{bmatrix}
2 & 1 \
1 & 2
end{bmatrix}
egin{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^{ ext{ex}} = -dfrac{3}{32}x^2 + dfrac{7}{8}x$ and
$sigma^ ext{ex} = -dfrac{3}{2}x + 7 = ar{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.]