PyJive workshop: Time-dependent analysis

PyJive workshop: Time-dependent analysis#

In this notebook, usage of the ExplicitTimeModule and NewmarkModule for time-dependent analysis of structures is explored. These modules are for solving dynamics problems in the time domain. The first of these is for explicit analysis with the central difference scheme, the second for implicit analysis with Newmark time integration.

In the notebook, we will also use the AccelerationModule for recording accelerations and velocities as additional output to the analysis. There is also an NLNewmarkModule in pyjive for implicit analysis of nonlinear dynamics problems, but we will not use it here.

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

# define a function we will use for plotting force and displacement signals
def plotForceDisp(globdat):
    plt.figure(figsize=(8,3))
    plt.subplot(1,2,1)
    plt.plot(globdat['lodi']['left']['load']['dx'])
    plt.plot(globdat['lodi']['right']['load']['dx'])
    plt.legend(['left','right'])
    plt.ylabel('F')
    plt.xlabel('time step')
    plt.subplot(1,2,2)
    plt.plot(globdat['lodi']['left']['disp']['dx'])
    plt.plot(globdat['lodi']['right']['disp']['dx'])
    plt.legend(['left','right'])
    plt.ylabel('u')
    plt.xlabel('time step')
    plt.show()

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

findfile("bar_explicit.pro")
findfile("bar_implicit.pro")
findfile("bar_harmonic.pro")
findfile("bar.geom")

Analysis 1: Wave propagation in a bar, explicit#

The first case concerns wave propagation in a bar. It is a bar problem, we have three models that can solve bar problems, the BarModel, the SolidModel, and the FrameModel. We use the last of these for its dedicated postprocessing routines, specifically its option to plot the normal force distribution. We first run the case with the explicit time module.

Task 1: Explicit analysis

Look at the code

  • The ExplicitTimeModule applies central difference scheme. Can you spot a difference with the equation from the lecture where we had \(\hat{\mathbf{M}}\mathbf{a}_{n+1} = \hat{\mathbf{f}}_n\quad\) with \(\quad\hat{\mathbf{M}} = \frac{1}{\Delta t^2}\mathbf{M} + \frac{1}{2\Delta t}\mathbf{C}\quad\) and \(\quad \hat{\mathbf{f}}_n =\frac{1}{\Delta t^2}\mathbf{M}\left(2\mathbf{a}_n-\mathbf{a}_{n-1}\right) + \frac{1}{2\Delta t}\mathbf{C}\mathbf{a}_{n-1} - \mathbf{K}\mathbf{a}_n+\mathbf{f}_n\)? How does this difference affect the functionality?

Look at the input files

  • What are the geometry and boundary conditions?

Run the cell and inspect the results

  • There are no time-dependent boundary conditions, yet this is a dynamics problem. What makes the problem time-dependent?

Investigate stability

  • Change the time step size and rerun. What is the critical time step size for stability?

props = pu.parse_file('bar_explicit.pro')
globdat = main.jive(props)
plotForceDisp(globdat)

Analysis 2: Wave propagation in a bar, implicit#

Another input file is provided for implicit analysis with the NewmarkModule.

Task 2: Implicit analysis

  • Run the cell and inspect the results

  • What happens if you increase the time step size?

  • What happens if you add numerical damping (increase \(\gamma\) adapt \(\beta\) to maintain \(\beta=\gamma/2\))?

  • The behavior in the first time steps is fundamentally different from that in explicit dynamics. What is the numerical source for this?

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

# run code with output suppressed
print('running pyjive with output suppressed...')
with contextlib.redirect_stdout(open(os.devnull, "w")):
    globdat = main.jive(props)
    
plotForceDisp(globdat)

Analysis 3: Harmonically loaded bar#

Now, instead of a constant load we apply an harmonic load.

Task 3: Harmonic loading

Run the first cell below.

  • Check how time-dependent boundary conditions are applied

  • Inspect the resulting velocity and accelerations

Run the second cell below.

  • What is changed in the case definition?

  • Would it be possible to obtain this result with the ExplicitTimeModule?

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

print('running pyjive with output suppressed...')
with contextlib.redirect_stdout(open(os.devnull, "w")):
    globdat = main.jive(props)
plotForceDisp(globdat)

# plot acceleration and velocity at the loaded end
plt.figure()
plt.plot(globdat['acc']['right']['disp']['dx'])
plt.plot(globdat['acc']['right']['velo']['dx'])
plt.plot(globdat['acc']['right']['accel']['dx'])
plt.legend(['disp','velo','accel'])
plt.ylabel('u/v/a')
plt.xlabel('time step')
plt.show()

props = pu.parse_file('bar_harmonic.pro')
props['model']['neum']['deltaTime'] = 1;
props['stepper']['deltaTime'] = 1;
props['model']['neum']['timeSignal']  = 'np.sin(0.01*t)**2'

print('running pyjive with output suppressed...')
with contextlib.redirect_stdout(open(os.devnull, "w")):
    globdat = main.jive(props)
    
plotForceDisp(globdat)
Contents