By Marco Morales
There are 6 configuration variables and they are the x,y, and theta for the cup and the x,y, and theta for the dice. The transformation matrix can be calculated from the world frame to the center of mass of the cup, which can be called g_WC by using x_cup,y_cup and theta_cup. Another transformation matrix can be calculated by using x_dice,y_dice and theta_dice and this is g_WD. The transformation matrix g_CD can be calculated by matrix multiplying the inverse of g_WD by g_WC. This will allow the dice frame to be described in the cup frame. The Kinetic energy of the system can be calculated by doing 0.5*(V_cup).T@I_cup@V_cup and adding it to 0.5*(V_dice).T@I_Dice@V_Dice which are the kinetic energies of the cup and dice respectively. The potential energy of the cup can be calculated by multiplying the mass of the cup by gravity and by the y component,which would be[1], of the product (g_WC@Matrix([0,0,0,1])). This is repeated for the potential energy of the dice excepting using g_WD instead. Subtracting the potential by the kinetic results in the Lagrangian. From here, performing d/dt(dL/dqot) - dL/dq will result in the Euler-Lagrange Equations. Since I want to make it easy to see impacts, I am applying a force in the positive y direction that cancels out the force due to gravity on the cup. Due to this, the force is acting in the y component with a magnitude of 9.81*mass of the cup. I am also applying a oscillating force in the x of the cup so that it mimics the cup shaking. There are 16 constraints in this system. Each corner of the dice has four constraints. Those being the top wall of the cup, the bottom wall of the cup, the left wall of the cup and right wall of the cup. More specifically, the x component of the dice is compared to the right and left walls and the y component of the dice is compared to the top and bottom walls. 2 conditions for each corner of the dice is the x component plus or minus the length of the cup and the other two are the y component of the dice plus or minus the width of the cup. These four conditions for each corner of the dice results in the 16 constraints. For the impact update, check the solutions calculated from sym.solve and use the solution where the lambda is greater than 0, or 1e-06. Then update the system with the new values for the velocities for x_cup,y_cup,theta_cup,x_dice,y_dice and theta_dice.
For my simulation, I applied forces in the y and x directions of the cup. This force was calculated in a way that it would cancel out the force of gravity acting on the cup and oscillate the system. Therefore, the whole simulation is a cup just moving in the x direction, at the start, and a dice inside the cup that is falling down. This would result in impacts as the cup is not moving vertically and the dice is falling down. As soon as the dice makes contact with the cup, the dice will bounce up and in turn will move the cup downward.
My simulation shows this exactly, as the simulation starts, the cup is moving only in the x direction and the dice starts to fall down. As the dice hits the bottom side of the cup, the dice bounces upward the cup starts to slowly move downward as some of the energy is being transferred to the cup that was not moving in the y axis at the start. Since the dice is continuously affected by gravity, it will keep going in the downward y direction and this results in the dice and cup system moving further down as time passes. The oscillating x force being applied on the cup results in the dice bouncing between the sides of the cup throughout the simulation.