#006 - SciPy in Structural and Civil Engineering, Part 3/3: Structural Dynamics

Exploring Structural Dynamics with SciPy: An In-Depth Journey into Engineering Analysis

Welcome back to the final installment of our exploration into using Python’s SciPy library for Structural and Civil Engineering applications.

In Part 1 and Part 2, our exploration was anchored in the static realm of structural engineering, using SciPy to solve linear algebra problems and perform optimization tasks. Now, in this final segment, we pivot towards the vibrant world of structural dynamics, delving into the ever-changing and moving aspects of engineering structures.

Update: November 19th, 2023 - I’m experimenting with an audio version of the newsletter, in podcast form, so if you prefer that type of format, make sure to let me know.

---

To quickly recap, SciPy is a Python library that provides modules for scientific computing, offering tools for optimization, linear algebra, integration, interpolation, and statistics. It extends the NumPy library, which offers numerical array operations.

SciPy

Structural Dynamics

Many books are dedicated to the topic of dynamics so in this post, we will just scratch the surface and hopefully ignite some curiosity.

Structural dynamics is a complicated topic and has caused me some consternation in my career. I’m certain this will continue. In terms of resources to help you wrap your head around the subject and its broad spectrum of applications, my go-to choice has always been Dynamics of Structures by AK Chopra. If you have another suggested resource, please let me know.

Structural Dynamics and SciPy

SciPy is well-suited for tackling problems in structural dynamics due to its powerful computational modules. Particularly, the integrate module which allows for solving ordinary differential equations (ODEs) that describe the motion of structures under dynamic loads.

Meanwhile, SciPy’s linalg module provides tools for linear algebra operations that are essential when dealing with the mass, damping, and stiffness matrices in dynamic analysis. Both of these capabilities make SciPy a valuable asset for engineers looking to analyze and design structures exposed to dynamic forces.

Modelling Dynamic Systems

Before we can analyze a structure's response to dynamic loads, we need to create a mathematical model of the system. This often starts with the simplest form, a single degree of freedom (SDOF) system. An SDOF system can be modelled as a mass attached to a spring and damper, representing the structure's mass, stiffness, and damping properties, respectively.

Single degree of freedom (SDOF) system

In SciPy, we can use arrays and functions to represent these properties and solve the system's motion equations.

Arrays in NumPy and SciPy are grids of uniformly-typed values, indexed by nonnegative integers, facilitating efficient storage, manipulation, and vectorized operations on large datasets, such as representing system states in motion equations at each timestep.

For more complex structures, a multi-degree of freedom (MDOF) system is required. Each degree of freedom represents an independent mode of motion of the structure. SciPy allows us to expand our models to MDOF systems by using larger matrices and more complex functions to describe these interdependent motions. This level of modelling is crucial for analyzing the behaviour of real-world structures under dynamic conditions. And this is the key complexity that makes 3D dynamic calculations so challenging. When you’re dealing with SDOF systems, it’s much easier to conceptualize the processes. Once you move into MDOF systems, things get complicated quickly and it can be rather challenging to track down cause and effect in an efficient manner.

The Direct Stiffness Method

The direct stiffness method is a foundational tool in structural analysis, enabling engineers to assess system responses to loads by breaking down complex structures into discrete elements. Its strength lies in transforming the complex interplay of forces and displacements into a manageable set of linear equations, typically resolved using matrix algebra.

Dynamic analysis extends the direct stiffness method by incorporating the time-dependent aspect of loads and helps solve the motion equations (often differential equations) over time steps, which isn't typically a focus of the direct stiffness method. However, for certain dynamic problems, especially those considering small vibrations around an equilibrium, the stiffness matrix derived from the direct stiffness method is a crucial component of the mass-stiffness-damping system that characterizes dynamic behaviour.

SciPy offers functions under its linalg module that are highly optimized for solving large systems of linear equations, eigenvalue problems, and matrix operations, which are core to the direct stiffness method. Even though SciPy is not specifically tailored for structural analysis; it provides the numerical and algebraic tools that can be adapted for such purposes. The direct stiffness method, with its explicit focus on structural systems, complements SciPy's capabilities by offering a structural-engineering context to the mathematical solutions that SciPy excels at.

