Search:
 

 
Examples:
 

 

Example - Pendulum

* The "\" symbol has been used at the end of lines to indicate that the line was too long for the display and has been
wrapped to the line below in order to fit within the width of the web page. These lines would need to be joined in the normal
program editor.

I. Introduction

This simulation demonstrates the dynamics of a simple pendulum. A periodic driving force, or a periodic shift in the centre of rotation can be added and the effects examined. A finite element array is used to model the situation.

II. Mathematics and Physics Involved

The equation of motion for a simple, undriven pendulum with damping proportional to velocity is where f is the angle (in radians) of the pendulum from the vertical, γ is the damping constant, and , where g is the acceleration due to gravity (9.81 ms -2 ) and l is the length of the pendulum.

Note that the mass of the pendulum is not relevant; Galileo famously discovered this result.

When the pendulum is driven by a force, the equation of motion changes to account for this. There are three examples of driving in the program: a periodic force, and force resulting from the motion of the suspension point horizontally and vertically. The equations of motion for each are

respectively, where a and f are the amplitude and frequency in each case.

When the frequency of the force is close to that of the pendulum, a phenomenon called resonance occurs: the angle the pendulum swings increases each time it is ‘hit' by the force. Note that frequency is much more important than amplitude.

III. How the Program Works

1) Initialisation of the Finite Element Array

To start with, a finite element array is set up to describe the behaviour of the undriven pendulum. The formula for the second differential of f is put into formu$ , and then the SetFormula() method is used to assign it to the finite element array.

     GLOBAL formu$
     
     'set up fea
     IF EXISTS ("Pendulum") THEN Pendulum.Delete()
     
     GLOBAL fea AS FiniteElmtArray
     fea.Name = "Pendulum"
     fea.Dimensions(1) = 1
     DIM fa AS FiniteElmtAttribute
     SET fa = fea.Add("phi")
     
     formu$ = "- (gamma% * fea.phi(1) + omega% ^ 2 * SIN (fea.phi\
     ()))"'basic formula (undriven)
     fa.Order = 2
     fa.SetFormula(formu$,2)
     fa.AllowableErrorFactor = 0.01
     fa.AllowableErrorMin = 0.01

2) Update of the Parameters

When the Start button is clicked, the values of the parameters are changed to reflect the new values in the text boxes. In addition, if any of the driving force options is on then the finite element formula is altered. 

SUB startproc()
   'initialise pendulum (eg change f.e.a.) before it is set in motion
   IF fea.Running = - 1 THEN fea.Stop()
   l% = VAL (f.lengthbox.Text)
   gamma% = VAL (f.dampbox.Text)
   IF f.damplab.Value = 0 THEN gamma% = 0
   omega% = SQR (9.81 / l%)
   
   IF f.none.Value = - 1 THEN 
     fea.phi.SetFormula(formu$,2)
   ELSE 
     A% = VAL (f.drampbox.Text)
     freq% = VAL (f.drfreqbox.Text)
     IF f.per.Value = - 1 THEN 
       formu2$ = "+(A%/l%)*omega%^2*cos(2*pi*freq%*fea.t)"'term \
       for driving force
       fea.phi.SetFormula(formu$ + formu2$,2)
     END IF 
     IF f.hor.Value = - 1 THEN 
       formu2$ = "+A%*cos(2*pi*freq%*fea.t)*cos(phi())"'term for \ 
       horizontal motion of pivot
       fea.phi.SetFormula(formu$ + formu2$,2)
     END IF 
     IF f.ver.Value = - 1 THEN 
       formu2$ = "-A%*cos(2*pi*freq%*fea.t)*sin(phi())"'term for \
       vertical motion of pivot
       fea.phi.SetFormula(formu$ + formu2$,2)
     END IF 
   END IF 
   
   fea.phi.SetInitialValue( PI / 4,0)
   fea.phi.SetInitialValue(0,1)'starts off stationary
   fea.AdvanceTo(0)
   f.Tag = "run"
 END SUB

3) The runproc() Procedure

Now, as long as the Stop or Quit buttons are not clicked, the program keeps executing the procedure runproc() . In this procedure the finite element array is advanced by 0.05 seconds, and the procedure drawpend() is called, which draws the current position of the pendulum, taking f as the argument. 

SUB runproc()
   DIM tim%'current time
   DIM tt%'time to which to advance fea
   tim% = Superbase.Now
   tt% = fea.t + 0.05
   fea.AdvanceTo(tt%)
   CALL drawpend(fea.phi())
   tim1% = Superbase.Now
   IF (tim1% - tim%) / 100 < 0.05 THEN 
     WAIT FOR 0.05 - (tim1% - tim%) / 100'in case program executing too \
     quickly
   END IF
 END SUB

4) Output

The following picture shows the program while it is running. The pendulum is undamped and undriven, which results in its swinging normally.

If, however, the driving is turned on, quite different results are obtained. The driving force (represented by the red line coming from the bob of the pendulum) can make the pendulum behave in a seemingly random fashion; it will sometimes make it spin ‘over the top' and sometimes swing in the usual fashion.