top of page

Dynamics Modeling/Simulation with Python and SymPy, Alternative to MatLab

Throughout my career as a mechanical engineer, I have always emphasized my competency as a programmer. These skills have lent themselves well to my specialization in dynamics, which has always been my favorite engineering subject - one where I had excellent teachers and mentors. Through these mentors, I learned a workflow combining hand derivations, Mathematica to perform symbolic computations and check my work, and MATLAB to implement numerical methods and simulation. I’ve applied this workflow alongside commercial finite element and multibody dynamics software to diverse problems, including shock and vibration load propagation, stability, and aeroelasticity analysis, experimental system identification and modal analysis, etc.


Souvenir Portugal
Portugal


Souvenir Morocco
Morocco


Souvenir Colombia
Colombia

A little over four years ago, I took a sabbatical from my job to take six months to travel Europe, North Africa, and South America. I wanted to keep my engineering skills sharp during my time away from work, and set a goal to learn Python. As a project oriented, “learn-by-doing” type, I looked for a way to extend my existing skill set in mechanism dynamics by applying that skill in what was at the time a brand new computer language for me. I learned that NumPy and SciPy would cover many of my MATLAB use cases, and similarly SymPy for Mathematica. Furthermore, each of these Python tools can be used free of charge, and may even be launched in a browser with Google Colaboratory without requiring any installation.



Worquick: A Python Notebook with Sympy Mechanics in Google Colaboratory
A Python Notebook with Sympy Mechanics in Google Colaboratory

The Python module that I have found to have the most utility to my particular line of work, and which I will be writing about here, is the sympy.physics.mechanics module. This module standardizes many of the processes that have challenged me over years in the derivation of equations of motion (EOM), such as obtaining the time derivatives to calculate velocity and acceleration in a mechanism with complex kinematics or applying Newton-Euler, Lagrange, or Kane’s methods. I find the module to have utility toward many classical dynamics problems, like pendulums, linkages, three-dimensional (3D) rigid body mechanics, etc. In my real job I’ve successfully applied it to calculating actuation loads in moving mechanical assemblies, and predicting force or acceleration measurements that are recorded by test equipment. As I’ve delivered models and analysis using these methods, I’ve found that this is a powerful tool while simultaneously being fairly obscure. My goal with this post is to try to pass along the experience and knowledge I have gained for others to apply to their own career.



Workflow


Aside from the code itself, I have found that mastering sympy.physics.mechanics has led to a standardization of my own workflow for analyzing a dynamics problem. The Python portion of this workflow may be summarized by the eight steps I have listed below, which I will elaborate on in the example to follow.


  1. Import the required modules

  2. Define constants and variables

  3. Define reference frames

  4. Define points

  5. Create bodies

  6. Create loads

  7. Derive EOM

  8. Simulate

There are, of course, additional steps in the dynamics problem solving process that exist outside of the Python workflow, or require more advanced or unique steps that are not so easily standardized. An example of the former is the generation of a diagram as a starting point to the EOM generation that effectively communicates the design intent, critical dimensions, or other knowledge necessary to starting the model. More advanced or unique steps may include implementation of 3D mechanics, systems with high dimensionalities, or systems with constraints. In these cases the standard transformations may need to be augmented with expressions that are obtained by hand or from a textbook, or in other modules not contained in the Sympy Mechanics module. I intend to cover these in a future post, whereas here I will focus on the basic workflow and Python syntax for a simple example.



Example


A simple, single degree-of-freedom pendulum is selected as an example to demonstrate the workflow as mentioned above. A Google Colaboratory notebook may be found at the link below containing detailed, line by line descriptions of the Python code, with links to the official documentation. 


For this article, we will stick to the highlights, however, you may view and experiment with the source notebook as you wish.


I have also created a companion video, which is hosted on the Mobile Multibody Dynamics YouTube page.


Simple Pendulum Simulations in MOMDYN for Varying Initial Condition 


