Before we start, a few resources. A general intro on structural engineering theory. Why FEA simulations almost always underpredict residual strain in bending. And a video explainer on the Rayleigh Ritz Method for a beam.
The hat function is introduced on page 109 of the textbook. The first step is to learn how to write this in LaTex, with the example of equation (10.8) from the book:
$$\varphi_i = \begin{cases} \frac{x-x_{i-1}}{x_i-x_{i-1}} & \quad x_{i-1} \leq x < x_i \\ \frac{x_{i+1} - x}{x_{i+1}-x_{i}} & \quad x_i \leq x < x_{i+1} \\ 0 & \quad x < x_{i-1} \text{or} x \geq x_{i+1} \end{cases} $$# This is the code:
$$\varphi_i =
\begin{cases}
\frac{x-x_{i-1}}{x_i-x_{i-1}} & \quad x_{i-1} \leq x < x_i \\
\frac{x_{i+1} - x}{x_{i+1}-x_{i}} & \quad x_i \leq x < x_{i+1} \\
0 & \quad x < x_{i-1} \text{or} x \geq x_{i+1}
\end{cases}
$$
We're assuming that both boundaries are fixed at zero. $n$ is the number of hat functions, that we will place at an evenly spaced position at $1/h, ... , n/h$, with $h=1/(n+1)$. Recreating the hat function in image 10.1:
For $1 \leq i \leq n$ we will get the following:
$$\varphi_i (x) = \begin{cases} \frac{x-(i-1)h}{h} & \quad \text{for } (i-1)h < x < ih \\ \frac{(i+1)h-x}{h} & \quad \text{for } ih \leq x < (i+1)h \\ 0 & \quad \text{otherwise} \end{cases} $$$$\frac{d\varphi_i (x)}{dx} = \begin{cases} \frac{1}{h} & \quad \text{for } (i-1)h < x < ih \\ \frac{-1}{h} & \quad \text{for } ih \leq x < (i+1)h \\ 0 & \quad \text{otherwise} \end{cases} $$In the previous question a) we found that $\textbf{A}$ is:
$$\textbf{A}_{i,j} = \int_0^1 \varphi_i \varphi_j dx$$Inverting a matrix can be done by hand, or coded with Python. This is and example of how to invert a matrix with NumPy:
import numpy as np
# Taking a 3 * 3 matrix
A = np.array([[6, 1, 1],
[4, -2, 5],
[2, 8, 7]])
# Calculating the inverse of the matrix
print(np.linalg.inv(A))
[[ 0.17647059 -0.00326797 -0.02287582] [ 0.05882353 -0.13071895 0.08496732] [-0.11764706 0.1503268 0.05228758]]
Next, let's replicate the matrix of 10.21 with LaTex. A useful guide on creating large brackets, parentheses etc. for matrices can be found here. This code:
$$\begin{pmatrix}
u_0 \\
\dot u_0 \\
u_h \\
\dot u_h
\end{pmatrix}
=
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
1 & h & h^2 & h^3 \\
0 & 1 & 2h & 3h^2
\end{bmatrix}
\begin{pmatrix}
a_0 \\
a_1 \\
a_2 \\
a_3
\end{pmatrix}
$$
Returns this:
$$\begin{pmatrix} u_0 \\ \dot u_0 \\ u_h \\ \dot u_h \end{pmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & h & h^2 & h^3 \\ 0 & 1 & 2h & 3h^2 \end{bmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{pmatrix} $$Next, let's make a matrix with SymPy. The documentation with syntax can be found here.
import sympy as sp
x, x_i, h, u_1, dot_u_i, u_j, dot_u_j = sp.symbols('x x_i h u_1 \dot{u_i} u_j \dot{u_j}', real=True)
M=sp.Matrix([
[1, 0, 0, 0],
[0, 1, 0, 0],
[1, h, h**2, h**3],
[0, 1, 0, 0],
])
print(M)
Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [1, h, h**2, h**3], [0, 1, 0, 0]])
The 1D beam equation in (10.29) is:
$$V = \int^L_0 \left( \frac{1}{2} EI \left( \frac{d^2u}{dx^2} \right) ^2 - u(x)f(x) \right) dx$$Where $u(x)$ is the displacement and $f(x)$ the work you do to move the beam.
The value for $EI$ is stiffness. $E$ is the modulus of elasticity which depends on the chosen material. $I$ is the second moment of inertia in the beam's section, and can be calculated
$E=N/m^2$
$I = m^4 (check)$
Solve for the displacement of the function.
a + phi are basis functions
Function 10.30 but without the time component $$u(x)= \sum a_i \psi (x) $$
Represent basis functions, because the Hermite function Choose overall interval, e.g. start with 0 to 1, than scale x. Start with 0.1 to test, later scale up.
Checkpoints:
Coefficients only apperar.
Because my undergrad is in structural engineering at a university of applied sciences, I sort of expected this week to be somewhat familiar. It couldn't be further from the truth! Below are a few examples of what I'd thought we'd cover this week.
During undergrad, I was taught to create a beam by hand, and draw diagrams with normal forces, shear forces, bending moments and deflection. Next, we'd use software that calculates it for us. You basically press a few buttons, and that's about it. I'd ask: 'What's behind this equation?' and 'What are the assumptions for this value?' and was told that I should go to Delft University. I'm really happy to finally have this opportunity within this class. It turns out, the whole point of this week is to create those solvers. And indeed, that's something completely different than what I've learned so far. Below are two solvers: IntermediateBeam and PolymerFEM 2D.
I found a package called IndeterminateBeam on Gitlab that does this in Python. To me, this is quite a revelation, and really helpful in further work.
Go topypi for quickly installing packages.
Install (includes sub-packages including sympy) https://pypi.org/project/indeterminatebeam/
This works in python.
Initially, the outputs were not showing up in my notebook, but it did work when I'd run the script from my terminal. In the browser, it shows controls, not visible in JupyterLab. ipympl can be installed for showing controls in Jupyter lab.
pip3 install ipympl
Then it shows up.
%matplotlib inline
from indeterminatebeam import Beam, Support, PointLoadV, PointTorque, UDLV, DistributedLoadV
# Create 7 m beam with E, I, A as defaults
beam = Beam(1.5)
# Create 9 m beam with E, I, and A assigned by user
beam_2 = Beam(9, E=2000, I =10**6, A = 3000)
from indeterminatebeam import Support
# Defines a pin support at location x = 0 m
a = Support(0, (1,1,1))
# Defines a roller support at location x = 0 m
# b = Support(1.5, (0,1,0))
# Defines a fixed support at location x = 7 m
# c = Support(7, (1,1,1))
# Assign the support objects to a beam object created earlier
beam.add_supports(a)
from indeterminatebeam import PointLoadV, PointTorque, DistributedLoadV
# Create a 1000 N point load at x = 1.5 m
load_1 = PointLoadV(-1000, 1.5)
# Create a 500 N/m UDL from x = 1 m to x = 1.5 m
# load_2 = DistributedLoadV(-200, (0, 1.5))
# Defines a 2 kN.m point torque at x = 3.5 m
# load_3 = PointTorque(2*10**3, 3.5)
# Assign the load objects to the beam object
beam.add_loads(load_1)
beam.analyse()
fig_1 = beam.plot_beam_external()
fig_1.show()
fig_2 = beam.plot_beam_internal()
fig_2.show()
# print(fig_1)
I found this former MIT student who created a
https://www.youtube.com/watch?v=1j_HdsVkglk https://polymerfem.com/full-finite-element-solver-in-100-lines-of-python/
First, I'll annotate the code example to understand its mechanics.
# Copyright 2022 Jorgen Bergstrom
# See also code from http://compmech.lab.asu.edu/codes.php
# Additions by Nicole Bakker (in the comments)
import numpy as np
import math
from matplotlib import pyplot as plt
def shape(xi):
x,y = tuple(xi)
N = [(1.0-x)*(1.0-y), (1.0+x)*(1.0-y), (1.0+x)*(1.0+y), (1.0-x)*(1.0+y)]
return 0.25*np.array(N)
def gradshape(xi):
x,y = tuple(xi)
dN = [[-(1.0-y), (1.0-y), (1.0+y), -(1.0+y)],
[-(1.0-x), -(1.0+x), (1.0+x), (1.0-x)]]
return 0.25*np.array(dN)
fig = plt.figure() # Adding this for labels
ax = fig.add_subplot(111) # Adding this for labels
print('create mesh')
# input
mesh_ex = 9
mesh_ey = 49
mesh_lx = 10.0
mesh_ly = 50.0
# derived
mesh_nx = mesh_ex + 1
mesh_ny = mesh_ey + 1
num_nodes = mesh_nx * mesh_ny
num_elements = mesh_ex * mesh_ey
mesh_hx = mesh_lx / mesh_ex
mesh_hy = mesh_ly / mesh_ey
nodes = []
for y in np.linspace(0.0, mesh_ly, mesh_ny):
for x in np.linspace(0.0, mesh_lx, mesh_nx):
nodes.append([x,y])
nodes = np.array(nodes)
conn = []
for j in range(mesh_ey):
for i in range(mesh_ex):
n0 = i + j*mesh_nx
conn.append([n0, n0 + 1, n0 + 1 + mesh_nx, n0 + mesh_nx])
print ('material model - plane strain')
E = 100.0
v = 0.48
C = E/(1.0+v)/(1.0-2.0*v) * np.array([[1.0-v, v, 0.0],
[ v, 1.0-v, 0.0],
[ 0.0, 0.0, 0.5-v]])
print('create global stiffness matrix')
K = np.zeros((2*num_nodes, 2*num_nodes))
q4 = [[x/math.sqrt(3.0),y/math.sqrt(3.0)] for y in [-1.0,1.0] for x in [-1.0,1.0]]
B = np.zeros((3,8))
for c in conn:
xIe = nodes[c,:]
Ke = np.zeros((8,8))
for q in q4:
dN = gradshape(q)
J = np.dot(dN, xIe).T
dN = np.dot(np.linalg.inv(J), dN)
B[0,0::2] = dN[0,:]
B[1,1::2] = dN[1,:]
B[2,0::2] = dN[1,:]
B[2,1::2] = dN[0,:]
Ke += np.dot(np.dot(B.T,C),B) * np.linalg.det(J)
for i,I in enumerate(c):
for j,J in enumerate(c):
K[2*I,2*J] += Ke[2*i,2*j]
K[2*I+1,2*J] += Ke[2*i+1,2*j]
K[2*I+1,2*J+1] += Ke[2*i+1,2*j+1]
K[2*I,2*J+1] += Ke[2*i,2*j+1]
print('assign nodal forces and boundary conditions')
f = np.zeros((2*num_nodes))
for i in range(num_nodes):
if nodes[i,1] == 0.0:
K[2*i,:] = 0.0
K[2*i+1,:] = 0.0
K[2*i,2*i] = 1.0
K[2*i+1,2*i+1] = 1.0
if nodes[i,1] == mesh_ly:
x = nodes[i,0]
f[2*i+1] = 20.0
if x == 0.0 or x == mesh_lx:
f[2*i+1] *= 0.5
print('solving linear system')
u = np.linalg.solve(K, f)
print('max u=', max(u), 'mm') # Adding units in mm
print('plotting displacement')
ux = np.reshape(u[0::2], (mesh_ny,mesh_nx))
uy = np.reshape(u[1::2], (mesh_ny,mesh_nx))
xvec = []
yvec = []
res = []
for i in range(mesh_nx):
for j in range(mesh_ny):
xvec.append(i*mesh_hx + ux[j,i])
yvec.append(j*mesh_hy + uy[j,i])
res.append(uy[j,i])
t = plt.tricontourf(xvec, yvec, res, levels=14, cmap=plt.cm.jet)
plt.scatter(xvec, yvec, marker='o', c='b', s=2)
plt.suptitle('FEA beam', fontsize = 12) # Added a title
ax.set_xlabel('x') # Added x label
ax.set_ylabel('y') # Added y label
plt.grid()
plt.colorbar(t, label="Displacement u",) # Added a label
plt.axis('equal')
plt.show()
create mesh material model - plane strain create global stiffness matrix assign nodal forces and boundary conditions solving linear system max u= 6.745400245721763 mm plotting displacement
Next, let's try to make this into our beam:
# Copyright 2022 Jorgen Bergstrom
# See also code from http://compmech.lab.asu.edu/codes.php
# Additions by Nicole Bakker (in the comments)
import numpy as np
import math
from matplotlib import pyplot as plt
def shape(xi):
x,y = tuple(xi)
N = [(1.0-x)*(1.0-y), (1.0+x)*(1.0-y), (1.0+x)*(1.0+y), (1.0-x)*(1.0+y)]
return 0.25*np.array(N)
def gradshape(xi):
x,y = tuple(xi)
dN = [[-(1.0-y), (1.0-y), (1.0+y), -(1.0+y)],
[-(1.0-x), -(1.0+x), (1.0+x), (1.0-x)]]
return 0.25*np.array(dN)
fig = plt.figure() # Adding this for labels
ax = fig.add_subplot(111) # Adding this for labels
print('create mesh')
# input
mesh_ex = 55 # discretization / number of dots in horizontal direction
mesh_ey = 12 # discretization / number of dots in vertical direction
mesh_lx = 50.0 # lenght of beam
mesh_ly = 10.0 # height of beam
# derived
mesh_nx = mesh_ex + 1
mesh_ny = mesh_ey + 1
num_nodes = mesh_nx * mesh_ny
num_elements = mesh_ex * mesh_ey
mesh_hx = mesh_lx / mesh_ex
mesh_hy = mesh_ly / mesh_ey
nodes = []
for y in np.linspace(0.0, mesh_ly, mesh_ny):
for x in np.linspace(0.0, mesh_lx, mesh_nx):
nodes.append([x,y])
nodes = np.array(nodes)
conn = []
for j in range(mesh_ey):
for i in range(mesh_ex):
n0 = i + j*mesh_nx
conn.append([n0, n0 + 1, n0 + 1 + mesh_nx, n0 + mesh_nx])
print ('material model - plane strain')
E = 100.0
v = 0.48
C = E/(1.0+v)/(1.0-2.0*v) * np.array([[1.0-v, v, 0.0],
[ v, 1.0-v, 0.0],
[ 0.0, 0.0, 0.5-v]])
print('create global stiffness matrix')
K = np.zeros((2*num_nodes, 2*num_nodes))
q4 = [[x/math.sqrt(3.0),y/math.sqrt(3.0)] for y in [-1.0,1.0] for x in [-1.0,1.0]]
B = np.zeros((3,8))
for c in conn:
xIe = nodes[c,:]
Ke = np.zeros((8,8))
for q in q4:
dN = gradshape(q)
J = np.dot(dN, xIe).T
dN = np.dot(np.linalg.inv(J), dN)
B[0,0::2] = dN[0,:]
B[1,1::2] = dN[1,:]
B[2,0::2] = dN[1,:]
B[2,1::2] = dN[0,:]
Ke += np.dot(np.dot(B.T,C),B) * np.linalg.det(J)
for i,I in enumerate(c):
for j,J in enumerate(c):
K[2*I,2*J] += Ke[2*i,2*j]
K[2*I+1,2*J] += Ke[2*i+1,2*j]
K[2*I+1,2*J+1] += Ke[2*i+1,2*j+1]
K[2*I,2*J+1] += Ke[2*i,2*j+1]
print('assign nodal forces and boundary conditions') # need to change this
f = np.zeros((2*num_nodes))
for i in range(num_nodes):
if nodes[i,1] == 0.0:
K[2*i,:] = 0.0
K[2*i+1,:] = 0.0
K[2*i,2*i] = 1.0
K[2*i+1,2*i+1] = 1.0
if nodes[i,1] == mesh_ly:
x = nodes[i,0]
f[2*i+1] = 20.0
if x == 0.0 or x == mesh_lx:
f[2*i+1] *= 0.5
print('solving linear system')
u = np.linalg.solve(K, f)
print('max u=', max(u), 'mm') # Adding units in mm
print('plotting displacement')
ux = np.reshape(u[0::2], (mesh_ny,mesh_nx))
uy = np.reshape(u[1::2], (mesh_ny,mesh_nx))
xvec = []
yvec = []
res = []
for i in range(mesh_nx):
for j in range(mesh_ny):
xvec.append(i*mesh_hx + ux[j,i])
yvec.append(j*mesh_hy + uy[j,i])
res.append(uy[j,i])
t = plt.tricontourf(xvec, yvec, res, levels=14, cmap=plt.cm.jet)
plt.scatter(xvec, yvec, marker='o', c='b', s=2)
plt.suptitle('FEA beam', fontsize = 12) # Added a title
ax.set_xlabel('Length of beam (mm)') # Added x label
ax.set_ylabel('Height of beam (mm)') # Added y label
plt.grid()
plt.colorbar(t, label="Displacement u",) # Added a label
plt.axis('equal')
plt.show()
create mesh material model - plane strain create global stiffness matrix assign nodal forces and boundary conditions solving linear system max u= 2.747833854475706 mm plotting displacement
Another option for FEM packages: http://www.juliafem.org/