ME211 Class Project-2 (Due midnight Friday, 12/15/06)

The goal of this project is to simulate the motion of an ice hockey puck in an elliptical rink. The project has two parts. The basic part is for a puck moving in one dimension, i.e., back and forth along a straight line between two ends. The advanced part is to generalize the motion to two dimensions, in the elliptical rink.

For simplicity, we neglect the size of the puck and treat it as a dot. The factors we try to simulate are surface friction and inelastic collisions with the boundary.

Part 1. 1D motion

Suppose the two ends of the linear rink are at x=0 and x=L. The motion of the puck under the friction force can be described by the following first-order ODEs,

.

The coefficient f represents the friction. The coefficient c is to simulate the fact that friction decreases as velocity (and also temperature) increases. The initial condition is v=v0 and x=x0, 0<x0<L. To simulate the collision, let v=-v* whenever x reaches 0 or L. Here v* (v) is the velocity right before (after) the collision. The coefficient of restitution <1 describes the energy loss during the collision.

Write a Matlab program to plot v(t) and x(t). Also find the final position of the ball, xf, when the ball stops (v=0). It is desirable to group different aspects of the task, such as ODE solver and collision, into different modules (or functions.)

Use the following parameter set as your test case, L=100 ft, x0=50 ft, v0=100 mi/h, f=1 m/s2, c=1.810-2s-1 and =0.55. You should choose other necessary simulation parameters in your program, such as the time step, to obtain a reasonably accurate final answer.

Part 2. 2D motion

Generalize your program to 2D. The rink boundary is an ellipse described by the equation x2/a2+y2/b2=1. The motion is described by

.

Here, v is the magnitude of the velocity, v=(vx2+vy2)1/2. The initial condition is (vx, vy)=(vx0, vy0) and (x,y)=(x0, y0). For the collisions, keep the velocity component parallel to the boundary unchanged and change the perpendicular component to v=- v*. Plot the trajectory of the ball in the (x,y)-plane. Use the following parameter set as your test case, a=100 ft, b=42.5 ft, x0=0, y0=0, vx0=50 mi/h, vy0=86 mi/h, f=1 m/s2, c=1.810-2s-1 and =0.55.

Grading

The project (the M-files and a short report) should be emailed to me211UR@gmail before the deadline. It will be graded on the following factors: (1) on-time turn-in (1 pt); (2) functioning programs and accurate results for the testing cases (15 pt); and (3) the written report, which should concisely present the algorithms and program structure, the results, the test of the programs, and the difficulties you have encountered (4 pt). Comments and suggestions on the whole class experience are also appropriate and welcome, but not mandatory.

Each group should finish the project independently. You are encouraged to talk to the instructor, but no one else.

Hints

  1. Time step size
  2. The step size depends on your choice of the ODE solver and the desired accuracy. You should choose your step size so that halving it will not change the results (within the desired accuracy.) This is the so-called convergence test.
  3. During collisions, the set step size will move the puck beyond the boundary (that’s how you know a collision will happen). When this happens, that time step needs to be reduced to move the puck right to the boundary and find its velocities there before performing further operations for the collision. One way of finding the right step size is to interpolate iteratively. For example in 1D, a certain step size t produces a new displacement x which moves the puck from x=xold to x= xold+x>L. A new step size tnew=t(L-xold)/x could be used to see whether it will move the puck to x=L. This process can be repeated until x is close enough to L.
  4. Similar care should be taken when deciding the final position of the puck due to friction. There, the set step size will reverse the puck velocity without collision (that’s how you know the puck will come to rest). When this happens, that time step needs to be reduced to find the exact time when it stops.
  5. Determining v|| and vin2D geometry during collisions

At each point (x,y) on the ellipse x2/a2+y2/b2=1, two vectors t, the unit vector in the tangential direction, and n, the unit vector in the tangential direction, can be defined. Their xy-components are

.

Then, v|| and vcan be related to vx and vythrough