Launching a Trebuchet

The Superior Siege Weapon

The trebuchet is a type of siege engine consisting of a lever-like throwing arm which launches projectiles via a sling. The throwing arm is rotated by the force of gravity acting on a massive counterweight fixed to the opposite end.

By I, Luc Viatour, CC BY-SA 3.0, Link.

The mechanics of the trebuchet are abstruse; requiring complex calculations of friction, torque, gravity, elasticity etc. They are out of scope for this article, but a thorough investigation is available in Trebuchet Mechanics (Siano, D., 2006) for the adventurous.

The following problem will assume Siano's "black box" model, where the mechanism for energy transfer to the projectile is unspecified. We are more interested in the equations of motion for the projectile once it has left the sling.

The System

Consider a trebuchet on a battlefield. Assume the terrain underneath and in front of the trebuchet is indefinitely flat. At time $t_0$, the projectile leaves the sling at height $h$ above the terrain, travelling at speed $v$ and angle $\theta$ above the terrain. Some time later, at $t_e$, the projectile hits the ground some distance away.

Choose a 2-dimensional cartesian plane as the coordinate system. The origin is fixed on the ground directly underneath the projectile release point. The $x_1$ axis is positive in the projectile's direction of travel, and the $x_2$ axis is positive above the ground.

From this coordinate system, the initial conditions become apparent; the projectile has position $x = (0, h)$, and velocity $\dot x = (v \cos \theta, v \sin \theta)$.

Equations of Motion

The naïve equations for the projectile's motion simply take gravity as the sole constant force,

$$\ddot x = - g \label{parabola}$$

which describes a parabola:

If it was that simple, the problem would be solved already. The analytic solution for equation \eqref{parabola} above is

$$x_2 = x_1 \tan \theta - \frac{g}{2}\frac{x_1^2}{(v_0 \cos \theta) ^2}\text{.} \label{parabola-solution}$$

Obviously reality is a little more complicated; there are other forces acting on the projectile as it flies through the air. An important force is drag.

Drag

The drag force due to air resistance dramatically changes the trajectory of the projectile. We will define it here as a function of the projectile's velocity,

$$F_D = -\dot x^2\frac{\rho A c_D}{2} \label{drag}$$

where:

• $\rho$ is the density of the fluid,
• $A$ is the reference area, and
• $c_D$ is the drag coefficient.

For simplicity, we will take $\rho$ and $A$ to be constant. The coefficient of drag, $c_D$, is a little more complicated, but assume it is always $> 0$. More on all of these values later.

The equation of motion, now taking drag into account, become:

$$\ddot x = -\dot {x}^2 \frac{0.5 \rho A c_D - g}{m} \label{drag-solution}$$

This new force acting on the projectile dramatically changes the trajectory of the projectile. The figure below shows the trajectory under drag, plotted with a Leapfrog integrator, compared to a parabola.

What's left is a system of ODEs which are more complicated than a simple parabolic trajectory. We will solve them using numerical methods.

The Problem

A trebuchet is capable of launching heavy projectiles over great distances. Supposing you had a particular target in mind, you'd need to be able to adjust the launch parameters in order to hit it.

To simplify things, assume that the projectile always launches at $\theta = 45{\unicode{xb0}}$ and $h = 15 \si{\metre}$. Additionally, fix the projectile as a smooth sphere with mass $90 \si{kg}$.

Suppose you have a target at distance $r$ in front of the trebuchet. Now determine what the initial velocity must be in order to have the projectile land at $(r, 0)$.

The problem becomes a boundary value problem, with boundary values $x(t_0) = (0, 15)$ and $x(t_e) = (r,0)$. With all the other information available, what's left is to determine what the initial velocity must be.

Physical Properties

Gravity

Assuming a flat earth approximation, acceleration due to gravity is a constant vector

$$g = (0, 9.81 \si{m/s^2}) \text{,} \label{gravity}$$

meaning gravity only accelerates the projectile in the $x_2$ dimension. It is set positive here because it is negated in \eqref{drag-solution}.

Projectile Shape

The projectile is described as a $90 \si{kg}$ sphere, but the equations of motion do not include the mass. They do however include the reference area; the cross section through the centre of the sphere which is essentially the area of a circle of equal radius. So the radius of the sphere is needed.

