disc-on-very-uneven-street

[moorepants]

From Peter:

I try to simulate a disc rolling on an uneven street without slipping. The street may be such that second contact points (CP2) happen. I 'numerically' look for possible CP2 at every step of the numerical integration. If there is one, I interrupt the integration, using solve_ivp's events keyword. Then I restart the integration with new initial conditions, particularly the new starting contact point is CP2.

I found this:

1. It is quite sensitive to the shape of the street. I am not sure, why, as the message of solve_ivp does not tell anything useful, for me anyway.

2. In order to find CP2 accurately, the max_step in solve_ivp must be set to small values, slowing down the integration. (As the system appears so well behaved to solve_ivp, no stiffness issues, I guess it selects large steps, hence the difference in integration time is even more visible)

3. The total energy jumps at the discontinuities. I understand this for the potential energy, if CP1 and CP2 are at different levels, but it is the kinetic energy which jumps much more. I do not know, which physical principle causes this to happen. (or did I make a big plunder in the description of my system.) It is not some mistake in calculating the total energy, the simulation also confirms these jumps.

I finally decided to adjust the rotational angular speed so that the total energy before and after the discontinuity are the same. (In between discontinuities, when the disc simply rolls on the street, the total energy is constant anyway)

If you want to run the simulation, it only runs about 30 sec on my iPad, but then the animation may take another minute.

[Peter230655]

Points 1 & 2 have been take care of, now that I use the keyword events in solve_ivp correctly. Soem helpful person on Stack Overflow explained it to me.

Point 3 is still unclear to me from a "physical" point of view.

[Peter230655]

Point 3: The explanation is this: The moment of inertia w.r.t. old contact point in general is different from the moment of inertia w.r.t. new contact point. As the kinetic energy must remain constant, this forces the angular speed to change discontinuously. Sorry it took me so long to find something this obvious!

[Peter230655]

Point 1. I have a better version, which does not exhibit this problem: Only difference is that I now apply the keyword events in solve_ivp correctly. The function event ( so that events=event) in solve_ivp) must return a negative real value, I used -1. if the event triggering event has not happened and a positive real value, I used +1. if it has triggered. ( a guy on stack overflow explained this to me )