Advanced user constraints ------------------------- .. container:: c-code .. highlight:: c .. container:: c-code 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. \ :math:`\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: \ :math:`h_{usr}(q) = 0`\ - The user has to provide, for a given configuration: - The value of the constraint: \ :math:`h_{usr}(q)`\ - The Jacobian matrix of the constrain: \ :math:`J_{usr}(q)`\ - The quadratic term: \ :math:`\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. .. figure:: figure/bielleManivelle_image.png :alt: Rod-crank-piston system Rod-crank-piston system: definition The model is composed of the following elements: - Three bodies - the rod \ :math:`OA`\ - Length: a = 150 mm - Center of mass \ :math:`G^1`\ located in the middle of the segment \ :math:`OA`\ - Mass: 0.3 kg - Inertia along axis \ :math:`\hat{I}_{3}`\ : 0.005 kg.m2 - the crank \ :math:`AB`\ - Length: b = 300 mm - Center of mass \ :math:`G^2`\ located in the middle of the segment \ :math:`AB`\ - Mass: 0.7 kg - Inertia along axis \ :math:`\hat{I}_{3}`\ : 0.003 kg.m2 - the piston \ :math:`B`\ - Center of mass \ :math:`G^3`\ located in \ :math:`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 \ :math:`O`\ : \ :math:`T\hat{I}_{3}`\ with \ :math:`T=0.5\,Nm`\ - Friction force resistive to the piston motion: \ :math:`-c\dot{x}\hat{I}_{1}`\ with \ :math:`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: - \ :math:`h_1 \triangleq a~\cos \alpha + b~\cos \left( \alpha + \beta \right) - x = 0`\ - \ :math:`h_2 \triangleq a~\sin \alpha + b~\sin \left( \alpha + \beta \right) = 0`\ - \ :math:`\alpha`\ : rod rotation; - \ :math:`\beta`\ : crank rotation; - \ :math:`x`\ : piston rotation. The Jacobian is given by \ :math:`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 \ :math:`\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 \ :math:`\dot{J}\dot{q}`\ is therefore \ :math:`\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 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. container:: c-code - 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: \ :math:`h_{usr}(q)`\ - Fill the jacobian matrix: \ :math:`J_{usr}(q)`\ - Include user_all_id.h. .. code:: c //... #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: \ :math:`\left(\dot{J}\dot{q}(q, \dot{q})\right)_{usr}`\ - Include user_all_id.h. .. code:: c //... #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. .. container:: c-code 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: .. code:: c //... #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: .. container:: c-code - 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) .. code:: c //... // 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| .. |Displacement and velocity of the piston| image:: figure/disp_vel_piston.png