Trebuchet simulation&some newbe question
Author Message Trebuchet simulation&some newbe question

The July issue of "Scientific American" describes the Trebuchet,
an ancient  (500 AD) war machine capable of hurling massive
projectiles an astonishing distance.

The magazine cover shows a car flying through the air after being
launched by a modern reconstruction of the machine. The article has
pictures of an upright piano subjected to a similar indignity.

The authors mention that a computer program was used to explore
performance parameters of the machine so I resolved to see what I
could do with the simulation problem.

As might be expected for such an ancient war engine, the basic
idea is quite simple, but I found that there are some subtle
nuances in the design when at last I got my program to run and
produce what I think are reasonable results.

The Trebuchet is an unequal arm lever pivoted on an elevated
bearing so the arm is free to swing in a vertical plane. A heavy
counterweight is fixed to the short end of the lever and a sling
(the David & Goliath type) is attached to the long end. The
mechanism is prepared to launch a projectile by pulling down the
long end of the lever, placing a projectile in the sling pocket and
hooking one of the sling ropes to a hook that is designed to
release when the arm reaches the proper angle. The other sling rope
remains attached to the arm.

My simulation revealed that designing this hook might be a mite
tricky along with fixing the relation between lengths of the sling
and the long arm of the lever.

After pursuing many dead end paths, I realized that the machine is
a compound pendulum where the overall system (arm, counterweight,
and projectile) rotates about the pivot with an angular
acceleration equal to the unbalanced torque divided by the moment
of inertia of the system.

The sling in turn rotates about the long end of the lever arm with
an angular acceleration given by a similar equation but in this
case the moment of inertia is constant; the mass (in slugs) of the
projectile times the square of the length of the sling rope.

To solve the problem first and second degree differential
equations relative to  both the arm and sling must be integrated to
obtain the angle of the lever from its acceleration and the angle
of the sling from a force balance on the projectile. I maintained
that simple Euler integration should be adequate because, given
small time increments in difference equations, all the derivatives
would change smoothly without inconvenient discontinuities. A
friend suggested that I use Runge-Kutta integration because the
double integration might introduce excessive error propagation.
My latest version of the simulation makes use of Adams predictor-
corrector integration which, I believe, produces a better result
but at the expense of complexity that I prefer to avoid here.

This brings me to my 'newbe' question about integration using J.
In fact I would appreciate it if someone could show me how to solve
the simulation in J. I have spent a good many hours trying the
iteration function f^:(n) g for integration subroutines, a loop
with label_loop. at the head of a list of J statements and
if. (200 >.N=.N+1) goto_loop.
end.
)
I'm missing something. It seems to me this should work but it
doesn't. I keep getting errors even when I simplify to the point of
just incrementing the counter, N. I'm using J for windows V2.05.

Here is my simulation program:
Please forgive me I'm using 'pigeon J' in order to transmit this
over the net so the APL program won't run until J symbols are
translated to corresponding APL characters. I have provided a
translation table at the end of the nomenclature list for those who
remain confirmed APL users. I have added one symbol :: to represent
the diamond APL symbol used to separate independent statements on
the same line.

Nomenclature list:
W1 = weight of counterweight (lbs)
m1 = mass of counterweight (slugs = W1% 32.2)
W2 = weight of projectile; m2 = W2%32.2
L1 = length short end of lever arm (ft); L2 = length of long end
SL = length of sling (ft); P = X-Y coordinates long end of lever
T = coordinates of projectile relative to P
S = coordinates of projectile  (S=P+T )
dT = integration time step
a = (dV/dT); V = (dS/dT) projectile acceleration & velocity
wp (for w-prime) the angular acceleration of the lever arm
w = angular velocity of the lever arm
wp2 = angular acceleration of sling and projectile
w2  = angular velocity of sling and projectile
Translation table:
% = divide  * = multiply o. = o (APL circle function)
=. is APL assignment arrow  _ is APL negation function
^ = * APL power function    \$ is rho, the APL shape fuunction
NB. = APL lamp symbol (comment)

Simulation program with typical parameter values:

 {del.} ZOT    (credits to the anteater in the comic strip BC)
 W1 =. 2000 :: W2 =. 20 :: L1=. 6 :: L2=. 30
 m1=.W1%32.2 :: m2=.W2%32.2 :: N=.1 :: dT=.0.01 :: w=.w2=.0
 A=. o.%6:: P=. L2*(1 _1)*1 2o.A:: T=. _8 _6 :: SL=.(+/T*T)^0.5
 A2 =. _1oT%SL :: S=. P+T:: a=.V=.0 0::
 C=. m1*L1*L1 :: I2=. m2*SL*SL
 TB=.1 17\$ 0,P,S,A,w,wp,A2,w2,wp2,V,a,T

 Ln0: I=.(C +m2*(s=.+/S*S)) :: s=. s^0.5 NB. Moment of inertia
 wp=. ((W1*L1*1oA) -W2*S) % I NB.System angular acceleration
 w =. w +dT*wp           NB. angular velocity
 A =. A +dT*w            NB. new lever arm angle
 P =. L2*(1 _1)*1 2 o A  NB. end of lever arm
 Z =. (w,wp) o.* s * 2 1 o A NB. angular velocity, acceleration
times vector s from origin to S
 V =.Z[1;] :: a =. Z[2;] NB.projectile velocity & acceleration
 f =. m2*a :: f=. (f>W2)*f-W2    NB.subtract gravity
force after lift-off
 wp2=.(-/f*T[2 1]) % I2  NB. vector cross product of force and
normal to vector T gives unbalanced
torque. Divide by sling inertia for
acceleration about point P.
 w2 =. w2 +dT*wp2        NB. angular velocity of sling
 A2 =. A2 -dT*w2         NB. new sling angle
 S =. P + T=. SL*(1 _1)*1 2 o A2 NB. new T & projectile coord.
                     NB. store new values in result matrix TB
 TB =. TB, (N*dT),P,S,A,w,wp,A2,w2,wp2,V,a,T
 goto Ln0 if. (160 > N=. N+1) NB. increment counter and repeat
calculations if limit is not
exceeded.
 TB[;6 9] =. TB[;6 9]* 180% o.1 NB. I like angles in degrees

I welcome coments and suggestions toward improving ZOT

Tue, 03 Mar 1998 03:00:00 GMT

 Page 1 of 1 [ 1 post ]

Relevant Pages

Powered by phpBB® Forum Software