Two springs & weights "piggy backed"

[davey110]

This simulation of 2 springs & weights “piggy-backed) together was suggested by Henko on Dutchman’s “Mass on Spring Resonance”. So I took up the “challenge”. It simulates the oscillations produced after pulling down & releasing the lower weight. I don’t know how accurate it is, but should give a general idea. It works by calculating the movement of the weights over a very small distance then stringing these distances together to see what happens. The tension in each spring is dependent on the other, making it almost impossible to program in small increments without any error. You would get a circular calculation. Therefore this simulation is probably not too accurate, although I think the general idea is correct.

I think it's because of this small error that the oscillations increase until both weights eventually fly off the top of the screen. I have simulated some "air resistance" to the velocities, which prevents this "runaway". With this friction, the system should gradually come to rest, but it appears that it will run forever. Possibly there's an error in the code that I have missed. Or maybe this is the best that can be done with this method. (I also wrote a program with only 1 spring using the same method, which seems to work perfectly.)

The time and distance units are arbitrary. Making the accelerations (a1 & a2) smaller simulates a shorter time period for variable t, which may increase the accuracy. The main variables can be adjusted near the top of the listing- weights (mass), spring constants, length of springs, etc

I have left in plenty of REM statements, which should help in following the logic. Delete any that you find too trivial.

Code: Select all

/* 
Suggestion ON FORUM: motion of 2 springs piggy-backed.

I don't guarantee the accuracy of this program, although it gives a reaonable demonstation.
Because the actions of the upper & lower springs both depend on each other, it is almost impossible to calculate using this method without calculus. So there is always a small error. The weights m1 and m2 eventually fly off the top of the screen. to preent this I have added a small amount of "air resistance" to contol increasing velocities. 

It also plots a graph of the movements while it is running. This can be REMed out if not wanted. 

All the constants can be adjusted as desired, to see the effect.

It does not give any spectacular results, but settles down to a steady pattern.

Assumptions:
1. Length of springs is calulated to CENTRE of weights.
2. Springs are always within their range of elasticity. 

NOTES:
1. The tensions in both springs are dependent on each other, which is almost impossible to program in small increments. You would get a circular calculation. Therefore this simulation is probably not too accurate, although I think the general idea is correct. 

2. The time and distance units are arbitrary. Making the acceleration smaller would simulate a shorter time period for variable t, which should increase the accuracy. The main variables can be adjusted near the top of the listing. 

3. I think it's because of this small error that the oscillations increase until both weights eventually fly off the top of the screen. I have simulated some "air resistance" to the velocities, which prevents this "runaway". With this friction, the system should gradually come to rest, but it appears that it will run forever. Or possibly there's an error in the code that I have missed. Or maybe this is the best that can be done with this method. [Programs "Spring-1- ONLY.sb" seems to work perfectly. 

VARIABLES: naming system:
Upper spring & weight are #1, bottom are #2. 
k       is Hooke's Law constant:  Force=k*change in length of a spring
all pxx  are positions - only used for placing of Sprites, not in calculations
psr1 -  relaxed posn of spring 1 (pixels from top of screen)
pb02 - initial hanging posn of ball 2
Lr1 -    relaxed length of spring 1
L02 -     initial sretched length of spring1 i.e. hanging stationary with weights (balls) 
L1, L2 -   current length of springs; only used to calc ps1 & 2 and pb1 & 2. 
m1, m2 -  mass (weight) of balls
ts1, ts2 -  tension in springs; ALWAYS a +ve number; direction depends if it's pulling up or  
             '  down on a ball. It must be added or subtracted accordingly. (UP is -ve, ie subtract). 
f1, f2 -   NET FORCE (sum of all forces) pulling on balls 1 & 2 (including weights); zero when stationary 
                'limit  a1 & a2 to prevent runaway
fh  -      "human" force pulling on ball 1 to start it off. 
*/

OPTION BASE 1
OPTION SPRITE POS CENTRAL
GRAPHICS ! GRAPHICS CLEAR 1,1,1
DRAW SIZE 2
DRAW COLOR 0,0,0
FILL COLOR 1,0,0
 
'g'  ===========  INITIAL SETTINGS, CONSTANTS =====================
'. ======== 1 IS TOP,  2 IS BOTTOM ===== '
    'HOOKES LAW: F=K*∆D === bigger K is stronger spring ========
    'relaxed posn's ie no tension in springs: psr1,2; f1, 2=0, ==================