This video, and the simulation plots in the notebook and subsequent sections of this post, are intended to show that even with something as simple as a single pendulum, you will see significant variations in response, often with very minimal changes in model properties, or in this case, initial conditions.


The derivation is based on the diagram as depicted below. Notably, the model is for a rigid pendulum link with length L (1 meter in the example), mass m, and moment of inertia J, with the latter calculated assuming a uniform rectangular beam of aluminum with 10 centimeter width and 1 centimeter thickness. As the beam is uniform, the center of mass is placed at half of the link length. The angle of the link relative to vertical is defined as θ, and the angular rate as ω. A uniform gravitational acceleration g is applied in the downward direction.



Worquick: Pendulum Illustration
Pendulum Illustration

As previously noted, the complete derivation including comments and descriptions of what each line item is doing is fairly extensive, however, when stripped down to strictly the sympy sections, the workflow defined earlier is as follows:



A few notable symbolic expressions that are derived in this process are the velocity and acceleration at the center of the link,


Speed equation snapshot

Acceleration equation snapshot

and derivation of kinetic energy


Derivation of Kinetic energy equation snapshot

which is used to determine a critical angular rate


critical angular rate equation snapshot

above which the pendulum will continue to rotate past the fully upright position, rather than falling back and oscillating. The end result of this derivation is the EOM,


equation of mouvement snapshot

which is reorganized into a first order ordinary differential equation (ODE),


ordinary differential equation (ODE) snapshot

which is subsequently simulated using an ODE solver, in our case the solve_ivp function from scipy. 

An initial example simulation is run with an initial angular rate of 1 rad/sec. This is a fairly benign initial rate, with the result being nearly sinusoidal motion with an amplitude of about 0.25 rad (about 15 deg).



Diagram Angle/Time
Diagram Angle/Time

A more substantial analysis is performed by varying the initial rate from 1 to 10 rad/sec, in very fine increments. The angle and rate plots in insets (a) and (b) below show the plot traces for each of these simulations, with the initial rate encoded in the colors ranging from blue (smallest) to red (largest). As the initial rates exceed the critical rate, which was determined to be 7.65 rad/sec in the notebook, the simulated results begin to show complete rotations around the vertical, instead of the oscillating behavior of the initial example. Inset (c) shows a phase plot, with the low initial rate cases plotting circles in the angle/rate domain, and the high initial rate cases veering off to the right. Frequency spectra for the angular rate are displayed in inset (b), where for low initial rates the oscillation rate is about 0.6 Hz. As the initial rate increases, but remains below the critical rate, the oscillation frequency decreases as the link will tend to pause at its angular extremes when it is above the horizontal. Above the critical rate, the observed frequency will again increase, as the spectral peak will track with the amount of time it takes for the link to complete one complete rotational cycle, rather than the natural frequency of oscillation.



Summary


The sympy.physics.mechanics package provides an excellent and free method to perform symbolic math for the derivation of EOM in mechanism dynamics analysis. This post discussed a standardized workflow including eight steps, with a link to a comprehensive notebook with line by line descriptions of the Python code. The author has applied this approach to mechanisms with low-to-medium complexity, including for linkage and gear trains used in actuation systems, and as a supplement to analyses that are performed in finite element software. A future post is planned to tackle more advanced usage of these tools, including analysis of 3D mechanisms, mechanisms with constraints, and automatic code generation.


 

Marcin Szyc Profile Picture Worquick

Eliott Radcliffe spent the first 15 years of his career working in aerospace, first with a fellowship at NASA where he worked in aeroacoustics research, then at the Johns Hopkins Applied Physics Laboratory. At the latter, he specialized in structural dynamics and programming in MATLAB, and held multiple formal leadership roles.

His career goal is to start a software company developing educational engineering applications for iOS and Android. He develops code in native Swift (iOS), and Kotlin (Android) programming languages, alongside Python and C#. Mobile Multibody Dynamics is available on the App Store and Google Play.

More information here : www.dc-engineer.com

 

Kommentare


bottom of page