User derivatives
----------------

.. container:: python

   .. highlight:: python

-  It is possible to add additional constitutive equations to the
   equation of motion

   -  This can be used to deal with multiphysics problem involving
      electrical circuits, pneumatics, hydraulics, …
   -  The additional state equations will be added to the set of
      mechanical equations

-  In practice

   -  The user must declare the state in a user model
   -  The user must compute the state derivative equation in the
      user_Derivatives function
   -  The state value can be accessed in any user function of the
      project

Back to the spring-pendulum example
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. figure:: figure/userDeriv_example_torque.png
   :alt: A torque coming from the electrical motor is applied on the crank

   Pendulum spring illustration with torque

The exercise consists in adding an electrical DC motor to the pendulum
example without the wall:

-  The motor will be connected to the revolute joint of the crank
-  The motor obeys the following relation

   -  Electrical circuit equation:
      \ :math:`U_{mot} = R_{mot}*i_{mot}+k_{\phi}*\omega_{mot} + L \frac{di_{mot}}{dt}`\ 
   -  Torque equation: \ :math:`T_{mot} = k_{\phi} * i_{mot}`\ 
   -  A reductor is inserted between the motor and the axle:
      \ :math:`\begin{matrix}  T_{axle} = \rho * T_{mot}\\  \omega_{axle} = \omega_{mot}/\rho  \end{matrix}`\ 

The parameter values are:

-  \ :math:`U = 48\ V`\ 
-  \ :math:`k_{\phi} = 48.6\ mNm/A`\ 
-  \ :math:`\rho = 200`\ 
-  \ :math:`R_{mot} = 4.49\ \Omega`\ 
-  \ :math:`L = 0.573\ mH`\ 

Step 1: Draw your multibody system
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

-  Open the PendulumSpring model in MBsysPad
-  Add a user model for the motor

   -  Add a scalar variable for each parameter
   -  Add a state variable for the current

-  Set the nature of the :

   -  crank joint to independent;
   -  pendulum joint to dependent;
   -  If you had modified the joint nature in the code in the previous
      section you have to update the code, modifications done in
      MBsysPad are overwritten.

.. figure:: figure/userDeriv_snapshot_userModel.png
   :alt: User state variable must be introduces via the user model dialog box

   snapshot of the user model gui for introducing state variable

..

   REMARK

   Several state variables can be introduced either by introducing
   several values for a single state variable (state vector) or by
   defining several parameters with “state” type.

..

   WARNING:

   You can have two state variables with the same name in two different
   user model. However it is recommended to set an unique name for each
   state variable.

Step 3: Write your user function
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. container:: python

   -  Edit the *user_Derivative* function and introduce the state
      equation for the current:

   .. code:: python

      def user_derivatives(mbs_data):
         #... 
         
         id_crank = mbs_data.joint_id['R2_crank']
         rho = mbs_data.user_model['Motor']['rho']
         U   = mbs_data.user_model['Motor']['U']
         Rmot= mbs_data.user_model['Motor']['Rmot']
         kphi= mbs_data.user_model['Motor']['Kphi']
         L   = mbs_data.user_model['Motor']['L']
         
         # The user derivative index are retrieved in the user model state value.
         # This is important if your model contain more than one user state.
         id_state =  mbs_data.user_model['Motor']['i'][1]
         
         omMot = -mbs_data.qd[id_crank]*rho
         iMot = mbs_data.ux[id_state]
         mbs_data.uxd[id_state] = (U-Rmot*iMot-kphi*omMot)/L
         
         #...

   ..

      REMARK

      The user state vector is accessed via the ux field of the MbsData
      class

      You can access to the ID of a state variable via the user model
      structure. If the state variable is a vector, you retrieve the ID
      of the first value of the vector in [1], second in [2], etc.

   -  Edit the *user_JointForces* function and introduce the torque
      equation:

   .. code:: python

      def user_JointForces(mbs_data,tsim):
         #... 
         
         id_crank = mbs_data.joint_id['R2_crank']
         id_state =  mbs_data.user_model['Motor']['i'][1]

         iMot = mbs_data.ux[id_state]
         rho = mbs_data.user_model['Motor']['rho']
         kphi= mbs_data.user_model['Motor']['Kphi']

         mbs_data.Qq[id_crank] = -rho*kphi*iMot
         
         #...
         return

   ..

      REMARK

      The dynamics of the motor required a small time-step (smaller than
      the default one) or an adaptive integrator.

      As result the simulation may be long. To reduce the simulation
      time you can multiply the motor inductance by a factor 10 without
      noticeable changes on the results

   The time history of the motor current can be directly plotted from
   the field ``ux`` of ``MbsDirdyn.results``.

..

   REMARK:

   To let the motor rotate freely, the wall must be removed (by
   modifying the external forces).

.. container:: python

   ..

      REMARK:

      The resolution of the motor dynamics requires a smaller time step
      than before. The computation time will increase.

      You must edit the *main* function and decrease it :
      ``mbs_dirdyn.set_options(dt0=1e-4)`` (except if you changed the
      inductance value).

Check the results
^^^^^^^^^^^^^^^^^

Plot the graph of the joint position (results ares available in
resultsR/ folder) and check your results with the following graph.
Please note that the wall force has been removed.

.. figure:: figure/userDeriv_curves_q1q2.png
   :alt: Plot the time history of the joint position of the pendulum spring example

   Time history of the joint position