PyJive workshop: Reproducing dynamics experiments

PyJive workshop: Reproducing dynamics experiments#

import numpy as np
import matplotlib.pyplot as plt

import contextlib
import os
from urllib.request import urlretrieve

import sys
pyjivepath = '../../../pyjive/'
sys.path.append(pyjivepath)

if not os.path.isfile(pyjivepath + 'utils/proputils.py'):
    print('\n\n**pyjive cannot be found, adapt "pyjivepath" above or move notebook to appropriate folder**\n\n')
    raise Exception('pyjive not found')

import main
from utils import proputils as pu
from names import GlobNames as gn

%load_ext autoreload
%autoreload 2

%matplotlib widget

**pyjive cannot be found, adapt "pyjivepath" above or move notebook to appropriate folder**

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
Cell In[1], line 14
     12 if not os.path.isfile(pyjivepath + 'utils/proputils.py'):
     13     print('\n\n**pyjive cannot be found, adapt "pyjivepath" above or move notebook to appropriate folder**\n\n')
---> 14     raise Exception('pyjive not found')
     16 import main
     17 from utils import proputils as pu

Exception: pyjive not found

In this exercise we will look into using pyJive to reproduce a hammer test performed on a rail beam. The goal of the experiment is to determine natural frequencies, mode shapes and time-domain acceleration response of a rail beam hit with an instrumented hammer. Accelerometers are placed along the beam and the location where the hammer hits can also be changed.

To make matters simple, we consider the following setup:

mesh

where we hit the hammer on one end of the beam and measure acceleration at midspan.

We consider the following material properties:

\[ E = 210\,\text{GPa}\quad\nu=0.3\quad\rho=7850\,\text{kg/m}^3 \]

and the following section properties:

\[ A=69.77\,\text{cm}^2\quad I = 2337.9\,\text{cm}^4 \]

and the span of the beam is:

\[ \ell = 3.259\,\text{m} \]

In this notebook the goal is to try to reproduce natural frequencies and the acceleration signal at midspan when we hit the hammer with a specific force signal. The experimental results we use here were obtained at the Macrolab of the Civil Engineering and Geosciences faculty of TU Delft.

Part 1: Natural frequencies#

First we will look at the natural frequencies. The experimental results for the first three frequencies are:

\[ \omega_1^\mathrm{experiment} = 95.4\,\text{Hz}\quad \omega_2^\mathrm{experiment} = 246.6\,\text{Hz}\quad \omega_3^\mathrm{experiment} = 529.1\,\text{Hz} \]

and we can also look at values coming from a simple 1D analytical solution, for reference:

\[ \omega_1^\mathrm{analytical} = 100.4\,\text{Hz}\quad \omega_2^\mathrm{analytical} = 276.7\,\text{Hz}\quad \omega_3^\mathrm{analytical} = 542.4\,\text{Hz} \]

We now try to reproduce these result with a 1D model in pyJive:

Task 1: Compute and compare natural frequencies

  • Set up rail_modal.pro and rail.geom with a FrameModel and a ModeShapeModule. Use an existing property file from a previous workshop as starting point, think about which of the workshops has a file that would provide the best starting point for this simulation

  • Adjust material and section properties to the values above. TIP: Make sure the units match, for instance by converting all properties to \(\text{Pa}\) and \(\text{m}\)

  • Try to reproduce the free-free boundary conditions of the experiment. How can that be achieved and what is the consequence of the BC choice when computing natural frequencies?

  • Discretize your mesh gradually. How many frame elements do we need for this problem?

  • In the block below we already provide a print call for the first 6 eigenfrequencies. We divide them by \(2\pi\) to convert them from \(\text{rad}/\text{s}\) to \(\text{Hz}\)

props = pu.parse_file('rail_modal.pro')

globdat = main.jive(props)

print(globdat[gn.EIGENFREQS][0:6]/2/np.pi)

Part 2: Implicit analysis#

In the experiment, the natural frequencies above are obtained by hitting the beam with a hammer at a number of locations and measuring the resulting acceleration at several other points of the beam. The hammer is instrumented with a load cell, from which it is possible to extract a load signal in \(\text{N}\). From the accelerometer measurements we can obtain acceleration signals in \(\text{m}/\text{s}^2\).

We leave the process to go from force and acceleration measurements to natural frequencies out of scope for this notebook. Here we will just attempt to reproduce one of the experiments:

  • Hit the right edge of the beam with the hammer

  • Measure accelerations at midspan

In the real experiment the accelerations gradually go back to zero due to damping. We can try to reproduce this here with numerical damping.