The direct stiffness method is powerful because it can be applied to a wide range of structure types and is scalable for large systems, making it the foundation for most structural analysis and design software.

Taking a first principles approach, the direct stiffness method is used to derive the system's stiffness matrix, which, along with the mass and damping matrices, would form the basis for writing the equation of motion. SciPy's linalg module could then be used to solve these equations, which are typically second-order ordinary differential equations in the context of structural dynamics.

For dynamic problems, the equation of motion for an undamped MDOF system can be expressed as:

\(\mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{K}\mathbf{u}(t) = \mathbf{F}(t)\)

and for a damped system…

\(\mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{C}\dot{\mathbf{u}}(t) + \mathbf{K}\mathbf{u}(t) = \mathbf{F}(t)\)

where:

\(\begin{align*} \mathbf{M} & : \text{Mass matrix} \\ \mathbf{C} & : \text{Damping matrix (only for damped system)} \\ \mathbf{K} & : \text{Stiffness matrix} \\ \mathbf{u}(t) & : \text{Displacement vector as a function of time} \\ \dot{\mathbf{u}}(t) & : \text{Velocity vector as a function of time} \\ \ddot{\mathbf{u}}(t) & : \text{Acceleration vector as a function of time} \\ \mathbf{F}(t) & : \text{External force vector as a function of time} \end{align*}\)

In this form, the equation can be numerically integrated over time using methods available in SciPy's integrate module.

Example: Pedestrian Footbridge Vibrations

In this example, let’s consider a pedestrian footbridge (as a SDOF) subjected to a rhythmic load due to the collective footsteps of people crossing it in unison. As joggers take to the bridge on their morning routines, the consistent thud of their steps imparts a fluctuating force onto the structure, setting it into a dynamic response. This force repeats with a regular cadence, akin to a metronome, providing a consistent push-and-pull at a frequency that could resonate with the bridge's own natural frequency.

We will examine the dynamic response of the bridge to such a harmonic force, simplified to have a constant amplitude and frequency. This type of loading is typical of organized events where participants maintain a synchronized pace, such as charity runs, marches or stadium crowds. Note that these harmonic loads could come from a variety of other sources (e.g. machine vibration, wind, seismicity, traffic).

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Define the system parameters
mass = 1.0  # Mass in kg
damping = 0.05  # Damping coefficient in kg/s
stiffness = 100.0  # Stiffness in N/m

# Define the external force function
def external_force(t):
    # A simple harmonic force with amplitude 10 N and frequency 1 Hz
    return 10 * np.sin(2 * np.pi * 1 * t)

# Define the system ODEs
def system_odes(t, y):
    # y[0] is displacement, y[1] is velocity
    displacement = y[0]
    velocity = y[1]
    
    # Calculate acceleration
    acceleration = (external_force(t) - damping * velocity - stiffness * displacement) / mass
    
    return [velocity, acceleration]

# Initial conditions: [initial_displacement, initial_velocity]
initial_conditions = [0.0, 0.0]

# Time span for the simulation
t_span = (0.0, 10.0)
t_eval = np.linspace(*t_span, 1000)  # Time points where the solution is computed

# Solve the ODE
solution = solve_ivp(system_odes, t_span, initial_conditions, t_eval=t_eval)

# Plot the results
plt.figure(figsize=(10, 4))
plt.plot(solution.t, solution.y[0], label='Displacement')
plt.plot(solution.t, solution.y[1], label='Velocity')
plt.title('SDOF System Response to Harmonic Force')
plt.xlabel('Time (s)')
plt.ylabel('Response')
plt.legend()
plt.grid(True)
plt.show()

The scipy.integrate.solve_ivp function is used to solve initial value problems for a system of ordinary differential equations (ODEs). Here's how the process works for the example provided:

1. The system_odes function defines this equation in terms of a first-order system by expressing it as two coupled first-order differential equations. The variable y is a list where y[0] represents the displacement (x), and y[1] represents the velocity. The acceleration, which is the second derivative of displacement (x), is calculated by rearranging the equation of motion as follows…

   \(\ddot{x} = \frac{F(t) - c \dot{x} - k x}{m} \)

