Tutorial 3.1: Solving the dynamic motion of a rod using Finite Differences#
In this tutorial you will learn how to solve the dynamic motion of a rod using Finite Differences. Here we will follow the example presented in the lectures.
Let’s consider the EOM of a rod of length \(L\):
with \(\rho\) the density, \(A\) the cross section area, \(E\) the young modulus and \(q\) an external (distributed) load.
The system is subject to boundary conditions:
\(u=0\) at \(x=0\)
\(EAu{'}=F(t)\) at \(x=L\)
The roadmap to finde the beam deflection using Finite Differences#
The steps needed to discretize the rod are listed next.
Discretize space into \(N+1\) points.
Discretize the EOM of the rod. To do this, you will have to approximate the 2nd order derivative with FD. Use a second order approximation. You can find the coefficients here (or google `finite difference coefficients’ and click the Wikipedia page).
Apply the boundary conditions. Use the definition of discrete derivatives and EOM to find the relation between the ghost points and the interior points.
Combine all equations into matrices.
Solve the resulting ODE system
While this may seem as a rather abstract problem, it is one encountered in a lot of civil engineering appilications. Take for example the pile driving process of an offshore wind turbine monopile, such as this Heerema hammering example
Let’s first install all required packages
# Load packages
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
import scipy as sp
from scipy.integrate import solve_ivp
Step 1: Discretize the domain#
L = 10
N = 10
x = np.linspace(0,L,N+1)
print(x)
[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]
Step 2: Discretize the EOM of the rod#
Using the coefficients for a centered FD scheme for 2nd order derivatives with 2nd order accuracy we have:
And for the 1st order derivative:
Replacing these expressions into the equation of motion we get the following discrete system:
Step 3: Apply boundary conditions#
At \(x=0\) we have that \(u(0)=0\), then
Since we know the values of \(u(0)\) for all times, we will not include them on the system. Using all this relation, we get the discrete equation for \(i=1\) as:
On the other side, we use the force relation
Thus
So, the discrete equations of motion for \(i=N\) is:
Step 4: Matrix form#
Summarizing we have the following discrete (in space) equations:
For \(i=1\):
For \(i=2,...,N-1\):
For \(i=N\):
This is equivalent to the following system:
And in a compact form:
Now we can define the matrices:
# Construct the matrices
M = np.eye(N)
K = np.zeros((N,N))
## Define some parameters (no physical sense)
E = 2e7
A = np.pi*2*0.07
rho = 4000
## Define the force
def q(x,t):
return 1e2*x*np.sin(3.0*t)
def F(t):
return 1e4*np.sin(0.3*t)
# Fill in stiffness matrix
K[0,0] = -2; K[0,1] = 1
K[N-1,N-2] = 2; K[N-1,N-1]=-2;
for i in range(1,N-1):
K[i,range(i-1,i+2)] = [1, -2, 1]
# Apply scaling
dx = L/N
M = rho*A*M
K = E*A/(dx**2)*K
# Fill in force
def Fvec(t):
Ft = np.zeros((N,))
for i in range(0,N):
xi = x[i+1]
Ft[i] = q(xi,t)
Ft[N-1] = Ft[N-1] + 2*F(t)/(dx)
return Ft
Step 5: Solve the ODE system#
Now, we just need to apply what we learned in the previous modules to solve the dynamic problem. The main novelty is that here we work with matrices.
Keen studens may notice the here specified solver uses the Radau method. The Radau method is one of many options in the list of Runge-Kutta integration methods, like for example Forward Euler.
While internal variations of the Radau method exist, the key take-away points are the full implicit nature and high order of accuracy, but at a large computational cost. The implicit nature makes the method well suited for very stiff problems. These problems normally require small time steps to capture the, due to the stiffness, fast accelerations and high-frequency effects. The use of an implicit method generally allows for a greater time step, then reducing the computational cost.
# Define the ODE function (given q, return dq)
def q_dot(t,q):
displ = q[0:N]
velos = q[N:2*N]
accel = np.linalg.inv(M).dot((Fvec(t)+K.dot(displ)+0.1*K.dot(velos))) #0.1 is damping factor, can change
dq = np.zeros((2*N,))
dq[0:N] = velos
dq[N:2*N] = accel
return dq
# Define initial state
q0 = np.zeros(2*N)
# Define time interval and time evaluation points
t0 = 0
tf = 100.0
nt = 50
tspan = np.linspace(t0,tf,nt+1)
# Solve (we use Radau method because the ODE is very stiff and it works best for this case)
sol = solve_ivp(fun=q_dot,t_span=[t0,tf],y0=q0,t_eval=tspan,method='Radau')
# Plot in local mode -- creates animated plot
import time
import pylab as pl
from IPython import display
for i in range(0,len(tspan)):
#fig,ax = plt.subplots()
display.clear_output(wait=True)
#display.display(pl.gcf())
plt.plot(np.linspace(0,L,N),sol.y[0:N,i])
plt.ylim((-0.015,0.015))
plt.xlabel("Beam position [m]")
plt.ylabel("Displacement [m]")
plt.title("Beam displacement animation")
plt.savefig(f".././images/Module3/w3_t1_gif/shot{i+1}.png")
time.sleep(0.05)
plt.show()
Now displaying the animated plots, gif conversion via https://ezgif.com/jpg-to-gif.
# Pick time value:
time_pick = 24
time_closest = sol.t.flat[np.abs(sol.t - time_pick).argmin()]
time_idx = np.where(sol.t == time_closest)[0]
disp_pick = np.reshape(sol.y[0:N,time_idx],N)
disp_full = np.concatenate(([0], disp_pick))
print(disp_full)
plt.plot(x,disp_full)
plt.xlabel("Beam position [m]")
plt.ylabel("Displacement [m]")
plt.title(f"Beam displacament at {time_closest:.0f} s");
[0. 0.00119706 0.00238766 0.00356537 0.00472378 0.00585652
0.00695731 0.00801991 0.0090382 0.01000615 0.01091787]
Exercise: Add friction to the system#
Solve the same problem, but considering friction in the rod for \(x<L/2\). In that case, the equation of motion reads:
with
Compare the solution using \(k_0=1.0e8\) and \(k_0=1.0e5\).