Medieval siege engineers might have used granite stones as trebuchet ammunition (most sources just say "stone".) A quick search turns up a density measurement for granite: $\num{2.65e3} \si{kg/\metre^3}$. So a $90 \si{kg}$ granite ball has a volume of $\frac{90}{\num{2.65e3}} = \num{3.396e-2} \si{\metre^3}$.

Plug the volume into the sphere volume equation and solve for the radius:

\eqalign{ \num{3.396e-2} &= \frac{4 \pi}{3} r^3\cr r^3 &= \frac{3 (\num{3.396e-2})}{4 \pi}\cr r &= (\frac{0.102}{4 \pi})^{1/3} \cr r &= 0.2009 \si{\metre} } \label{radius}

A radius of $r = 0.2009 \si{m}$ gives a reference area for the drag force equation $A = \pi r^2 = 0.1268 \si{m^2}$.

Air Density

Air density in reality is not a constant value. It depends on the temperature, humidity and altitude. We will assume for simplicity that the air density during the trebuchet launch remains constant for the entire flight of the projectile.

A search gives the standard sea level air density as $1.225 \si{kg/m^3}$.

Coefficient of Drag

The drag coefficient depends on the Reynolds number of the air relative to the velocity of the projectile. The Reynolds number for the sphere is given by

$$Re = \frac{2r \rho |\dot x|}{\mu} \label{reynolds}$$

where $\rho$ is the same air density as above, and $\mu$ is the dynamic viscosity of air. Set $\mu = 1.789 \times 10^{-5} \si{kg/m.s}$, the dynamic viscosity of air at standard sea level.

Now this is where it gets a little complicated. The drag coefficient is not a well defined function of $Re$; it varies in a way that makes it hard to describe. You could use this reference (Cimbala, J., 2012) to build a big table of $c_D$ vs. $Re$ and perform linear interpolations to get the correct $c_D$ for the given $Re$ value. Instead take it as a given that for this simulation $10^5 \le Re \le 10^6$, and simplify the calculation for the drag coefficient as follows:

$$c_D = \begin{cases} 0.47, & \text{if}\ Re < \num{2e5} \cr 0.1, & \text{otherwise.} \end{cases} \label{drag-coefficient}$$

The sharp drop at $Re = \num{2e5}$ simulates the transition from laminar to turbulent flow.

Solving the BVP

Finally, we get to solving the boundary value problem. The equations given above are good candidates for the shooting method; you can expect $v$ to be an increasing function of $r$.

For integration I chose to use the Leapfrog method. The integrator continue to step forward in time until $x_2 \le 0$, at which point it will return the projectile's range $x_1$.

I used the Secant method method to perform root finding. Two initial values for the launch speed, $v_1$ and $v_2$, are chosen and the ranges $r_1$ and $r_2$ are computed by performing the integration described above for each initial speed. The difference between the true range $r_t$ and the computed ranges are then taken.

While these differences are sufficiently far apart, calculate a new initial launch speed

$$v_{new} = v_2 - (r_2 - r_t) \frac{v_2 - v_1}{r_2 - r_1}\text{.} \label{launch-speed}$$

Then set $v_1 = v_2, v_2 = v_{new}$ and repeat.

When $r_1 - r_2 < \epsilon$, then the root has been found. Take the midpoint $v_1 - \frac{v_1 - v_2}{2}$ as the answer.

Summary

This article had a little bit of everything; ODE integrators, root finding, and boundary value problems. Ultimately, it was a demonstration of the utility of numerical methods for real-life situations. Often, equations will become difficult or impossible to solve analytically; computers are needed to find the answer.

If a trebuchet engineer had just taken the easy route and used the parabolic model to estimate projectile trajectories, his attempts to hit the target would have fallen short. After accounting for drag and using the shooting method to solve the BVP, we get closer to the correct answer.

So, if we're releasing at $\theta = 45 \unicode{xb0}$ and $h = 15 \si{\metre}$, what is the initial speed required to launch a $90\si{kg}$ projectile exactly $300\si{m}$?

$53.528 \si{m/\second}$, or $192.7 \si{km/\hour}$ which is impressively fast.

There are a few things that will improve the efficiency and get the required launch speed down. The first is to change the launch angle; due to drag, $45 \unicode{xb0}$ is not necessarily the most efficient launch angle. Secondly, by raising the initial height $h$, a greater distance could be achieved for the same amount of work. To decrease drag, move the trebuchet to a higher elevation where the atmosphere is thinner. Also, switch the projectile material to something denser, like a lead ball; a smaller reference area would result in less drag.

