Advanced user constraints¶

The present tutorial is about an advanced use of the user constraints. The difference with the basic “Modelling Features” tutorial is the non-null quadratic term, i.e. \(\left(\dot{J}\dot{q}(q, \dot{q})\right)_{usr}\neq0\).

As a reminder, user constraints are defined as follows:

The user constraint is expressed in the form: \(h_{usr}(q) = 0\)
The user has to provide, for a given configuration:
- The value of the constraint: \(h_{usr}(q)\)
- The Jacobian matrix of the constrain: \(J_{usr}(q)\)
- The quadratic term: \(\left(\dot{J}\dot{q}(q, \dot{q})\right)_{usr}\)
The user constraints are implemented using the user_cons_hJ and user_cons_jdqd functions.

Rod-crank-piston example¶

The system to analyse is a simple rod-crank-piston system as it can be found in engines.

Rod-crank-piston system: definition¶

The model is composed of the following elements:

Three bodies
- the rod \(OA\)
  - Length: a = 150 mm
  - Center of mass \(G^1\) located in the middle of the segment \(OA\)
  - Mass: 0.3 kg
  - Inertia along axis \(\hat{I}_{3}\): 0.005 kg.m2
- the crank \(AB\)
  - Length: b = 300 mm
  - Center of mass \(G^2\) located in the middle of the segment \(AB\)
  - Mass: 0.7 kg
  - Inertia along axis \(\hat{I}_{3}\): 0.003 kg.m2
- the piston \(B\)
  - Center of mass \(G^3\) located in \(B\)
  - Mass: 0.35 kg
Three joints
- A revolute joint between the base and the rod
- A revolute joint between the rod and the crank
- A prismatic joint between the base and the piston
A ball cut between the crank and the piston
Forces
- External torque applied on the rod in \(O\): \(T\hat{I}_{3}\) with \(T=0.5\,Nm\)
- Friction force resistive to the piston motion: \(-c\dot{x}\hat{I}_{1}\) with \(c=5\,Ns/m\)
- No gravity in the system

The goal is to analyze the evolution over 5sec of the system with the piston starting at the top-dead-center configuration.

Although this system has three generalized coordinates, there is only one DOF, i.e. the piston translation. Two constraints are therefore necessary to link the rotations of both rod and crank to the piston translation.

The constraints can thus be written as:

\(h_1 \triangleq a~\cos \alpha + b~\cos \left( \alpha + \beta \right) - x = 0\)
\(h_2 \triangleq a~\sin \alpha + b~\sin \left( \alpha + \beta \right) = 0\)
\(\alpha\): rod rotation;
\(\beta\): crank rotation;
\(x\): piston rotation.

The Jacobian is given by

\(J = \left( \begin{array}{rrr} - a \sin\alpha - b \sin \left( \alpha + \beta \right) & -b \sin \left( \alpha + \beta \right) & -1 \\ a \cos\alpha + b \cos \left( \alpha + \beta \right) & b \cos \left( \alpha + \beta \right) & 0 \end{array} \right)\)

Its time derivative is equal to

\(\dot{J} = \left( \begin{array}{rrr} - a \dot \alpha \cos \alpha - b (\dot \alpha + \dot \beta) \cos(\alpha + \beta) & -b (\dot \alpha + \dot \beta) \cos(\alpha + \beta) & 0 \\ - a \dot \alpha \sin \alpha - b (\dot \alpha + \dot \beta) \sin(\alpha + \beta) & -b (\dot \alpha + \dot \beta) \sin(\alpha + \beta) & 0 \\\end{array} \right)\)

And the quadratic term \(\dot{J}\dot{q}\) is therefore

\(\begin{array}{ccl} \dot{J}\dot{q} & = & \left( \begin{array}{rrr} - a \dot \alpha \cos \alpha - b (\dot \alpha + \dot \beta) \cos(\alpha + \beta) & -b (\dot \alpha + \dot \beta) \cos(\alpha + \beta) & 0 \\ - a \dot \alpha \sin \alpha - b (\dot \alpha + \dot \beta) \sin(\alpha + \beta) & -b (\dot \alpha + \dot \beta) \sin(\alpha + \beta) & 0 \\\end{array} \right) \left( \begin{array}{c} \dot{\alpha}\\ \dot{\beta}\\ \dot{x} \end{array} \right) \\ \\ &=& \left( \begin{array}{c} - a \dot{\alpha}^2 \cos \alpha - b (\dot \alpha + \dot \beta) \dot{\alpha} \cos(\alpha + \beta) - b (\dot \alpha + \dot \beta) \dot{\beta} \cos(\alpha + \beta)\\ - a \dot{\alpha}^2 \sin \alpha - b (\dot \alpha + \dot \beta) \dot{\alpha} \sin(\alpha + \beta) - b (\dot \alpha + \dot \beta) \dot{\beta} \sin(\alpha + \beta)\\ \end{array} \right) \\ \\ &=& \left( \begin{array}{c} - a \dot{\alpha}^2 \cos \alpha - b (\dot \alpha + \dot \beta)^2 \cos(\alpha + \beta)\\ - a \dot{\alpha}^2 \sin \alpha - b (\dot \alpha + \dot \beta)^2 \sin(\alpha + \beta)\\ \end{array} \right) \end{array}\)