2. External Force: external_force is a function that computes the harmonic external force applied to the system as a function of time. This force is what drives the system's dynamics.

3. Initial Conditions: The state of the system at the initial time is provided (initial_conditions), which includes the initial displacement (0 m) and initial velocity (0 m/s).

4. Time Span: The time span (t_span) over which the solution should be computed is defined, along with specific time points (t_eval) where the solution is desired.

5. Solving the ODE: The solve_ivp function takes the system ODEs, the time span, initial conditions, and other optional arguments to solve the ODEs. It uses numerical integration methods to solve these equations. The default method is 'RK45', which is a flexible and adaptive step-size Runge-Kutta method, suitable for a wide range of problems.

6. Results: solve_ivp returns an object that includes the time points (solution.t) and the solution at those points (solution.y). For this system, solution.y[0] contains the displacement of the mass over time, and solution.y[1] contains the velocity.

Dynamic Response of a Single-Degree-of-Freedom System Subjected to Harmonic Excitation: The plot illustrates the time-history of displacement (in meters) and velocity (in meters per second) of the mass, revealing the system's behaviour under a sinusoidal force over a 10-second interval.

SciPy's solver works by stepping through the time span, at each step calculating the state of the system using the provided ODE function. The solver adjusts its step size based on the required accuracy and the behaviour of the solution, which is especially useful for stiff problems or problems requiring high precision over long timescales.

In summary, scipy.integrate.solve_ivp is used to numerically simulate the time evolution of the system under the given initial conditions and external force, producing the dynamic response of the system.

You can expand on this by incorporating more complex loading, non-linear behaviour, or multiple degree-of-freedom systems. I hope to discuss larger systems in the future.

Commercial Software and SciPy

Commercial structural analysis software offer a range of computational methods for dynamic analysis and while I do rely on these packages, I am adamant that SciPy offers some distinct advantages.

1. Flexibility and Customization: SciPy allows for custom analysis and the development of algorithms tailored to your specific engineering problems, which can be more adaptable than the often rigid frameworks of commercial packages.

2. Integration with Python Ecosystem: Being a part of the Python ecosystem, SciPy can be seamlessly integrated with other libraries for data manipulation, visualization, and machine learning, creating a powerful toolchain for complex data analysis and engineering computations.

3. Open Source and Community Support: SciPy is open-source, meaning it is freely available without the licensing costs associated with many commercial packages. SciPy is accessible to engineers and students around the world, fostering a more inclusive and global engineering community. This is one of my favourite aspects. No licencing/support headaches.

4. Transparency and Auditability: SciPy’s open-source nature allows for full transparency in the algorithms and numerical methods used.

5. Education and Knowledge Sharing: SciPy promotes a deeper understanding of the underlying numerical methods since users may need to implement or modify algorithms themselves, encouraging a more profound learning experience and knowledge sharing within the engineering community. This is a blessing and a curse.

6. Automation: Utilizing SciPy within Python scripts enables the automation of recurring tasks, analytical processes, and report generation. This level of automation may not be readily achievable in commercial software that primarily operates through graphical user interfaces (GUIs). However, it's important to note that while SciPy itself doesn't come with a GUI, you have the option to develop one if needed.

7. Scalability: SciPy ‘s computational models can be scaled up with relative ease to handle larger data sets or more complex simulations without significant changes to the codebase. This is a big deal.

In conclusion, our series has examined SciPy‘s potential to supplement your structural engineering process. With a focus on real-world application, we've demonstrated that through tools like linear algebra, integration, optimization, and interpolation. SciPy empowers engineers to tackle complex, dynamic systems with a more autonomous level of efficiency, flexibility, and control.

As we continue to leverage Python's powerful libraries, the scope for creative solutions and process refinement in your engineering workflows continue to broaden.

If you are interested in learning more about Python for structural or civil engineering, come and sign up at flocode.dev. We are building big things.

Your feedback and experiences, especially any unique applications or insights related to this series, are highly valued. Sharing these can foster a richer learning environment for all.

See you in the next post!

James 🌊