Dynamic Exports

If the generalized coordinates that you define are independent, the dynamic equations governing the system response constitute a set of ordinary differential equations (ODEs).

If the coordinates are dependent, the ODEs for the dynamic response can only be solved simultaneously with the nonlinear algebraic equations from the position-level kinematic constraints; the resulting system equations constitute a set of differential-algebraic equations (DAEs). With MapleSim/Multibody, you can represent the constraint reactions by their actual forces and moments, or indirectly using Lagrange multipliers.

 Baumgarte Constraint Stabilization One method of converting these DAEs to ODEs is to replace the position constraints with the acceleration constraints, which are then numerically integrated simultaneously with the ODEs from the dynamic equations. However, during the integration process, the accumulation of numerical errors leads to violations in the position constraint equations. (Visually, the cotree joints will float apart.) Baumgarte proposed a method to stabilize these constraints, by combining the position, velocity, and acceleration constraints into a single expression:   $\mathrm{\chi }+2\mathrm{\alpha }\mathrm{\Psi }+{\mathrm{\beta }}^{2}\mathrm{\Phi }=0$   which can be written as a linear equation in terms of the accelerations:        ${\mathrm{\Psi }}_{p}\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}=e-2\mathrm{\alpha }\mathrm{\Psi }-{\mathrm{\beta }}^{2}\mathrm{\Phi }$ $=\mathrm{\epsilon }$   By integrating this expression for the accelerations, the Baumgarte parameters, α and β, act to stabilize the constraints at the position  level. These parameters are set using the command:   > Model:-SetBaumgarte([alpha,beta]);   where Model is the module returned from the GetMultibody command, and α and β must be given numerical values prior to numerical solution. Typical values range from 1 to 10, and depend upon the characteristics of the particular system under study.   As explained below, Baumgarte stabilization is automatically used in the exports for constrained dynamic systems. If the Baumgarte parameters are not set, default values of 0 are assumed (that is, no constraint stabilization).

Unconstrained Dynamics (ODEs in Independent Coordinates)

If the generalized coordinates, q, and generalized speeds, p, are independent, the dynamic equations take the "unconstrained" form:

$M\frac{\mathrm{dp}}{\mathrm{dt}}=F$

where $M$ is the mass matrix and $F$ contains external loads and quadratic velocity terms. The full set of system equations are formed by combining these dynamic equations with the kinematic transformations, giving 2n ordinary differential equations in 2n unknowns:

$[\begin{array}{cc}M& 0\\ 0& 1\end{array}][\begin{array}{c}\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ \frac{ⅆq}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\end{array}]=[\begin{array}{c}F\\ h\end{array}]$

These unconstrained dynamic equations, and their variables, are summarized in the table below.

 Expression BuildEQs Export Description $M\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}-F=0$ GetDynEQs() Dynamic equations in n independent coordinates and speeds $M$ xM Mass matrix (n x n) $F$ vF Right-hand side of dynamic equations (n x 1 column matrix) $[\begin{array}{cc}M& 0\\ 0& 1\end{array}][\begin{array}{c}\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ \frac{ⅆq}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\end{array}]=[\begin{array}{c}F\\ h\end{array}]$ GetSysODEs() Complete set of 2n ODEs for systems modeled by n independent coordinates and speeds $[\begin{array}{cc}M& 0\\ 0& 1\end{array}],[\begin{array}{c}\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ \frac{ⅆq}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\end{array}],[\begin{array}{c}F\\ h\end{array}]$ GetAYBSysODEs() Returns the three matrices comprising the system ODEs.  Note: The column matrix Y corresponds to the derivative of state variables, vXdot.

Constrained Dynamics (DAEs with Constraint Reactions)

If the generalized coordinates, q, and generalized speeds, p, are dependent, the "augmented" dynamic equations can be written with the constraint reactions appearing explicitly (AugType option set to Reaction):

$M\frac{\mathrm{dp}}{\mathrm{dt}}+{C}^{T}f=F$