m1=10  'mass (weight) of balls
m2=10   
k1=0.1   'spring constant: tension/stretch  (mass/pixels) (Hooke's law: F=k*∆L)
k2=0.1   
Lr1=100 
Lr2=100     'relaxed length, i.e. lying flat on a table, no forces in effect
fh=10       'human pull applied
speed=.01   'delay for control of display speed. 
RANDOMIZE
fh=20       ' "human" pull applied to start it off. 
r1=RND(.9)
fh=fh*r1+3  'Random variation in human starting pull. REM out if not wanted.
' dbp        DEBUG PAUSE.   See DEF at end of listing

'r'    
SPRITE "s2" BEGIN 20, Lr2    ' Spring 2
   FOR y=1 TO Lr2 STEP 10
     DRAW LINE 1, y TO 20, y+5     ' down to right
     DRAW LINE 1, y+10 TO 20, y+5  ' up to right.            
   NEXT y
SPRITE END
SPRITE "s2" AT 500, Lr1+Lr2/2 ! SPRITE "s2" SHOW
'-------------
SPRITE "s1" BEGIN 20, Lr1    ' Spring 1 
   FOR y=1 TO Lr1 STEP 10
     DRAW LINE 1, y TO 20, y+5     
     DRAW LINE 1, y+10 TO 20, y+5  
   NEXT y
SPRITE END
SPRITE "s1" AT 500,Lr1/2 ! SPRITE "s1" SHOW
'-------------
SPRITE "m1" BEGIN 60,60   'Mass 1  
  FILL CIRCLE 30, 30 SIZE 30
SPRITE END 
SPRITE "m1" AT 500, Lr1 ! SPRITE "m1" SHOW     'relaxed posn
'--------
SPRITE "m2" BEGIN 60,60   'Mass 2
  FILL CIRCLE 30,30 SIZE 30
SPRITE END 
SPRITE "m2" AT 500, Lr1+Lr2 ! SPRITE "m2" SHOW   'relaxed posn

201 
'' ============== Inital Positions, Relaxed on table =========
'       ========== As set above =======
DRAW TEXT "Relaxed " AT 200, Lr1+Lr2
PAUSE 2

202
'g' ======== start with INITAL STRETCH, stationary: L01, L02,ps01, pb01 etc ==========
        ' =======CALCULATE SPRING 2 FIRST because it affects Spring 1 ==========
        
ts2=m2    'tension/forces on springs, hanging stationary
f2=m2-ts2 'NET forces=0 on ball: UP = - and DOWN = + (because screen counts top down
          'net force on ball=0; DOWN-UP; Tension is always a +ve number
strch2=m2/k2  'adds to Lr2
L02=Lr2+strch2 ' stretch;  f= k8d,------- L=L+d ---------  

 ' ======== NOW DO SPRING 1 & THEN PLACE SPRING 2 & m2 SUSPENDED FROM IT ===== 
102
'g'start with INITIAL STRETCH, stationary: L01, L02,ps01, pb01 etc ==============

ts1=m1+m2          'tension/forces on springs, hanging stationary
f1=m1+ts2-ts1      'net force=0 on ball: UP is - and DOWN is + (because screen counts top down
                        'net force on balls;  =0 when stationary  
strch1=(m1+m2)/k1  'stretch of spring; adds to current length
L01=Lr1+strch1     ' f=k*d, so L=L+d ---------   

ps01=L01/2        'spring pos'n is CENTRE of spring
ps02=L01+L02/2  
pb01=L01
pb02=L01+L02

SPRITE "s1" RESIZE 20, L01 'resize sprites to hanging (with weights), stationary size
                   ‘ RESIZE is a permanent resizing. Don’t use in a loop 
SPRITE "s1" AT 500, ps01 
SPRITE "s2" RESIZE 20, L02.  'stretch spring to initial length
SPRITE "s2" AT 500, ps02  
SPRITE "m1" AT 500, pb01  
SPRITE "m2" AT 500, pb02 

GRAPHICS CLEAR 1,1,1 
DRAW TEXT "Hanging, stationary" AT 200,pb02  ‘ so you can see what’s happening
PAUSE 1    ‘ so you can read it  
L1=L01  !  L2=L02 'as well as ts1, ts2 above:  use in next section (both f1, f2=0)