The animation below shows the shooting method in action; converging upon the correct launch speed for the given launch angle and boundary values.

Code

Below is the Python code used to solve the BVP described in this article. To use it, call shoot with the start and end points of the required trajectory, the launch angle, the time step and the required accuracy. It will return the speed required to achieve the required trajectory, for example:

#!/bin/python2
import math

theta = 45 * (math.pi / 180) # 45 degrees in radians
dt = 1.0e-2
e = 1.0e-3

xa = [0.0, 15.0]
xb = [300.0, 0.0]

speed = shoot(xa, xb, theta, dt, e)


From here, you can also call the integrate function and pass in the optional traj variable. This variable will be populated with the trajectory of both $x$ and $\dot x$ as calculated by the integrator.

traj=[]
integrate(xa, xb[1], speed, theta, dt, traj)

for t in traj:
print ' '.join(map lambda x: str(x), t)


Shooting Method Code

#!/bin/python2
import math
G = [0.0, -9.81] # m/s^2

PROJ_MASS = 90 # kg
PROJ_DENSITY = 2.65e3 # kg/m^3
PROJ_VOLUME = PROJ_MASS / PROJ_DENSITY # m^3
PROJ_RADIUS = (PROJ_VOLUME / (4.0/3.0 * math.pi)) ** (1.0 / 3.0) # m
PROJ_CS = math.pi * PROJ_RADIUS ** 2 # m^2

AIR_DENSITY = 1.225 # kg/m^3
AIR_VISCOSITY = 1.789e-5 # kg/(s.m)

DRAG_FACTOR = 0.5 * AIR_DENSITY * PROJ_CS # kg/m

# Limit integration time
M_MAX = 1000
N_MAX = 100000

def mag(x):
""" Compute the magnitude of a vector. """
return math.sqrt(reduce(lambda i,j: i+j**2, x))

def reynolds_number(dx):
""" Get the Reynolds number as a function of current velocity. """
v = mag(dx)

def drag_coefficient(dx):
""" Get the drag coefficient as a function of current velocity. """
Re = reynolds_number(dx)

if Re < 2e5:
return 0.47
else:
return 0.1

def acc(dx):
""" Calculate the acceleration as a function of the current velocity."""
cd = drag_coefficient(dx)
df = (DRAG_FACTOR * cd) / PROJ_MASS
return map(lambda g,v: -v**2 * df + g, G, dx)

def step(x,dx,dt):
""" Advance x by its rate of change dx """
for i in xrange(0,len(x)):
x[i] += dx[i] * dt

def integrate(x, h, speed, theta, dt, traj=None):
""" Perform leapfrog integration. Stop when the
projectile hits the ground. """
n = 0
x = list(x)
dx = [math.cos(theta) * speed, math.sin(theta) * speed]
d2x = acc(dx)

# Record starting conditions
if traj is not None:
traj.append(x + dx)

# Half-step velocity
step(dx, d2x, dt/2)

while x[1] > h:
step(x,dx,dt)

# Get acceleration
d2x = acc(dx)

if traj is not None:
# Temporary half step velocity
# to align vt with xt
tempdx = list(dx)
step(tempdx, d2x, dt/2)
traj.append(x + tempdx)

# Do full step
step(dx,d2x,dt)
n += 1
if n > N_MAX:
raise Exception('Too many steps (x={})'.format(x))

# return the distance travelled by the projectile
return x[0]

def shoot(xa, xb, theta, dt, e):
""" Use the shooting method to find the correct launch speed. """
# Pick a starting speed
speed1 = 50
speed2 = 55
h = xb[1]

# Shoot for xb
r1 = integrate(xa, h, speed1, theta, dt)
r2 = integrate(xa, h, speed2, theta, dt)

m = 0
while abs(r1 - r2) > e:
f2 = r2 - xb[0]
speed3 = speed2 - f2 * (speed2 - speed1)/(r2 - r1)
speed1 = speed2
speed2 = speed3

r1 = r2
r2 = integrate(xa, h, speed2, theta, dt)

m += 1
if m > M_MAX:
break
return speed1 - (speed1 - speed2)/2