{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Elementary Examples in Inverse Dynamics\n", "\n", "Falk Mielke, FunMorph, Universiteit Antwerpen\n", "\n", "You can download this notebook [here](user://notebooks/ID_LMX0b_FForces.ipynb?download&target=_blank)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is one part of a series of blog posts on Inverse Dynamics.\n", "- [Prologue A: Wrench Town Rock](../id_lmx0a_wrenches)\n", "- [Prologue B: The Fictitious Force Awakens](../id_lmx0b_fforces)\n", "- [Elementary Example 1: Stuttering Motor (Euler Force)](../id_lmx1_stutter)\n", "- [Elementary Example 2: Needle On A Vinyl (Centrifugal Force)](../id_lmx2_vinyl)\n", "- [Elementary Example 3: Radial Slider (Coriolis Force)](../id_lmx3_slider)\n", "- [Application: N-Link Inverse Dynamics](../id_lmx4_nlink)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as NP # numerics\n", "import sympy as SYM # symbolic operations\n", "import sympy.physics.mechanics as MECH # some physics/mechanics tools\n", "import matplotlib.pyplot as MPP # plotting\n", "\n", "import WrenchToolbox as WT # contains an improved \"Wrench\" object\n", "\n", "# printing of formulas\n", "SYM.init_printing(use_latex='mathjax', pretty_print = False)\n", "MECH.init_vprinting(use_latex='mathjax', pretty_print = False)\n", "\n", "# shorthand to get the basis of a reference frame\n", "GetCoordinates = lambda rf: [rf[idx] for idx in rf.indices]\n", "\n", "# transforming a sympy vector to a matrix and back\n", "MatrixToVector = lambda mat, rf: sum([elm * base for elm, base in zip(mat, [rf.x, rf.y, rf.z]) ])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The custom Wrench Toolbox [can be downloaded here.](http://mielke-bio.info/share/WrenchToolbox.py?download&target=_blank)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The Fictitious Force Awakens\n", "\n", "In the previous notebook, I introduced my notation to solve inverse dynamics for systems of rigid bodies. \n", "\n", "Based on a simple example, I demonstrated the general procedure, calculated the equations involved in Balance of Forces and Balance of Moments, and introduced Wrenches to naturally combine both balances. \n", "These balances directly emerge from the **conservation of linear and angular momentum**, so they can be regarded as a physical first principle." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I closed the example by pointing at a mismatch between the expected Coriolis Force (zero) and the calculation according to a naïve application of the Wikipedia formula. \n", "I will continue this here by calculating Balances in the rotating reference frame. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The main purpose of this notebook series is to see how **fictitious forces** enter the balance equations in inverse dynamics.\n", "My idea is to use the information about the inertial properties of the system under study to dissect the net joint force and moment into (i) biological components and (ii) components that stem from passive, fictitious forces. \n", "I want to know whether my piglets are on the dark side (greedily exploiting fictitious forces to reduce energy demand of their movements) or on the light side (resisting the fictitious force). \n", "\n", "There are three fictitious forces: the **Euler** force, the **Centrifugal** force, and the **Coriolis** force.\n", "Starting in the next notebook with the Euler force, I intend to construct elementary examples in which one can see the fictitious forces in isolation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fictitious forces only emerge when executing the balance equations in a non-static reference frame. \n", "This is the sole reason that makes them \"fictitious\", because reference frames are, and one might argue that the static frame is the only one that is relevant. \n", "The consequences of fictitious forces are real, though; they can be felt and will cause trouble, not only when you try to walk drunk over a merry-go-round.\n", "\n", "The balance equations I previously demonstrated were all in the static reference frame $\\{S\\}$. So to get to the fictitious forces, we'll have to repeat the procedure in the \"body\" frame $\\{B\\}$. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following line will reload most of the [preparation](WrenchExampleReload.py?download&target=_blank) so that we do not have to repeat it." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# %load WrenchExampleReload.py\n", "%run WrenchExampleReload.py\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall the problem:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This time, we'll move into the body frame $\\{B\\}$. \n", "But to be able to do that, we'll revisit and extend the kinematics." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Kinematics Ambiguity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll use shortcuts for selecting which derivative of position vectors applies. \n", "- p is the position\n", "- v is the velocity (first time derivative of position)\n", "- a is the acceleration (second time derivative of position)\n", "\n", "But these are only indices to lists that are assembled below." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# shorthands for the derivative levels\n", "p = 0 # zero'th derivative: the position vector\n", "v = 1 # first derivative: the velocity\n", "a = 2 # second derivative: the acceleration" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear Kinematics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The situation is more complex than before, where we only had to use the three derivatives for a single $x$ and $\\omega$. \n", "Instead of constructing longer and longer variable names for many variables (e.g. eightteen of the kind x_S_joint_ddot), I'll give this some structure by using multi-level dictionaries and lists." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# prepare an empty dictionary\n", "x = {}\n", "\n", "# start with the position vectors\n", "x['S'] = {pt: [pointmass.pos_from(refpoint).express(world).to_matrix(world)] \\\n", " for pt, refpoint in [('O', origin), ('J', joint), ('P', pointmass)] \\\n", " }\n", "x['B'] = {pt: [pointmass.pos_from(refpoint).express(body).to_matrix(body)] \\\n", " for pt, refpoint in [('O', origin), ('J', joint), ('P', pointmass)] \\\n", " }\n", "\n", "# automatically get time derivatives, in loops\n", "# loop reference frames\n", "for frame in x.keys():\n", " # loop reference points\n", " for refpoint in x[frame].keys():\n", " # loop differentials\n", " for diff_nr in range(2): \n", " x[frame][refpoint].append(SYM.simplify(x[frame][refpoint][-1].diff(t)))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "one can access them easily, for example by writing x['B']['O'][p], which gives the *point mass (x) position vector (p) relative to the origin (O) in the body frame (B)*." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 1.0 0]\n", "[0 1.0 0]\n" ] } ], "source": [ "# printout of position vectors in different frames\n", "SYM.pprint(PlugAllIn(x['B']['O'][p]).T)\n", "SYM.pprint(PlugAllIn(trafo['B']['S']*x['S']['O'][p]).T)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the position vactors are equal in both reference frames, except for the mandatory rotation (transformation). \n", "But look at their derivatives:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 0 0]\n", "[-1.0 0 0]\n" ] } ], "source": [ "# in contrast, the time derivatives depend on reference frame\n", "SYM.pprint(PlugAllIn(x['B']['O'][v]).T)\n", "SYM.pprint(PlugAllIn(trafo['B']['S']*x['S']['O'][v]).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is a first interesting observation ****: \n", "- The position vectors are identical in both reference frames.\n", "- However, the velocities are not identical. The derivatives in the body frame are blind to the change of the body frame position/orientation over time. \n", "- Of course this is equally true for the accelerations.\n", "\n", "Hence, there is an **ambiguity** in how to get kinematics vectors in the body frame: one can use zero vectors which are correct \"per definition\", or one could take the equally correct vectors from the static frame and transform them. \n", "\n", "The latter one, for distinction, will be denoted by a *tilde*:\n", "\n", "$\\tilde{x}_{B} = R_{BS}x_{S}$\n", "\n", "\n", "with $R_{BS}$ being the transformation matrix. \n", "The time derivatives of $\\tilde{x}_{B}$ are $\\dot{\\tilde{x}}_{B}$ and $\\ddot{\\tilde{x}}_{B}$, accordingly." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Angular Kinematics\n", "\n", "To a certain extent, this is analogous to linear kinematics. \n", "However, angular kinematics do not require a reference point. \n", "If the frame turns, it turns, and all points fixed to the frame turn the same way (i.e. if they don't move relative to each other). " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 0 0]\n", "[0 0 0]\n" ] } ], "source": [ "# the same dictionary magic as above\n", "ω = {}\n", "\n", "# angular positions\n", "ω['S'] = [SYM.Matrix([,,[φ]])]\n", "ω['B'] = [SYM.Matrix([,,])]\n", "\n", "# time derivatives\n", "for frame in ω.keys():\n", " for diff_nr in range(2): \n", " ω[frame].append(SYM.simplify(ω[frame][-1].diff(t)))\n", "\n", "# quick check: which of these matches the settings we entered for our reference frames:\n", "SYM.pprint(PlugAllIn(ω['S'][v] - body.ang_vel_in(world).to_matrix(world)).T)\n", "SYM.pprint(PlugAllIn(ω['B'][v] - body.ang_vel_in(body).to_matrix(body)).T)\n", "# answer: both!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can also check for the same ambiguity as above:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 0 0]\n", "[0 0 1.0]\n" ] } ], "source": [ "# In contrast, the time derivatives depend on reference frame\n", "SYM.pprint(PlugAllIn(ω['B'][v]).T)\n", "SYM.pprint(PlugAllIn(trafo['B']['S']*ω['S'][v]).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Second observation ****: \n", "- The angular kinematics are analogously ambiguous, \n", "- however, here it is not so clear that the rotation of the object around its COM is really zero. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As with linear measures, I keep $\\tilde{\\omega}_{B}$, $\\dot{\\tilde{\\omega}}_{B}$ and $\\ddot{\\tilde{\\omega}}_{B}$ for distinction:\n", "\n", "$\\tilde{\\omega}_{B} = R_{BS}\\omega_{S}$\n", "\n", "\n", "The $\\omega_{B}$ (the zero one) might be completely artificial. \n", "Most textbooks I have seen just write $\\omega_{S}$ (and actually mean: $\\dot{\\tilde{\\omega}}_{B} = R_{BS}\\dot{\\omega}_{S}$) for the kinematic equations. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But which ones of the two *should* we take? And is the choice identical to that for the linear motion parameters?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Chosing the Correct Vectors\n", "\n", "Let's birefly check with the **Coriolis Force,** as defined by [Wikipedia](https://en.wikipedia.org/wiki/Centrifugal_force?target=_blank#Force) (but with explicit dots for derivatives):\n", "\n", "$F_{cor} = -2m \\dot{\\omega} \\times \\dot{x}$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use $\\dot{{x}}_{B}$ and $\\dot{{\\omega}}_{B}$:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 0, 0]])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F_cor = -2*m*(ω['B'][v]).cross(x['B']['O'][v])\n", "PlugAllIn(F_cor).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or we use $\\dot{\\tilde{x}}_{B}$ and $\\dot{\\tilde{\\omega}}_{B}$:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 2.0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 2.0, 0]])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F_cor = -2*m*((trafo['B']['S'] * ω['S'][v]).cross(trafo['B']['S'] * x['S']['O'][v]))\n", "PlugAllIn(F_cor).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are a third and fourth option: the \"mixed ones\".\n", "\n", "One with $\\dot{{x}}_{B}$ but $\\dot{\\tilde{\\omega}}_{B}$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 0, 0]])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F_cor = -2*m*(trafo['B']['S'] * ω['S'][v]).cross(x['B']['J'][v])\n", "PlugAllIn(F_cor).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And yet another with $\\dot{\\tilde{x}}_{B}$ and $\\dot{{\\omega}}_{B}$" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 0, 0]])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F_cor = -2*m*(ω['B'][v]).cross(trafo['B']['S'] * x['S']['O'][v])\n", "PlugAllIn(F_cor).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Looks like we can rule out the second option. \n", "However, some of the others are certainly zero by coincidence, because our example is simple. \n", "\n", "This shows that there is definitely more to work through. *Spoiler*: number three (the first mixed strategy, $\\dot{{x}}_{B}$ and $\\dot{\\tilde{\\omega}}_{B}$) is most likely the correct option. \n", "But I will only demonstrate that three notebooks ahead." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I focused here on Coriolis because it is the decisive one: Euler- and Centrifugal Force are indifferent to the observation  above because they only require a position vector for $P$. \n", "\n", "But there is a general way of how to derive fictitious forces, which others have explained better than I could." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# The Five-Term Acceleration Equation\n", "\n", "Obviously, the faulty Coriolis force gave me a trouble. I did some more web research and found [this notebook by Frank Owen](http://www.aoengr.com/Dynamics/RotatingReferenceFrame.pdf?target=_blank), who runs a company that does \"engineering consultancy\". \n", "\n", "That document derives the accelerations which appear in a rotating reference frame with sufficient mathematical rigor. \n", "He takes a position of an arbitrary point on a rotating, not necessarily rigid body. \n", "The author then deconstructs the position vector into\n", "- one vector in $\\{S\\}$ from the origin to the center of mass of the body\n", "- a second vector in $\\{B\\}$ from the center of mass to the point.\n", "\n", "Deriving that combined vector in all its elements with the product rule gives all the fictitious accelerations. \n", "Because that procedure is sufficiently rigorous, it is obvious that there can't be any other fictitious forces (answering a question that I've been asking myself early one, *why are these three supposed to be all?*). \n", "The outcome of the maths are vectorial formulas that can be easily applied in the sympy setting we established above: \n", "\n", "(1) The acceleration of the rigid body COM (abbreviated C, with distance $OC$ from the origin) where the reference frame is attached.\n", "\n", "$a_{A} = a_{com, S} = \\frac{d^2(OC)}{dt^2}$\n", "\n", "\n", "(2) The movement of a focus point in that reference frame.\n", "\n", "$a_{xy} = a_{com, B} = = \\frac{d^2(CP)}{dt^2}$\n", "\n", "(where $P$ is a point different from the COM; this term is zero in our simplified example)\n", "\n", "\n", "(3) Coriolis acceleration\n", "\n", "$a_{C} = -2m\\cdot (R_{BS} \\dot{\\omega}_{S}) \\times \\dot{x}_{S} = -2m\\dot{\\tilde{\\omega}}_{B} \\times \\dot{x}_{B}$\n", "\n", "\n", "(4) tangential acceleration (corresponding to Euler Force)\n", "\n", "$a_{t} = -m\\cdot (R_{BS} \\ddot{\\omega}_{S}) \\times x_{S} = -m\\ddot{\\tilde{\\omega}}_{B} \\times {x}_{B}$\n", "\n", "\n", "(5) normal acceleration (corresponding to Centrifugal Force)\n", "\n", "$a_{n} = -m\\cdot (R_{BS} \\dot{\\omega}_{S}) \\times \\left( (R_{BS} \\dot{\\omega}_{S}) \\times \\dot{x}_{S} \\right) = -m\\dot{\\tilde{\\omega}}_{B} \\times \\left( \\dot{\\tilde{\\omega}}_{B} \\times \\dot{{x}}_{B} \\right)$\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is more explanation worth reading on the particularly non-intuitive Coriolis force in a [separate notebook](http://www.aoengr.com/Dynamics/CoriolisAcceleration.pdf?target=_blank).\n", "\n", "**What is the difference between Owen's formulas and the Wikipedia one?**\n", "\n", "The difference is that Owen keeps track of all the reference frames and definitions. \n", "From his document, it is clear which \"$r$\", \"$x$\" and \"$\\omega$\" he uses. He also adds accurate sketches that visualize the outcome. \n", "\n", "I wrote about these two different possible $x_B$  and $\\omega_B$  above: one that is zero due to the definition of the reference frame, and one that is the transformed $\\tilde{\\omega}_{B} = R_{BS} \\omega_{S}$. \n", "Seeing them in the acceleration formulas now, you might be tempted to forget about $\\omega_{B}$, for example, and use $\\tilde{\\omega}_{B}$ all over the place (*yes*, I actually did that mistake, and it cost me hours to find out). \n", "However, we actually need both: \n", "recall that, for the Balance Equation, we have to use $x_B$ and $\\omega_{B}$, and not the transformed $R_{BS} x_{S}$/$R_{BS} \\omega_{S}$.\n", "\n", "**Take home message:** \n", "Make sure, in this and all that follows, in what reference frame parameters are defined and whether they were transformed. \n", "Oh, yes, and never trust Wikipedia." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Overview: Fictitious Forces\n", "\n", "We'll assume for now (and try to disprove) that the third option from above is correct. \n", "Also, because the choice should be consistent for all three fictitious forces, I will apply the same logic to the others as well.\n", "\n", "This gives:\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The **Euler force**\n", "\n", "$F_e = -m\\ddot{\\tilde{\\omega}}_{B} \\times {x}_{B}$\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 0, 0]])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F_e_B = -m*(trafo['B']['S'] * ω['S'][a]).cross(x['B']['J'][p])\n", "PlugAllIn(F_e_B).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "the **Centrifugal Force:**\n", "\n", "$F_{cf} = -m\\dot{\\tilde{\\omega}}_{B} \\times \\left(\\dot{\\tilde{\\omega}}_{B} \\times {x}_{B} \\right)$\n" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 1.0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 1.0, 0]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F_cf_B = -m*(trafo['B']['S'] * ω['S'][v]).cross((trafo['B']['S'] * ω['S'][v]).cross(x['B']['J'][p]))\n", "PlugAllIn(F_cf_B).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "the **Coriolis Force:** (never forget the factor two!)\n", "\n", "$F_{cor} = -2m\\dot{\\tilde{\\omega}}_{B} \\times \\dot{{x}}_{B}$\n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 0, 0]])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F_cor_B = -2*m*(trafo['B']['S'] * ω['S'][v]).cross(x['B']['J'][v])\n", "PlugAllIn(F_cor_B).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are the values we expect to see for our initial example. \n", "And for those of my readers who noticed: indeed, I intentionally haven't discussed the reference point for linear motion parameters here, because 'J' and 'O' coincide. \n", "I'll save that for a later notebook." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Balance Of Wrenches\n", "## In the Inertial Frame\n", "\n", "We repeat the calculation of Balance of Wrenches for the inertial frame, because (if correct) that result must be correct, independent of fictitious forces. \n", "It is always a good starting point. \n", "\n", "Because the balance equations are based on a first principle, **they must give the same result, no matter which reference frame was used for calculation**. We can and should compare later." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(1) Assemble all external wrenches that change the momentum of the rigid body. \n", "Here, this will only be the joint wrench we are interested in. " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[f_{JS1} f_{JS2} f_{JS3} m_{JS1} m_{JS2} m_{JS3}]\n" ] } ], "source": [ "force_components_S = SYM.symbols('f_{JS1:4}', real = True)\n", "moment_components_S = SYM.symbols('m_{JS1:4}', real = True)\n", "\n", "joint_wrench_S = WT.Wrench.FromComponents(world, joint, force_components_S, moment_components_S)\n", "\n", "SYM.pprint(joint_wrench_S.Matrix().T)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(2) Collect the dynamic wrench." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}\\operatorname{sin}\\left(t\\right) & - \\operatorname{cos}\\left(t\\right) & 0 & 0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[sin(t), -cos(t), 0, 0, 0, 0]])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dynamic_force_S = m*x['S']['O'][a]\n", "dynamic_moment_S = I_S['P']*ω['S'][a] + ω['S'][v].cross(I_S['P']*ω['S'][v])\n", "\n", "# assembling the dynamic wrench\n", "dynamic_wrench_S = WT.Wrench.FromMatrices(world, joint, dynamic_force_S, dynamic_moment_S \\\n", " )\n", "# display the outcome\n", "PlugAllIn(dynamic_wrench_S.Matrix()).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(3) Solve the balance equations for the unknown joint wrench." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[sin(t) -cos(t) 0 0 0 0]\n" ] } ], "source": [ "# get the equations (one per component)\n", "equations_S = dynamic_wrench_S.Equate(joint_wrench_S)\n", "\n", "# desired outcome variables\n", "components_S = [*force_components_S, *moment_components_S]\n", "\n", "# extract the solution\n", "solutions = {}\n", "for param, sol in SYM.solve(equations_S \\\n", " , components_S \\\n", " ).items():\n", " solutions[param] = sol\n", "\n", "# transform the result into a wrench\n", "W_JS = WT.Wrench.FromMatrix(world, joint, SYM.Matrix([[solutions[cmp]] for cmp in components_S]) )\n", "SYM.pprint(PlugAllIn(W_JS.Matrix()).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Same result as before. \n", "Keep it stored, e.g. by writing it on your hand or forehead." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Naïve Balance in the Body Frame\n", "The procedure is the same in the body frame. \n", "As mentioned somewhere, always double check that you get the reference points and reference frames right!\n", "\n", "As you will see, it's not enough to change all 'S' to 'B', the outcome will not match." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[f_{JB1} f_{JB2} f_{JB3} m_{JB1} m_{JB2} m_{JB3}]\n" ] } ], "source": [ "force_components_B = SYM.symbols('f_{JB1:4}', real = True)\n", "moment_components_B = SYM.symbols('m_{JB1:4}', real = True)\n", "\n", "joint_wrench_B = WT.Wrench.FromComponents(body, joint, force_components_B, moment_components_B)\n", "\n", "SYM.pprint(joint_wrench_B.Matrix().T)\n" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 0, 0, 0, 0, 0]])" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dynamic_force_B = m*x['B']['O'][a]\n", "dynamic_moment_B = I_B['P']*ω['B'][a] + ω['B'][v].cross(I_B['P']*ω['B'][v])\n", "\n", "# assembling the dynamic wrench\n", "dynamic_wrench_B = WT.Wrench.FromMatrices(body, pointmass, dynamic_force_B, dynamic_moment_B \\\n", " )\n", "# display the outcome\n", "PlugAllIn(dynamic_wrench_B.Matrix()).T" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 0 0 0 0 0]\n" ] } ], "source": [ "# get the equations (one per component)\n", "equations_B = dynamic_wrench_B.Equate(joint_wrench_B)\n", "\n", "# desired outcome variables\n", "components_B = [*force_components_B, *moment_components_B]\n", "\n", "# extract the solution\n", "solutions = {}\n", "for param, sol in SYM.solve(equations_B \\\n", " , components_B \\\n", " ).items():\n", " solutions[param] = sol\n", "\n", "# transform the result into a wrench\n", "W_JB = WT.Wrench.FromMatrix(body, joint, SYM.Matrix([[solutions[cmp]] for cmp in components_B]) )\n", "SYM.pprint(PlugAllIn(W_JB.Matrix()).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our Centrifugal Wrench is gone." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding Fictitious Wrenches\n", "\n", "Recall the definition of fictitious forces: they have to be introduced when going to the non-static reference frame in order for the balance equations to be correct. \n", "\n", "Here you go. This is exactly what we have to do. \n", "Of course we use wrenches, and as elaborated above we use certain variants of the kinematic elements (namely $\\tilde{\\omega}_{B}$ and $x_{B}$). \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here comes **Euler**, the \"tangential\" one (which actually means \"not normal\", as he certainly was):" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 0 0 0 0 0]\n" ] } ], "source": [ "euler_wrench = WT.Wrench(body, pointmass \\\n", " , MatrixToVector(-m*((trafo['B']['S'] * ω['S'][a]).cross(x['B']['O'][p])), body) \\\n", " )\n", "\n", "SYM.pprint(PlugAllIn(euler_wrench.Matrix()).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here comes the normal component, also called **centrifugal**, so hold your seat:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 1.0 0 0 0 0]\n" ] } ], "source": [ "centrifugal_wrench = WT.Wrench(body, pointmass \\\n", " , MatrixToVector(-m*( \\\n", " (trafo['B']['S'] * ω['S'][v]).cross( \\\n", " (trafo['B']['S'] * ω['S'][v]).cross(x['B']['O'][p]) \\\n", " ) \\\n", " ), body) \\\n", " )\n", "\n", "SYM.pprint(PlugAllIn(centrifugal_wrench.Matrix()).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here the one that everyone remembers from, but noone understands at school, **Coriolis**:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 0 0 0 0 0]\n" ] } ], "source": [ "coriolis_wrench = WT.Wrench(body, pointmass \\\n", " , MatrixToVector(-2*m*((trafo['B']['S'] * ω['S'][v]).cross(x['B']['J'][v])), body) \\\n", " )\n", "\n", "SYM.pprint(PlugAllIn(coriolis_wrench.Matrix()).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note a little computational trick that I used above in initializing the wrenches:\n", "- I define them at the point mass (i.e. the COM).\n", "- Because they physically apply there, the moment about that specific point is zero.\n", "- the WT.Wrench constructor sets the moment to zero by default, so I don't need to specify.\n", "- By translating the wrench to other points, a moment can appear.\n", "\n", "But pay attention: defining wrenches by only the forces requires you to set the correct point of application and the correct reference frame. \n", "That is easier said than done.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In real problems, one will be less interested in the ones above. We are actually after the **joint wrench**:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[f_{JB1} f_{JB2} f_{JB3} m_{JB1} m_{JB2} m_{JB3}]\n" ] } ], "source": [ "force_components_B = SYM.symbols('f_{JB1:4}', real = True)\n", "moment_components_B = SYM.symbols('m_{JB1:4}', real = True)\n", "\n", "joint_wrench_B = WT.Wrench.FromComponents(body, joint, force_components_B, moment_components_B)\n", "\n", "SYM.pprint(joint_wrench_B.Matrix().T)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then there is the left handside with the **dynamic wrench**:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 0, 0, 0, 0, 0]])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dynamic_force_B = m*x['B']['O'][a]\n", "\n", "dynamic_moment_B = I_B['O']*ω['B'][a] + ω['B'][v].cross(I_B['O']*ω['B'][v])\n", "\n", "# assembling the dynamic wrench\n", "dynamic_wrench_B = WT.Wrench.FromMatrices(body, pointmass, dynamic_force_B, dynamic_moment_B \\\n", " ).Translate(joint)\n", "\n", "\n", "PlugAllIn(dynamic_wrench_B.Matrix()).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And, finally, the **balance equations**:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0 -1.0 0 0 0 0]\n" ] } ], "source": [ "# balance equations\n", "equations_B = dynamic_wrench_B.Equate(joint_wrench_B + euler_wrench + centrifugal_wrench + coriolis_wrench)\n", "\n", "\n", "# the usual solving procedure\n", "components_B = [*force_components_B, *moment_components_B]\n", "\n", "solutions = {}\n", "for param, sol in SYM.solve(equations_B \\\n", " , components_B \\\n", " ).items():\n", " solutions[param] = SYM.simplify(sol)\n", "\n", "\n", "# and here comes the body frame joint wrench\n", "W_JB = WT.Wrench.FromMatrix(body, joint, SYM.Matrix([[solutions[cmp]] for cmp in components_B]) )\n", "SYM.pprint(PlugAllIn(W_JB.Matrix()).T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Still different from the result above ... ... because we need to transform!" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}\\operatorname{sin}\\left(t\\right) & - \\operatorname{cos}\\left(t\\right) & 0 & 0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[sin(t), -cos(t), 0, 0, 0, 0]])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "PlugAllIn(W_JB.express(world).Matrix()).T" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}\\operatorname{sin}\\left(t\\right) & - \\operatorname{cos}\\left(t\\right) & 0 & 0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[sin(t), -cos(t), 0, 0, 0, 0]])" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "PlugAllIn(W_JS.Matrix()).T" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0 & 0 & 0\\end{matrix}\\right]$" ], "text/plain": [ "Matrix([[0, 0, 0, 0, 0, 0]])" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# because the \"PlugAllIn\" could cheat us, here is the ultimate check:\n", "PlugAllIn((W_JS.express(body) - W_JB ).Matrix()).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Check!**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Summary\n", "\n", "I have encountered several types of definitions of fictitious forces during my initial attempts.\n", "\n", "The first type is a practical one, which goes along the lines of \"*fictitious forces arise because an observer who is captured in the rotating frame is blind to the movement of the frame*\". \n", "Yes, that is a valid description, but a disappointing definition: What else might the observer be blind to? How can we calculate something we are blind for? And why can't I not always switch perspective in cases where I'm in charge of the data?\n", "\n", "The second type of definition is a technical one: \"*fictitious forces need to be introduced to make balance equations work in the rotating reference frame*\". \n", "That is also a valid point, yet it leaves similar questions as above: Which ones do I have to include when? Where do the formulas come from? \n", "\n", "I have hinted at these aspects above. \n", "They certainly help to develop an intuition about fictitious forces, but for me they were not satisfactory. \n", "You can find these two categories scattered all over the [Wikipedia articles about fictitious forces](https://en.wikipedia.org/wiki/Fictitious_force?target=_blank). \n", "Despite my rants herein, I'm not disappointed about Wikipedia, it is just sad that this confirmed my prior on that \"Encyclopedia\". \n", "But I also see myself as a source of the problem: whereas for trained physicists it is often obvious which and how vectors add up, my humble, untrained mind needs some more explicit guidance.\n", "\n", "\n", "Hence, I found a third definition most convincing, yet it is more like a prescription than a definition. It is the one which I found in the **workbooks by Frank Owen** (five-term acceleration, see references). \n", "- Conceive a position vector to a point on a rigid body as a superimposition of the position vector of the COM (in static frame), and the point relative to COM (in body frame). \n", "- Calculate component-wise time derivatives of that superimposed vector; use product rule, and make sure to incorporate time changes of reference frame basis vectors. \n", "- In the outcome, two acceleration components are intuitive cartesian superimpositions; the others depend on cross products of angular and linear motion parameters. The latter ones are irritatingly labelled \"fictitious\". \n", "- A corollary rule of thumb for those three: angular motion parameters ($\\omega$ and derivatives) stem from the static frame, whereas linear motion parameters ($x$ and derivatives) are best seen in the body frame. \n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "This notebook has solved some puzzles which I encountered when starting to learn about fictitious forces. \n", "Make sure to keep track of your reference frames and reference points. \n", "Good (i.e. explicit, structured, well commented) programming style can facilitate this. \n", "And remember that, even if you think it is obvious which parameters you took in which transformation, make sure that others can reproduce your line of thought.\n", "\n", "With these tools at our hands, we can advance to elementary examples. I'll start by **isolating Euler**, i.e. with an example that should show the Euler force in isolation, in the [next notebook](../id_lmx1_stutter)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# References\n", "\n", "- **Dumas**, R., Aissaoui, R., & de Guise, J. A. (**2004**). *A 3D generic inverse dynamic method using wrench notation and quaternion algebra.* Computer methods in biomechanics and biomedical engineering, 7(3), 159-166. [https://doi.org/10.1080/10255840410001727805](https://doi.org/10.1080/10255840410001727805?target=_blank)\n", "\n", "- **Dumas**, R. (**2019**). *3D Kinematics and Inverse Dynamics* , [MATLAB Central File Exchange](https://www.mathworks.com/matlabcentral/fileexchange/58021-3d-kinematics-and-inverse-dynamics?target=_blank). Accessed July 1, 2020. \n", "\n", "- **Lynch**, K. M., & Park, F. C. (**2017**). *Modern Robotics.* Cambridge University Press. ISBN 9781107156302. [http://www.modernrobotics.org](http://www.modernrobotics.org?target=_blank) \n", "\n", "- **Meurer** A., Smith C.P., Paprocki M., Čertík O., Kirpichev S.B., Rocklin M., Kumar A., Ivanov S., Moore J.K., Singh S., Rathnayake T., Vig S., Granger B.E., Muller R.P., Bonazzi F., Gupta H., Vats S., Johansson F., Pedregosa F., Curry M.J., Terrel A.R., Roučka Š., Saboo A., Fernando I., Kulal S., Cimrman R., Scopatz A. (**2017**). *SymPy: symbolic computing in Python.* PeerJ Computer Science 3:e103 [https://doi.org/10.7717/peerj-cs.103](https://doi.org/10.7717/peerj-cs.103?target=_blank)\n", "\n", "- **Moore**, J., McMurry, R., Milam, B. (**2016**). *Simulating Robot, Vehicle, Spacecraft, and Animal Motion.* SciPy 2016 Tutorial. [https://www.youtube.com/watch?v=r4piIKV4sDw](https://www.youtube.com/watch?v=r4piIKV4sDw?target=_blank). Accessed July 1, 2020.\n", "\n", "- **Owen**, F. (undated). *Rotating reference frame and the five-term acceleration equation*. Alpha Omega Engineering, Inc. [web document](http://www.aoengr.com/Dynamics/RotatingReferenceFrame.pdf?target=_blank). Accessed July 1, 2020." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.3" } }, "nbformat": 4, "nbformat_minor": 4 }