Calculating impact forces, or the theory of bouncing bodies

MATHML

This page contains some formulas in MATHML coding. In fact this is my first attempt at MATHML. MATHML is too elaborate to write by hand. So the actual formulas are generated from a kind of ‘shorthand’ code (in the source code of this page) by means of a piece of Javascript called ASCIIMathML.js. You must have Javascript enabled to see this page correctly.

The math is displayed correctly if you use a modern browser (Netscape 7.1, Mozilla, Firefox), but you will need some math fonts:

  • In Linux, you need the ttf version of some LaTeX fonts. In Debian, you can get them by installing the package latex-xft-fonts. Your browser may complain that it wants more math fonts (like the Mathematica fonts), but for this page, they are not necessary.
  • In Windows when using Mozilla et al., if you haven’t viewed MATHML pages before, you will be shown a link to a font installer.

If you use MS Internet Explorer, you may have to install MathPlayer first; this will provide IE with MATHML capability.In Debian Linux (the newest version, Sid), you must set the environment variable MOZ_DISABLE_PANGO to 1; otherwise the Debian browsers (‘iceweasel’ or ‘iceape’, as Debian now prefers to call them) will make a complete mess of MATHML display. Setting this environment variable may not be necessary (and may very well be harmful) in other distributions of Linux, in particular FC6.

Printing MATHML works fine on Windows. On Linux systems you need the newest Firefox/Iceweasel (3.0) for MATHML printing to work.

1. Introduction

In high-school physics it is taught that in a collision, the change in momentum that a body undergoes is equal to the collision force times the duration of the collision:

Δ(mv) = FΔt

or more accurately:

Δ(mv) = ∫F(t)dt

You can do some nice sums with this (‘assuming that the duration of the collision is 1 millisecond, what is the force?’) but in general, elementary textbooks give no hint about calculating F and Δt separately.

When I was an undergraduate physics student many years ago, I did some experiments with bouncing balls. I tried to analyse the problem theoretically, but did not succeed. I could, however, analyse the bouncing of shapes of uniform cross section, like a rectangular block.

Assume a rectangular body falls to the ground vertically. What happens? See the picture below, which shows the various stages of the ‘bounce’. We assume that the block is perfectly elastic (deformable without energy loss), and the ground perfectly un-yielding (infinitely ‘stiff’, and infinitely massive), so as to resist any impact. Actually by an ‘un-yielding’ surface I mean: a surface which will not allow any mechanical energy to be transferred upon it, having infinite ‘mechanical impedance’. Think of a rectangular rubber eraser falling down upon a stone slab. We disregard air resistance.

bounce-t

At moments A and B, the block is still in free fall. It does not ‘feel’ any force acting upon it; it is in an unstressed state (indicated by green). C is the moment when it hits the ground.

At the moment C, the top of the block ‘does not know’ that the bottom has already hit the ground. The top goes right on falling! Not for long, however. The bottom of the block does feel the impact. It goes into a stressed (and compressed) state, indicated by red in the picture. A shock wave moves up through the block at the speed of sound (i.e. the speed of sound in the solid, which can be much higher than the speed in air – and normally very much higher than the ‘fall speed’ of the block itself).

At the moment E, the compression wave has reached the top of the block. Now every particle in the block is momentarily standing still; all kinetic energy has been transformed into potential (i.e. elastic) energy. At this moment, the bottom of the block does not ‘know’ that the compression wave has reached the top; this information could not have reached it. But the top knows; it starts to ‘un-compress’ – and thus moving upwards. A ‘decompression wave’ starts to move downwards through the block – again at the speed of sound.

At the moment G, the decompression wave has hit the ground. The block is now again unstressed, and moving upwards (as a whole). At H, the block is free of the ground, and again in ‘free fall’ (as the physicists say; of course it is not falling in the common-sense meaning of the word, but moving upwards).

By putting some numbers into this, and using the formula for the speed of (longitudinal) sound in solids, it is easy to deduce, in the case of a perfectly elastic block:

  1. How long the ‘bounce’ takes.
  2. With which force the block hits the ground (e.g. as a function of the initial height from which the block falls).
  3. How much greater this ‘impact force’ is than the ‘static force’– the force that the ground feels when the block is just resting upon it.
  4. Which stress is experienced by the block (and also which strain, in the case of a perfectly linearly-elastic block), and of course also how much greater this is than the stresses in the block when it it is resting on the ground statically.
  5. That this really makes sense: all the numbers ‘add up’.