First we download the load and acceleration signals from the experiment and take a look at them:

# download input files (if necessary)
def findfile(fname):
    url = "https://gitlab.tudelft.nl/cm/public/drive/-/raw/main/hammer-test/" + fname + "?inline=false"
    if not os.path.isfile(fname):
        print(f"Downloading {fname}...")
        urlretrieve(url, fname)

findfile("load_signal.dat")
findfile("accel_signal.dat")

# read acceleration signal
accel = np.loadtxt('accel_signal.dat')

# store time increment and number of steps
dt = accel[1,0]-accel[0,0]
maxstep = accel.shape[0]

# plot the load signal, for reference
load = np.loadtxt('load_signal.dat')
fig, axs = plt.subplots(1,2,figsize=(6,4),layout='constrained')
axs[0].plot(accel[:,0],load)
axs[0].set_xlabel(r'time [s]')
axs[0].set_ylabel(r'load [N]')
axs[0].set_title('Full load signal')
axs[1].plot(accel[0:40,0],load[0:40])
axs[1].set_xlabel(r'time [s]')
axs[1].set_ylabel(r'load [N]')
axs[1].set_title('Zoomed-in view')
plt.show()

# plot the acceleration signal, for reference
plt.figure()
plt.plot(accel[:,0],accel[:,1])
plt.xlabel(r'time [s]')
plt.ylabel(r'acceleration [m/s$^2$]')
plt.title(r'Acceleration at midspan')
plt.show()

We will therefore be looking at an interval of about \(0.8\,\text{s}\) and model it with \(20000\) time steps. Now we move on to the actual task.

Task 2: Compute an acceleration time signal

  • Set up rail_implicit.pro. You will now need a node group at midspan. You can create that directly in the properties file by setting xtype = mid when specifying the node group, or by dividing the beam into two members in rail.geom and explicitly defining a node at midspan

  • Set up a NeumannModel at the right edge of the beam. Make sure to set values = [-1.0] and timeSignal = load_signal.dat to have pyJive read the load signal directly from the file

  • Use NewmarkModule to run the simulation implicitly. Use the values for dt and maxstep we stored above

  • Start with \(\gamma=0.5\) and \(\beta=0.25\) and introduce damping by increasing \(\gamma\) to try to reproduce the experimental signal. Make sure \(\beta\) is also high enough.

  • To plot the accelerations, include an AccelerationModule in your properties file, and set groups to the node group you created for the node at midspan

  • Reflect on what would happen if we were to use the numerical acceleration signal to compute the frequencies using FFT. Would you expect good results? What would be the influence of the numerical damping?

props = pu.parse_file('rail_implicit.pro')

globdat = main.jive(props)

Part 3: Mode shapes#

For this part we do not have experimental results, but we can compare mode shapes obtained with the numerical model with those coming from an analytical solution. For a free-free beam, the first three mode shapes can be computed from:

\[ a_n(x) = c_n\left[\sin\beta_nx+\sinh\beta_nx + \alpha_n\left(\cos\beta_nx+\cosh\beta_nx\right)\right] \]

where \(c_n\) is an arbitrary scaling factor,

\[ \alpha_n = \left(\displaystyle\frac{\sin\beta_n\ell-\sinh\beta_n\ell}{\cosh\beta_n\ell-\cos\beta_n\ell}\right) \]

and:

\[ \beta_1\ell=4.730041\quad\beta_2\ell=7.853205\quad\beta_3\ell=10.995608 \]

We can then compare these shapes with the ones we get with FEM:

Task 3: Compare mode shapes with analytical solution

  • Go back to rail_modal.pro and run the simulation again

  • Extract the modes from globdat using the function provided below and plot them

  • Compute the analytical solution for the modes and plot those together with the numerical predictions. Scale the analytical solution accordingly. TIP: a simple estimate for \(c_n\) can be obtained by aligning both solutions at \(x=0\)

  • Discretize your mesh gradually. What do you observe?

props = pu.parse_file('rail_modal.pro')

globdat = main.jive(props)

def get_mode_shape (mode, globdat):
    nset = globdat[gn.NSET]
    ds   = globdat[gn.DOFSPACE]    
    x = []
    dy = []
    
    for n in range(len(nset)):
        coords = nset[n].get_coords()
        x.append(nset[n].get_coords()[0])
        dy.append(globdat[gn.MODALSHAPES][mode][ds.get_dof(n,'dy')])

    order = np.argsort(x)

    return np.array(x)[order], np.array(dy)[order]
Contents