where f contains the reaction loads in the cotree joints that enforce the kinematic constraints, and C is the corresponding coefficient matrix.  These constrained dynamic equations, and their variables, are summarized in the table below.

 Expression BuildEQs Export Description $M\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}+{C}^{T}f-F=0$ GetDynEQs() Dynamic equations in n dependent coordinates and speeds $M$ xM Mass matrix (n x n) $F$ vF Right-hand side of dynamic equations (n x 1 column matrix) $f$ vReactions Reaction forces and moments in the cotree joints that enforce the kinematic constraint equations (m x 1 column matrix) $C$ xC Coefficient matrix of constraint reactions (m x n) $M\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}+{C}^{T}f-F=0$ $\mathrm{\Phi }\left(q,t\right)=0$ $\frac{ⅆq}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}=h\left(p,q,t\right)$ GetSysEQs() Complete set of 2n+m DAEs in terms of 2n coordinates and speeds and m constraint reactions $[\begin{array}{ccc}M& 0& {C}^{T}\\ 0& 1& 0\\ {\mathrm{\Psi }}_{p}& 0& 0\end{array}][\begin{array}{c}\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ \frac{ⅆq}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ f\end{array}]=[\begin{array}{c}F\\ h\\ \mathrm{\epsilon }\end{array}]$ GetSysODEs() Conversion of DAEs to 2n+m ODES, by replacing position constraints with Baumgarte form of acceleration constraints $[\begin{array}{ccc}M& 0& {C}^{T}\\ 0& 1& 0\\ {\mathrm{\Psi }}_{p}& 0& 0\end{array}],[\begin{array}{c}\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ \frac{ⅆq}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ f\end{array}],[\begin{array}{c}F\\ h\\ \mathrm{\epsilon }\end{array}]$ GetAYBSysODEs() Returns the three matrices comprising the system ODEs

Constrained Dynamics (DAEs with Lagrange Multipliers)

Alternatively, you can write the augmented dynamic equations with the constraint reactions represented implicitly by Lagrange multipliers (AugType option set to Lagrange):

$M\frac{\mathrm{dp}}{\mathrm{dt}}+{{\mathrm{\Psi }}_{p}}^{T}\mathrm{\lambda }=F$

where λ are the Lagrange multipliers (one for each constraint equation) and ${\mathrm{\Psi }}_{p}$ is the Jacobian matrix of the velocity constraint equations with respect to the generalized speeds p.  These constrained dynamic equations, and their variables, are summarized in the table below.

 Expression BuildEQs Export Description $M\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}+{{\mathrm{\Psi }}_{p}}^{T}\mathrm{\lambda }-F=0$ GetDynEQs() Dynamic equations in n dependent coordinates and speeds $M$ xM Mass matrix (n x n) $F$ vF Right-hand side of dynamic equations (n x 1 column matrix) $\mathrm{\lambda }$ vLambda Lagrange multipliers that indirectly represent the reaction forces and moments in the cotree joints that enforce the kinematic constraint equations (m x 1 column matrix) $M\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}+{{\mathrm{\Psi }}_{p}}^{T}\mathrm{\lambda }-F=0$ $\mathrm{\Phi }\left(q,t\right)=0$ $\frac{ⅆq}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}=h\left(p,q,t\right)$ GetSysEQs() Complete set of 2n+m DAEs in terms of 2n coordinates and speeds, and m Lagrange multipliers $[\begin{array}{ccc}M& 0& {{\mathrm{\Psi }}_{p}}^{T}\\ 0& 1& 0\\ {\mathrm{\Psi }}_{p}& 0& 0\end{array}][\begin{array}{c}\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ \frac{ⅆq}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ \mathrm{\lambda }\end{array}]=[\begin{array}{c}F\\ h\\ \mathrm{\epsilon }\end{array}]$ GetSysODEs() Conversion of DAEs to 2n+m ODES, by replacing position constraints with Baumgarte form of acceleration constraints $[\begin{array}{ccc}M& 0& {{\mathrm{\Psi }}_{p}}^{T}\\ 0& 1& 0\\ {\mathrm{\Psi }}_{p}& 0& 0\end{array}],[\begin{array}{c}\frac{ⅆp}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ \frac{ⅆq}{ⅆ\phantom{\rule[-0.0ex]{0.2em}{0.0ex}}t}\\ \mathrm{\lambda }\end{array}],[\begin{array}{c}F\\ h\\ \mathrm{\epsilon }\end{array}]$ GetAYBSysODEs() Returns the three matrices comprising the system ODEs