I give the derivations below, but any reader could easily do this by him/herself. All the essential physics has already been mentioned above.

2. Duration of the impact

The most remarkable thing about the impact of a rectangular block (or a cylinder) is that the region behind the shock front (red in the picture) is uniformly stressed. I remember how amazed I was when I discovered this. I found it by some complicated reasoning involving Fourier series. For the time being the reader should just take this on faith. This also means that the stress at the bottom of the block during the collision, and therefore the force of the impact (stress × cross-section of the block) is constant during the impact.

The speed of sound in a solid, σ, is σ=Eρ where E is the elasticity (Young’s modulus) of the solid and ρ is its density. The duration of the collision is:

Δt=2Lσ

where L is the length (or height) of the block, because the sound wave has to run through the block twice. A typical value for the speed of sound in a solid, e.g. in steel, or in wood along the grain, is 5000 meters/second; in soft materials it is lower, e.g. in natural rubber it is about 1600 m/s.

3. Force of collision

If the mass of the block is m, its falling speed is v, and the collision is 100 % elastic, so that after the bounce the upward speed is also v, the change in momentum is 2mv, and the impact force is

Fimpact=2mvΔt=mvσL

If the block started falling from rest, from a height h (the height of the bottom of the block above the ground), the speed of impact is given by conservation of energy:

12mv2=mgh

where g is the acceleration of gravity (about 10 m/sec2). So:

vimpact=2gh

So the impact force can be written as

Fimpact=mσL2gh

This is interesting by itself; if the falling body is some kind of rod of a given density and cross-section, and we double its length, its mass will also be doubled – so the impact force remains the same! It just lasts twice as long.

We can also calculate the time that it takes for the block to fall from height h; it follows from the definition of g:

tfall=vimpactg=2hg

When the block is just resting on the ground, the force felt by the ground is

Fstatic=mg

So,

FimpactFstatic=σL2hg=tfalltsound

where tsound is the time it takes for a sound wave to run through the block once (i.e. half the impact time). This ratio of ‘dynamic force’ to ‘static force’ could also be called the extra ‘g force’ that the falling object feels due to the impact (because, when just resting on the ground, it feels a force of ‘one  g’).

4. Elastic energy

At the moment E, the whole block is stressed. Because we assume perfect elasticity, the block will be strained (shortened) by an amount ΔL given by

EΔLL= stress =FimpactA

where A is the cross-section of the block. So:

ΔL=LFimpactAE

So we can regard the block as a spring with ‘spring constant’ S (=force/displacement) given by

S=AEL

The elastic (potential) energy of a such a spring, when the displacement is ΔL, is, according to another well-known high school physics formula:

P=12S(ΔL)2

By substituting our earlier results into this formula, we get the satisfying result:

P=Fimpact2L2AE=m2σ2ghLAE=m2ghLAρ=mgh

because LAρ=m, the mass of the block. So really all the original potential energy of position has been transformed into elastic energy. This should give some confidence in the consistency of this theory.

5. Preliminary discussion

The key result is the equation in section 3:

FimpactFstatic=tfalltsound

This by itself gives already a few practical results:

  • Impact forces (caused, e.g., by a fall) can be much greater than static forces, because the speed of sound in solids is so high (and therefore tsound so small). This is a fact of experience (balancing a brick on your head is fairly easy and painless, in contrast to having it fall on your head from a height of several meters), but now we can explain this fact.
  • The longer you fall, the harder you hit the ground. This, also, is obvious, but isn’t it nice to have this proved scientifically?
  • If you want to survive a fall, try to make tsound as long as possible. It is better to fall with your feet down than horizontally. Don’t curl yourself up into a ball. Even better, wrap yourself in thick layers of soft stuff with low speed of sound (like NASA did with the Mars Pathfinder). Note: this interesting page gives more elaborate tips about surviving a fall, including the advice to fall feet first, but wisely cautions that it is better to avoid falling in the first place.