Step 1: Draw your multibody system¶

Draw the rod-crank-piston system in MBsysPad. Naming the joints is strongly recommended. Set both rod and crank rotations as dependent.

There is no specific action to do in MBsysPad for the user constraint.

Step 2: Generate your multibody equations¶

You have to generate the equations of motion of the system. You also need to regenerate the user models files.

Step 3: Write your user function¶

First, edit the user_JointForces function to add the external torque and the friction force to the system. A reminder about editing this function is available in the basic tutorial.
Edit the user_cons_hJ function:
- Enter the expression of the constraint: \(h_{usr}(q)\)
- Fill the jacobian matrix: \(J_{usr}(q)\)
- Include user_all_id.h.

//...
#include "user_all_id.h"
//...

void user_cons_hJ(double *h, double **Jac, MbsData *mbs_data, double tsim)
{
    //...

    // declare and define current variables
    double a_length = 0.15;
    double b_length = 0.3;
    double alpha = mbs_data->q[R3_rod_id];
    double beta = mbs_data->q[R3_crank_id];
    double x = mbs_data->q[T1_piston_id];

    // define the value of the constraints
    h[1] = a_length*cos(alpha) + b_length*cos(alpha+beta) - x;
    h[2] = a_length*sin(alpha) + b_length*sin(alpha+beta);

    // define the value of the jacobian matrix
    Jac[1][R3_rod_id] = -a_length*sin(alpha) - b_length * sin(alpha+beta);
    Jac[1][R3_crank_id] = -b_length * sin(alpha+beta);
    Jac[1][T1_piston_id] = -1;

    Jac[2][R3_rod_id] = a_length*cos(alpha) + b_length * cos(alpha+beta);
    Jac[2][R3_crank_id] = b_length * cos(alpha+beta);

    //...
}

Edit the user_cons_jdqd function:
- Fill the quadratic term: \(\left(\dot{J}\dot{q}(q, \dot{q})\right)_{usr}\)
- Include user_all_id.h.

//...
#include "user_all_id.h"
//...

void user_cons_jdqd(double *jdqd, MbsData *mbs_data, double tsim)
{
     //...

     // declare and define current variables
     double a_length = 0.15;
     double b_length = 0.3;
     double alpha = mbs_data->q[R3_rod_id];
     double beta = mbs_data->q[R3_crank_id];
     double alpha_dot = mbs_data->qd[R3_rod_id];
     double beta_dot = mbs_data->qd[R3_crank_id];

     // define the value of the quadratic term
     jdqd[1] = -a_length*pow(alpha_dot,2)*cos(alpha) - b_length*pow(alpha_dot+beta_dot,2)*cos(alpha+beta)
     jdqd[2] = -a_length*pow(alpha_dot,2)*sin(alpha) - b_length*pow(alpha_dot+beta_dot,2)*sin(alpha+beta)

     //...
}

Step 4: Run your simulation¶

Before running the simulation, you must modify the main script in order to indicate that there is a user constraint.

Open the main.c file and:

Include mbs_set.h.
add the mbs_set_nb_usrc after the program loading and before the partitioning in order to set the number of user constraints:

//...
#include "mbs_set.h"
//...
    //loading mbs_data

    // Set the number of user constraints
    int N_usr_c = 2;
    mbs_set_nb_userc(mbs_data, N_usr_c);

    //partitioning

// ...

The script must also call the coordinate partitioning before running the simulation. If it is not already there, add the instructions to call the coordinate partitioning module:

Open the main.c file and add the following command
- After the project loading (mbs_load() function)
- Before the time integration (mbs_run_dirdyn() function)

//...
// initialisation of structure for the coordinate partitioning
mbs_part = mbs_new_part(mbs_data);

// options for the coordinate partitioning
mbs_part->options->rowperm = 1;
mbs_part->options->verbose = 1;

// compute the coordinate partitioning
mbs_run_part(mbs_part, mbs_data);

// free the memory
mbs_delete_part(mbs_part);

//...

REMARK:

Some options are available for the coordinate partitioning process, see the documentation.

Check the results¶

Plot the graph of the displacement and the velocity of the piston and check your results with the following graph. Displacement and velocity of the piston