Examples

The example below shows an embedded MapleSim model of a planar slider-crank mechanism.  The form of the commands below assume that you are using a Worksheet attached to the MapleSim model through Add Apps or Templates > Templates > Worksheet. In this example we are going to name our multibody exports module SliderCrank.

 > $A≔\mathrm{MapleSim}:-\mathrm{LinkModel}\left(\right):$
 >
 ${"Analyzing system..."}$
 ${"Performing constraint analysis..."}$
 ${"The system has 1 degree\left(s\right) of freedom. It is modeled using 3 generalized coordinate\left(s\right) coupled by 2 algebraic constraint\left(s\right)."}$
 ${"Performing a dynamic analysis using an augmented reaction formulation - system variables shown below:"}$
 ${"vQ"}{,}\left[\begin{array}{c}{\mathrm{P1_s}}{}\left({t}\right)\\ {\mathrm{R1_theta}}{}\left({t}\right)\\ {\mathrm{R3_theta}}{}\left({t}\right)\end{array}\right]{,}{"vP"}{,}\left[\begin{array}{c}\frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{P1_s}}{}\left({t}\right)\\ \frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{R1_theta}}{}\left({t}\right)\\ \frac{{ⅆ}}{{ⅆ}{t}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{R3_theta}}{}\left({t}\right)\end{array}\right]{,}{"vPdot"}{,}\left[\begin{array}{c}\frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{P1_s}}{}\left({t}\right)\\ \frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{R1_theta}}{}\left({t}\right)\\ \frac{{{ⅆ}}^{{2}}}{{ⅆ}{{t}}^{{2}}}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\mathrm{R3_theta}}{}\left({t}\right)\end{array}\right]{,}{"vReactions"}{,}\left[\begin{array}{c}{\mathrm{Fx}}{}\left({t}\right)\\ {\mathrm{Fy}}{}\left({t}\right)\end{array}\right]$
 ${"Dynamic analysis complete."}$ (1.6.1)
 > ${\mathbf{SliderCrank}}{\mathbf{:-}}{\mathbf{xM}}{\mathbf{;}}$
 $\left[\begin{array}{ccc}\frac{5}{2}& 0& \frac{1}{4}{}\mathrm{sin}{}\left(\mathrm{R3_theta}{}\left(t\right)\right)\\ 0& \frac{3}{400}& 0\\ \frac{1}{4}{}\mathrm{sin}{}\left(\mathrm{R3_theta}{}\left(t\right)\right)& 0& \frac{129}{1000}\end{array}\right]$ (1.6.2)
 > ${\mathbf{SliderCrank}}{\mathbf{:-}}{\mathbf{vF}}{\mathbf{;}}$
 $\left[\begin{array}{c}-\frac{1}{4}{}{\left(\frac{ⅆ}{ⅆt}{}\mathrm{R3_theta}{}\left(t\right)\right)}^{2}{}\mathrm{cos}{}\left(\mathrm{R3_theta}{}\left(t\right)\right)\\ -\frac{2943}{10000}{}\mathrm{cos}{}\left(\mathrm{R1_theta}{}\left(t\right)\right)\\ \frac{981}{400}{}\mathrm{cos}{}\left(\mathrm{R3_theta}{}\left(t\right)\right)\end{array}\right]$ (1.6.3)