More interesting is the case of a section of a skyscraper, say of 50 floors each of 4 meters in height, which suddenly finds itself in free fall because some miscreants, by flying a huge aeroplane into it, blew away a floor below, or weakened it beyond its carrying capacity by softening the steel columns keeping the structure up. The following analysis is somewhat simplified (we should take into account both the elasticity of the falling block and of the bottom half of the skyscraper) but I believe it is essentially correct.

We would have (for a fall of 4 meters):

tfall=2×410= 0.9 s

while

tsound=50×45000= 0.04 s

The ratio of the two values is about 23 (¹). After falling through the height of one floor the section would experience a dynamic stress equal to about 23 times the static stress (that the building was designed for). As no engineer would design a structure with a safety factor of that magnitude (safety factors of about 3 are common in buildings), further failure of the structure would certainly occur.

¹ The actual number depends a little bit on assumptions, like the number of floors in the falling section of the building. The 2002 paper by Bazant and Zhou estimates a ratio of 31 (formula 1 in the paper is almost identical to the one in the present article; it just adds the static force to the total force).

The ‘overstress’ will communicate itself almost instantly (at the speed of sound in the steel supports) through the whole structure. As a result, the building may very well fail in many places at once. Its fragments would then tumble to the ground at free-fall speed. An experimental confirmation of near-simultaneous breakage at many spots in a structure can be seen here.

Many ‘9/11 conspiracy theorists’ say that the collapse of the WTC towers ‘would have violated elementary physics’ unless demolition charges had been placed throughout the structure. There are thousands of websites that claim this (²), but I suspect their authors do not really know much about elementary physics.

²Just look up ‘9/11 physics wtc’ (without the quotes) on Google. Almost all results point to sites believing in ‘controlled demolition’. An honourable exception is this site.

(8 August 2006, revised 31 August 2006, 3 November 2006)

6. Further development of the theory: the ‘mass-spring chain model’

The theory so far is very rudimentary. Several things are wrong with it, like:

  • The reader is asked to accept something on faith.
  • The ground is supposed to be infinitely unyielding.
  • The theory is one-dimensional; for instance it does not take into account that objects which are compressed lengthwise also expand sideways (Poisson’s ratio).
  • The theory only handles objects of uniform cross-section.

WARNING. From here on this article will become slightly more difficult. However, if you read on to the end, you will be REWARDED with an applet that shows the actual bouncing in progress!

6.1 Modelling the problem

Addressing the first point first, we can gain some more insight by making a ‘model’ of the falling block. We model the block as a chain of massless springs and point masses. A model in which we have one spring topped by one mass is obviously very inaccurate (at least it conflicts grossly with the theory sketched above), because the collision will just be one half-cycle of a harmonic oscillator, with the compression of the spring (and therefore the force on the ground) being a half-sinusoid as a function of time. Above, I asserted that the force would be constant during the collision.

hbounce2-t

In the picture above, the block is modeled as a chain of four masses and four springs. In general we would have N masses and N springs, and we may expect the model to become more accurate as N is increased. To save space, I turned the picture ninety degrees, so the collision is horizontal. The chain is shown moving with velocity -v, about to impinge upon the perfectly unyielding block B.

(22 September 2006)

Our chain must be a model of a solid block with unstressed total length L, mass M, and spring constant (stiffness) S. Let us use lower case letters for the corresponding quantities in the individual springs and masses. Then

m=MN

 

l=LN

 

s=SN

The last expression may not be immediately obvious, but a little thought will show that it is true. This shows, incidentally, that when you go to smaller and smaller dimensions, matter becomes harder and harder!

The instantaneous positions (x-coordinates) of the individual masses are x1,x2,…,xN, the lowest-numbered mass being the one closest to B.

The only forces acting on the masses (and therefore their accelerations) come from the springs (we neglect gravity here; that was another reason for drawing the picture horizontally). Mass N only experiences force from one spring (from the left); the other masses are acted upon by two springs.

Let us call the extension of spring k (the amount that it is longer than its unstressed length) ek. We have:

ek=xk-xk-1-l   (for k>1)

For the leftmost spring, e1=0 when x1>l. Otherwise, e1=x1-l (a negative value).

A spring with stiffness s and extension e pulls with a force se (if e is positive; otherwise it pushes, of course). So the net force on mass k, and its acceleration x..k, are related (by good old Newton’s second law, and using Newton’s notation: a double dot above means double differentiation with respect to time) as follows:

mx..k=-s(ek-ek+1)

(please check the signs; this is pretty crucial). Of course, for k=N,  ek+1 is zero, because there is no spring to the right.

Now let us shift the time coordinate so we have t=0 at the moment the chain just touches B (i.e. the instant C in the first picture). From then on (at least for a while) we will have e1=x1-l.

Let’s also shift the space coordinates as follows:

yk=xk-kl

At time t=0 therefore, yk=0 for all k. Also, because we shifted each xk by a constant, x.k=y.k and x..k=y..k. We also have at all time (check this!):

ek=yk-yk-1

which is a little bit simpler than the previous expression for ek.

Now let us treat the yk, all together, as a vector. I.e. we use y→ as shorthand for {y1,y2,…,yN}. Then, by substituting one thing into another (the reader may need some scrap of paper to verify this), we get

d2dt2y→=-smAy→=-N2(SM)Ay→

where A is a tridiagonal matrix:

A=(2-1-12⋱⋱⋱-1-12-1-11)

On the diagonal of this matrix, all elements Akk are equal to 2, apart from the last one, which is 1 (corresponding to ‘only one spring’)! The elements just above and below the diagonal are equal to -1. All other elements are zero.

Now it just a simple (!) matter of solving this equation system. It is a set of N second-order differential equations, so we also need 2N initial conditions. They are: at t=0: yk=0  and y.k=-v   for all k.

The key is to find a way to completely diagonalise the matrix A. Amazingly, this is possible (at least in principle) without resorting to numerical methods. Some searching on the Internet revealed a recent paper by Wen-Chyuan Yueh, which proves that our matrix has eigenvalues λk given by (for 1≤k≤N) :

λk=2+2cos(2kπ2N+1)

You can check this formula by using a pocket calculator and a numerical ‘diagonaliser’ (this is a good one, on-line on the Internet). The formula shows that all eigenvalues are distinct, real, and positive. They become smaller with increasing k. Yueh also gives formulas for the eigenvectors.

Let’s work out an example for the first non-trivial case: N=2. This can be done even without Yueh’s formulas. The explanation will be lengthy and tedious, but it will serve as a guide to find the solution for the ‘general N’ case.

(4 October 2006)

6.2 The case N=2

In the case N=2 the problem can be solved by hand, although it takes some effort. The matrix A is

(2-1-11)

and the eigenvalues are the two roots of the equation det(A-λI)=0. We get:

λ12=3±52

This does not look like the cosine expressions of Yueh’s formula, but they are really the same. Think of the close association between the ‘golden ratio’ and the regular pentagon (for the case N=2, the expression 2N+1 in Yueh’s formula becomes 5). It is in fact useful to introduce here the ‘golden ratio constant’ φ=5-12, or 0.61803399, which can be used to simplify many expressions for the case under consideration (N=2). φ has many ‘magical’ properties (which can be checked easily), like φ-1=1+φ,  2+φ=1+φ,  φ-2=2+φ, etc.

The first eigenvalue can be written as λ1=2+φ, and the first eigenvector as:

u→1={1,-φ}f1(t)

Where f1 is some function.

The second eigenvalue is λ2=1-φ, and the second eigenvector is

u→2={φ,1}f2(t)

So

{1,-φ}f..1(t)=-smA{1,-φ}f1(t)=-N2(SM)(2+φ){1,-φ}f1(t)

Now for simplicity we assume SM=1 (by some scaling of the time coordinate), so we have (because N=2)

f..1(t)=-4(2+φ)f1(t)

This has a solution (which is, as required, zero for t=0)

f1(t)=sin(2t2+φ)=sin(2tφ)

Similarly

f..2(t)=-4(1-φ)f2(t)

with solution

f2(t)=sin(21-φt)=sin(2φt)

The solution that we seek, {y1,y2}, is some linear combination of the eigenvectors (or rather eigenfunctions) {1,-φ}f1(t) and {φ,1}f2(t). So:

y1=αsin(2tφ)+βφsin(2φt)

 

y2=-αφsin(2tφ)+βsin(2φt)

α and β have to be chosen so that y.1=y.2=-v at t=0. That means:

φ-1α+φ2β=-v2

 

-α+φβ=-v2

This yields (using the ‘magical’ properties of φ):

α=-(v2)φ3+φ=-0.0854102v

 

β=-(v2)2+φ2-φ=-0.9472136v

to fit the initial conditions. The resulting functions y1(t) and y2(t) can now be plotted with gnuplot (actually we are plotting the x1 and x2 now):

The view here is vertical again. In our timescale (with SM=1) the bounce takes approximately 2.63427 time units. The top mass starts from height 2. It descends more steeply than the bottom mass, which starts at height 1, and whose motion is a kind of ‘bathtub shape’.

We notice that the motion is not time-symmetrical. The ‘floor’ of the bathtub is not horizontal. Also, a numerical calculation shows that at the moment when the bottom mass reaches height 1 again (so that the bottom spring becomes free of the ground), its velocity is not exactly +v, nor are the speed and height of the top mass exactly +v and 2 respectively. Close, but not quite. In fact, at the moment of lift-off, the velocity of the center of gravity of the two masses is 0.9738149v. In other words, the coefficient of restitution of the ‘two masses and two springs’ system is 0.9738149. The ‘missing energy’ is in the vibration (driven by the top spring) of the two masses about their common center of gravity. An example of the ‘arrow of time’ in the macroscopic world of classical mechanics?

I wonder if the motion will become more symmetrical, and the coefficient of restitution closer to one, when we increase the number of masses and springs. Maybe they won’t! That would be pretty bad for the simple theory of sections 1-5.

(9 October 2006)

6.3 The general case for arbitrary N

Now we try to calculate the motion of the ‘ball and spring’ model for arbitrary N, by analogy to the above method for the N=2 case.

Let the eigenvectors u→k of matrix A be stored in a matrix U (so that column k of U is the k-th eigenvector; ukj is the j-th component of the k-th eigenvector).

Let the square roots of the eigenvalues be stored in a vector ω→ with components ωj. The eigenvalues are real and positive, so they have real square roots.

The ‘weights’ (quantities like α and β in the example above) are in vector w→ with components wj. w→ is determined by the initial conditions:

y.k(0)=∑jNωjwjukj=-v   for all k

This can solved by computing a vector d→ with components dj, which is the solution of

Ud→=v→  (with v→={-v,-v,…,-v})

and then putting

wj=djNωj    for all j.

Now the movement of mass  k can be described by

yk(t)=∑jwjukjsin(Nωjt)

I put this in a small C program, using GSL (the ‘GNU Scientific Library’) for the matrix operations. Now we can compute and plot the other cases. It turns out that we didn’t need to worry. Already from N=3 we see that the picture has become much more symmetrical, the bathtub of the bottom mass has become much flatter, and the middle mass also describes a bathtub, but of half the width.

When we go to higher numbers of masses and springs, all the sinusoids combine to produce a clear triangular pattern:

On the sides of the triangle, we see the ‘compression wave’ moving upward until it meets the descending top of the chain, then moving down just as quickly. This is exactly like the behaviour of the solid block, as asserted in paragraphs 1-5.

These numerical experiments also produce some other things: as N increases, the restitution coefficient becomes closer to 1 (although not completely uniformly), and the duration of the collision approaches 2 (in our time units). The case N=1 is the trivial case mentioned at the beginning of paragraph 6.1. It takes π time units. It cannot ‘lose’ energy to internal degrees of freedom (because it doesn’t have any), so its restitution coefficient is 1.

N duration Coeff. of restitution
1 3.141592654 1.0000000
2 2.634270792 0.9738149
3 2.429536835 0.9706905
4 2.332784360 0.9741725
5 2.276606027 0.9736748
6 2.232998859 0.9742083
7 2.204197093 0.9756464
8 2.182579247 0.9759477
16 2.100138195 0.9805711
32 2.055544323 0.9851961
64 2.031227683 0.9892320
100 2.021683613 0.9913945
1000 2.003566250 0.9976620

The time approaches 2? In what units? Well, of course in units of the time the compression wave takes to run through the chain (or the block of which the chain is a model). See the next section (still has to be written, unfortunately).

Comments
Simplue WordPress theme, Copyright © 2013 DicasLivres.org Simplue WordPress theme is licensed under the GPL.