203
'c'
' add (human) pull on bottom weight m2 =================
' initial human pull on ball added
        ' calc new t=1 posns etc from t=0 data
ts2=m2+fh 
f2=m2+fh-ts2    ' now includes fh    ' THESE ALL NEEDED IN NEXT SECTION
ts1=m1+ts2      ' (or ts1=m1+m2+fh since ts2=m2+fh)
f1=m1+m2+fh-ts1 '(or f1=m1+ts2-ts1)
L2=L02+fh/k2     'new L2, L1  both increased by fh/k. This follows Hookes Law f=k*d,  d=f/k  
L1=L01+fh/k1

ps1=L1/2  !   ps2=L1+L2/2.   'new positions
pb1=L1    !   pb2=L1+L2

 SPRITE "s1" AT 500, ps1 SCALE L1/L01   ‘ scales on screen to current length 
SPRITE "s2" AT 500, ps2 SCALE L2/L02    'new L2.  
SPRITE "m2" AT 500, pb2 
SPRITE "m1" AT 500, pb1
'dbp
GRAPHICS CLEAR 1,1,1 
DRAW TEXT "Added pull" AT 200, L1+L2
PAUSE 1 'before starting main  
                                      
'======================= NOW THE FUN PART - Make it all move ==================
204
'b' now REMOVE human pull fh; calculate a, d, new ts === MAIN PART OF PROGRAM ===
GRAPHICS CLEAR 1,1,1 
DRAW TEXT "Pull released; now running" AT 150, L1+L2
PAUSE 1
'start with  current values at t=0
' calculate new values at t=1 
   ' fh=0        Just a reminder
f2=m2-ts2       ' ts2 (from above) includes fh, so f2 now does as well. 

 'These 3 equations are used: (1)a=f/m  (2)d=(a*t^2)/2 where initial t=1. (3) v=a*t
 'Also Hookes Law  f=k*d,  [ d=f/k ]
 'All positions AND TENSIONS are still unchanged UNTIL m1 STARTS TO MOVE! ' f2 from previous
 'NET FORCE on m2 will be UP (-ve) at start (t=0), so it will accelerate up
 'Tensions are POSITIVE #s. Direction depends what it's pulling on. 
'dbp
FOR  t=1 TO 30000  '=================== MAIN LOOP =================
    ' f2 now begins to accelerate m2 upwards.
   a2=f2/m2     ' since ∆t=1    for t=1, eqn (1) above, will be NEG UP
   v2=v2+a2     ' since ∆t=1;   v2+(eqn 3)
  v2=v2*(1-.0001*v2^2)   'a small amount of "air friction" to prevent runaway.
   d2=a2/2    ' distance moved  d=v*t; time interval ∆t=1. eqn (2)
   ts2=ts2+k2*d2 
   IF ts1<0 THEN ts1=0
   L2=L2+d2     ' d2 will be -ve (UP) AT t=1
' ----------------- now d1 ----
  f2=m2-ts2     ' new net f2 - reduced by ∆L2.  F=K*D; d2 is -ve i,e. UP
  f2=f2+k*d2
          ' ======= NOW DO m1,L1 etc, including effect of m1 ========
   f1=m1+ts2-ts1  ' still using old ts2
   a1=f1/(m1+ts2)
   v1=(v1+a1)
  v1=v1*(1-.00005*v1^2)   'a small amount of "air friction" to prevent runaway.
  ‘ d1=v1/2
   d1=a1/2
   ts1=ts1+k1*d1
   IF ts1<0 THEN ts1=0  ' spring doesn't compress
   L1=L1+d1       'at t=1, d1 is -ve 

   ps1=L1/2     
   ps2=L1+L2/2      'new curr posns. 
   pb1=L1
   pb2=L1+L2
   
   SPRITE "s1" AT 500, ps1 SCALE L1/L01  'scale sprite for display
   SPRITE "s2" AT 500, ps2 SCALE L2/L02
   SPRITE "m1" AT 500, pb1
   SPRITE "m2" AT 500, pb2

GOSUB prgraph  'draws graph, REM out this line if you dont want it. 

'dbp
'PAUSE 1
   PAUSE speed  'display speed - set at top of listing. 
NEXT t       '============ END OF MAIN LOOP ================== 

DRAW TEXT "========= DONE =========" AT 100,550
END          

