moorepants / learn-multibody-dynamics

Interactive computational book on multibody dynamics
https://moorepants.github.io/learn-multibody-dynamics/
Other
121 stars 28 forks source link

Energy and Power Chapter #158

Closed moorepants closed 1 year ago

github-actions[bot] commented 1 year ago

Preview this pull request here https://moorepants.github.io/htmlpreview/?https://github.com/moorepants/learn-multibody-dynamics/blob/preview-pr-158/index.html

moorepants commented 1 year ago

Merging this, will take comments after posting. I didn't get everything I want in it but it gives enough to calculate energy.

moorepants commented 1 year ago

Great feedback @wwolfie I'll work on it.

Peter230655 commented 1 year ago

Above eq(218) you write: The rotational energy...and is defined as: (218) I believe, eq(218) expresses the total kinetic energy of this rigid body?

Peter230655 commented 1 year ago

In the paragraph above eq(214), it seems to me you leave the term 'energy' undefined. I guess, you have to, as there are so many forms of energy. Further down, you give examples of forms of energy. Would it make sense to add something like this: Energy has many forms, below specific examples are given ?

Peter230655 commented 1 year ago

Is the 'minus sign' in eq(219) just a convention, or is there more to it?

Peter230655 commented 1 year ago

Above eq(221) you write: A linear spring generates a conservative force F=kx between two points P and Q. Would this also be true, if k was a non linear term, depending only on the distance P-Q, und not on the specific way to get from P to Q?

Peter230655 commented 1 year ago

In fig 52, you label the unit vectors of N as n_1, n_2. In the code below, (of course) you orient frame A by rotating it around N.z. Would it be more consistent to label the unit vectors in fig 52 as n_x, n_y?

Peter230655 commented 1 year ago

Above the code where you define T_A you write: The torques on the thigh and calf will include a passive stiffness an damping to represent muscle tendons and tissue effects with coefficients kk and ck I think, this is a typo?

Peter230655 commented 1 year ago

Below thw graphs of the conservative motion you write: Checking whether energy remains constant is a useful for sussing out whether your model is likely valid It is probably fair to say, that every example I have played around with was set up wrongly initially! I now look at the total energy of every example I try yo set up.

Peter230655 commented 1 year ago

You write somewhere below Conservative Simulation with Ground Spring: This mismatch in the energies are due to numerical inaccuracies associated with the foot penetration not aligning precisely with the integrator’s time step. If you tighten the simulation tolerances and simulate a small enough time steps, the total energy should come closer to constant over the collision. This is the nature of numerical simulation for collision of very stiff systems

I had a lot of trouble with this issue on my 'bouncing ball': As a freely falling ball is simple to integrate, I guess solve_ivp will take large time steps. So, my ball simply fell trough the floor, until I forced small steps with max_step. ( and it took me quite a while to come up with this...) Should you show the improvement you will get with the measures you mention?

moorepants commented 1 year ago

Is the 'minus sign' in eq(219) just a convention, or is there more to it?

I copied from Kane's book. I think it must just be convention. Maybe I should remove the negative sign.

Would this also be true, if k was a non linear term, depending only on the distance P-Q, und not on the specific way to get from P to Q?

Yes, but a nonlinear term will give a different energy function. I just kept it simple.

It is probably fair to say, that every example I have played around with was set up wrongly initially! I now look at the total energy of every example I try yo set up.

If I couldn't only impart this wisdom to the students, but I think the only way to gain it is to struggle your self with simulations.

Should you show the improvement you will get with the measures you mention?

I spent a little time trying to improve the simulation with the settings of the IDA solver, but didn't get improvements in time. It needs a closer look. If you have some tips, I can update it. But now I have to move on to other things.

I fixed the other notes, but still need to update the figure with nx and ny.

Peter230655 commented 1 year ago

I spent a little time trying to improve the simulation with the settings of the IDA solver, but didn't get improvements in time. It needs a closer look. If you have some tips, I can update it. I have virtually no experience with the tae solver. What little I did, solve_ivp, with configuration constraints converted to speed constraint and with method='Radau' gave more reasonable results. Should I try to change your simulation to using solve_ivp?

moorepants commented 1 year ago

Should I try to change your simulation to using solve_ivp?

No I don't really want to do that. In the holonomic constraint sim chapter I basically teach them that you should use the DAE solver to get the correct solution. So I'd like to stick to that but it means we need to have a DAE solver that also behaves nicely for stiff systems. Here is the explanation of the solver we are using: https://computing.llnl.gov/projects/sundials/ida I probably need to give it a Jacobian and fiddle with the order and such.

moorepants commented 1 year ago

But you are welcome to show me what you get with solve_ivp and Radua, maybe you'll convince me to switch explanations.

Peter230655 commented 1 year ago

I will do it, hopefully over the weekend and let you know whatever I will have managed.

Only time I used dae was when , last year, you got to the chapter holonomic constraints. I 'complicated' your example a bit, and compared solver_ivp to dae (with IDA). Total energy was always 'more wrong' with dae.

Peter230655 commented 1 year ago

I may have found a small mistake in your simulation.

You set the spring force on the foot as Kf = kf*zp**(3/2)

But I think, then the spring energy should be Ef = 2/5 kf zp(5/2) (you seem to have set it as 0.5 kf zp2)

NB: Even this correction does not seem to make solve_ivp work, but I have not given up, yet.

Peter230655 commented 1 year ago

maybe another small mistake?

In the simulation you set T_A = kk(q3 - pi/2) Should the energy then not be E_A = 0.5 kk * (q3 - p/2)2 - (pi/2)2 ?

moorepants commented 1 year ago

Yes, thanks. Those are errors. Will fix.

Peter230655 commented 1 year ago

Would it make sense to mention in the code that the first equation of eq(216) was used to get these terms?

moorepants commented 1 year ago

I could write out the integrals write before I introduce them.

moorepants commented 1 year ago

Or solve the integrals with sympy.

Peter230655 commented 1 year ago

Solving them with sympy looks more attractive to me - I did it the old-fashioned way.

moorepants commented 1 year ago

Fixing the energy equations fixed things:

image

But I wasn't quite sure why I needed the negative sign on the foot potential energy.

Peter230655 commented 1 year ago

I still play around with solve_ivp, not too successful so far. I am unsure, I understand what you mean by the minus sign. You subtract the spring energy of the foot? You still use kf = 1.e7, kk = 10?

moorepants commented 1 year ago

I'm not subtracting the energy, but the integral gives a negative sign that isn't needed.

Peter230655 commented 1 year ago

So, sympy's symbolic integration is wrong? I have never tried this. Still messing around with solve_ivp :-)

Peter230655 commented 1 year ago

I think, I give up on solve_ivp for your jumper. While the configuration constraint is kept pretty well (order of magnitude = 1.e-8) the error in the speed constraint is in order of the speed. So, I guess the kinetic energy cannot be right , and the total energy is far from being constant. Even to get this, I had to set atol = rtol = 1.e-12

In order to be in 'my world', I set up the equations of motion using sympy mechanics.

If I calculated more than about 0.3 sec, I had to set kk = 1000 ( you have it at kk = 10), else the jumper would 'flip over', makes sense!

--> When I retire and get my desk top computer, I have to learn about dae!