Model questions and clarifications #1

[jonathanolson]

After reviewing some of the model, I want to ask for a few clarifications (including documentation). If something in here is correct and not documented in the model, it should ideally be added (it will help future developers).

It looks like height/velocity (on the Pendulum) isn't used (see https://github.com/phetsims/pendulum-lab/issues/31). 'acceleration' is documented as "acceleration value of pendulum", but it seems to be angular acceleration. I'd recommend renaming to angular acceleration (or alpha, to match omega as angular velocity).

The step first uses the typical 1D motion formula to determine the new angle, and then computes angular acceleration for the next step based on the new angle.

Friction looks almost correct, but I believe it should be proportional to the velocity squared (not just the velocity). Drag force is proportional to Av^2 (A: surface area, v: velocity). Assuming constant density (and a sphere to simplify), A is proportional to r^2 (radius), and r is proportional to m^(1/3) since volume will be proportional to mass. That means the drag force is proportional to v^2m^(2/3), and with F=ma we have our acceleration proportional to v^2/m^(1/3). The current sim model is using v/m^(1/3).

Then angular velocity is stepped forward with an average of the old and new acceleration values. Does this help ensure the accuracy of the simulation in some way?

I'm unsure of Pendulum's updateAccelerationVector. Since the magnitude of the net force on a pendulum is -mgsin(theta) (see the "Force" derivation for http://en.wikipedia.org/wiki/Pendulum_%28mathematics%29#math_Eq._1), the magnitude of acceleration is -gsin(theta), which in our model is (angular acceleration)length. We seem to be adding an additional factor to the magnitude by using the formula length*sqrt( (angular acceleration)^2 + (angular velocity)^2 ), and I can't find justification for it (here's the code):

accelerationMagnitude = this.length * Math.sqrt( this.acceleration * this.acceleration + omegaSq * omegaSq )

Also notably, the result is the magnitude of the acceleration vector, while the input 'acceleration' is the angular acceleration.

For understanding, it would also help to have a function like the following in Pendulum:

function setAcceleration( frictionContribution ) {
  this.acceleration = -this._gravityProperty.value / this.length * Math.sin( this.angle ) + frictionContribution;
}

so that this formula (which is essentially copied with different ways of referencing gravity/length/angle) from both Pendulum.js and PendulumLabModel.js.

If I made mistakes above (or have conceptual misunderstandings), please let me know. Additionally, I appreciate the variable-timestep model from the start!

[andrey-zelenkov]

Renamed Pendulum.acceleration to Pendulum.alpha and change corresponding comments. Also added Pendulum.setAlpha function.

Сoefficient of power of the velocity for drag force depends on Reynold's number (http://courses.ncssm.edu/aphys/labs/vlab1/theory.pdf) and can vary from 1 to 2. For laminar flow (not very high velocity) usually takes F∝v^1. Friction formula was taken from the original simulation and there was proportionality to the first power of velocity, but if necessary we can change it to the second. Original code (_global.drag is value set by friction slider):

this.alpha = -_global.g / this.length * Math.sin(this.theta) - _global.drag / this.mPowThird * this.omega;

Formulas for calculating parameters of motion of the pendulum and friction were taken from the original simulation. Please let us now if anything should be changed. Here is code for calculating acceleration value:

if (this.accelerationShown) {
    var accArrow = this.clip.pendulum_mc.vectorHolder_mc.accArrow_mc;
    var alpha = this.pendulum.alpha;
    var omegaSq = omega * omega;
    var alphaMag = R * Math.sqrt(alpha * alpha + omegaSq * omegaSq);
    var angleInDeg = Math.atan2(omegaSq, alpha) * 180 / 3.141593;
    accArrow._xscale = accArrow._yscale = 20 * alphaMag;
    accArrow._rotation = -angleInDeg;
} // end if

[jonathanolson]

The mass part of the drag computation was actually correct, since I was noting the force, and the old code was noting acceleration (div by m^1/3 => multiply by m^2/3).

http://fy.chalmers.se/~f7xiz/TIF081C/pendulum.pdf notes that for many pendulums, the Reynolds number is sufficient to have drag be proportional to the square of the velocity, so I've implemented that. It's quite noticeable, so if it's undesirable I can revert it back to the small-angle (or low Reynolds number) approximation of proportional to the velocity.

Model should now be pretty robust, and accepts variable time-steps so we avoid a number of issues.

[veillette]

reopening this issue. There is a problem with having solely a quadratic drag. At low velocity, the quadratic drag becomes very ineffective at dissipating energy, (since the velocity becomes so low, the square of it becomes even smaller). For instance, a cart given an initial push can still reach infinity, if it is only subject to quadratic drag.

In the current simulation, even at the highest friction level, the pendulum keep swinging even after a few minutes. The amplitude of oscillation is inversely proportional to time whereas the amplitude should be exponentially decaying with time.

In real life the linear and quadratic drag are both present, and the quadratic term dominates at high speed, and the linear term dominates a low speed. The linear term is responsible for bringing it to a halt in an exponential fashion.

The paper quoted above includes both type of drags and it should be straightforward to include both type of frictions (both would have coefficients that would be proportional to the value from the drag slider). Any thoughts? @jonathanolson

[jonathanolson]

I definitely don't have objections to including both a linear and quadratic term.

[veillette]

The linear and quadratic drag was implemented. Closing