'y'===================================================
DEF dbp ! DEBUG PAUSE ! ENDDEF     'very handy for debugging
'-------------------------------------------------
prgraph:          'optional subroutine - for printing graph of motion
'DRAW SIZE 1
 IF t/100=INTEG(t/100) THEN DRAW LINE t, 150 TO t, 600 
 FILL COLOR 0,0,0
 FILL CIRCLE t/2,pb1 SIZE 2       'plot graph points
 FILL CIRCLE t/2,pb2 SIZE 2
RETURN  

If any errors or typos, please let me know.

[Henko]

Hi Davey,

You did quite a job on the “two body/two spring” object.
I had a look at your method, the math, and the code, and have some remarks/sugestions.

If the kinematics are done correctly, the system should remain within its limits, also without the correction due to air drag.

Most important is an erronous formula 2, to calculate the distance d, which an object travels in a given time interval dt (one second in the program).
The formula should read : d=v*dt+0.5*a*dt^2 or d=v+a/2 if dt=1 (instead of d=a/2)
Hence the distances as calculated in the program are far to low.

The program deserves that the time step is not fixed upon 1 second, but be flexible.
All that is necessary is:
- define dt=<value of time step> somewhere in the beginning of the code (may be dt=1)
- use the timestep in the formulas where appropiate.

I suggest the following calculation model for each time step:

a = F/m (exclusive weight and air drag)
dV = a*dt
Vnew = Vold + dV
Vaverage = (Vold + Vnew)/2
dS = Vaverage * dt

And the new status of the object is Vnew for speed and a new position augmented with dS.

Other suggestion: you could incorporate the effect of weight and air drag as optional.
Both are forces, hence they must be added to the first formula, which would take the form:

F = <calculation of string force>
If weight=1 then F += m*g ( g being the gravitational constant, 9.81 on earth)
If air_drag=1 then F += c*V^2 ( c being a drag coefficient)
a = F/m
( += or -= depending on positve direction of forces).

I hope that the remark are usefull, and that you can have a go with it.

[davey110]

Thank you for your detailed comments. They are very helpful. I have to look at them in detail. At 1st glance I’m not sure if I have explained my method clearly.

1) The velocity & acceleration are per “t”, not per actual second. At least that’s what I intended. Therefore d=.5*a*t^2 reduces to d=a/2, since t=1. So if you reduce the acceleration (a) then (t) is effectively representing a smaller time interval, whatever it may be. “t” itself is not defined except as the time interval we are looking at. But after reading your comments, I now realize that the only way to change the acceleration is to put a “scaling factor” in a=F/m which I think would equate to a shorter t. I missed that, but I’m not sure if it’s the right answer.
[While working on this program I wrote one with only 1 spring & weight. It works perfectly & doesn’t run away, so I thought I must have it right. So I duplicated it & basically just “hung” it under the 1st one, but also had to be integrated into the movements of the upper one.]

2) Trying to keep things simple, I have assumed that weight = mass, and that mass is the downward force regardless of g. As a result many of the values are totally fictitious. e.g. the spring constant k.

3) The incorrect formula is REMed out (unless I missed one). But I now see a problem with d=v/2; I’ll have to look at that again. It could me my error. I think it should be something containing d=∆v/2, since v is not constant over the time interval. I think that’s what I was aiming for but didn’t quite get it.

4). The air drag can be REMed out to see the runaway. It is square law so should make sense. I used a graphing calculator app to get the response I was looking for, for + and - velocities. It’s not necessarily realistic, but stopped the runaway.

And that’s as close as I’ve been able to get it. Hopefully your comments will help me improve it. Especially to stop the runaway without resorting to air drag. Thank you.

[davey110]

As a start, I added an “acceleration factor”, af, making the equation a=af*f/m in both places where the equation is used. With values of around af=0.1 and also REMing out the air drag in both places [v=v*(1-.0001*v^2)], the whole thing plays out with a very nice looking resonance. However, I still have to study your suggestions to see if this pattern is realistic. That will take a bit more time.

[davey110]

As a start, I added an “acceleration factor”, af, making the equation a=af*f/m in both places where the equation is used. With values of around af=0.1 and also REMing out the air drag in both places [v=v*(1-.0001*v^2)], the whole thing plays out with a very nice looking resonance. However, I still have to study your suggestions to see if this pattern is realistic. That will take a bit more time.