# User derivatives¶

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¶

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: \(U_{mot} = R_{mot}*i_{mot}+k_{\phi}*\omega_{mot} + L \frac{di_{mot}}{dt}\)

Torque equation: \(T_{mot} = k_{\phi} * i_{mot}\)

A reductor is inserted between the motor and the axle: \(\begin{matrix} T_{axle} = \rho * T_{mot}\\ \omega_{axle} = \omega_{mot}/\rho \end{matrix}\)

The parameter values are:

\(U = 48\ V\)

\(k_{\phi} = 48.6\ mNm/A\)

\(\rho = 200\)

\(R_{mot} = 4.49\ \Omega\)

\(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.

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¶

Edit the *user_Derivative* function and:

Include user_model.h and user_all_id.h;

Introduce the state equation for the current:

```
//...
#include "user_model.h"
#include "user_all_id.h"
//...
void user_Derivative(MbsData *mbs_data)
{
UserModel *um = mbs_data->user_model;
double rho = um->Motor.rho;
double U = um->Motor.U;
double Rmot= um->Motor.Rmot;
double kphi= um->Motor.Kphi;
double L = um->Motor.L;
int id_crank = R2_crank_id;
int id_state = um->Motor.i[1];
double omAxe = mbs_data->qd[id_crank];
double omMot = -omAxe * rho;
double 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 mbs Data class

You should be able to access to the ID of a state variable via the user model. The um->Motor.i[1] should return you the id of the state variable in the mbs_data->ux (and uxd) vector. You can then access to the state value with this command:

`mbs_data->ux[ um->Motor.i[1] ]`

.If the state variable is a vector of several values, you retrieve the ID of the nth value with um->Motor.i[n]. In [0], it should be the length of the vector.

Edit the

*user_JointForces*function and:Include user_model.h and user_all_id.h;

Introduce the torque equation:

```
//...
#include "user_model.h"
#include "user_all_id.h"
//...
double* user_JointForces(MbsData *mbs_data, double tsim)
{
// Force due to the spring in the pendulum-slider joint
//...
// Force due to the electrical motor
UserModel *um = mbs_data->user_model;
double rho = um->Motor.rho;
double kphi = um->Motor.Kphi;
id = R2_crank_id;
int id_state = 1;
double imot = mbs_data->ux[id_state];
mbs_data->Qq[id] = - rho * kphi * imot;
//...
}
```

The time history of the motor current (and its derivate) is stored in
your resultR/ folder in the *dirdyn_ux.res* and *dirdyn_uxd.res*
files .

REMARK:

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

WARNING:

The resolution of the motor dynamics requires a smaller time step than before. You must edit the

mainfunction and decrease it :`mbs_dirdyn->options->dt0 = 1e-4;`

.This must be done because the integrator use constant time step.

Another solution is to use another integrator as described in the tips and tricks.

### 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.