Physics Lab Report

Ph2A Win 2020 Numerical Analysis Lab

Max Yuen

Mar 2020

(use g = 9.8m/s2 for all problems.)

Background Many physics problems cannot be solved directly by hand or analytically. We resort to numerical methods to give us approximations to the problem. In this lab you will learn the Euler method, which allows you to solve Newton’s laws of motion. This is done by treating the velocity as a piecewise linear function with many time intervals and during interval the acceleration is assumed to be uniform. This allows us to use the kinematic equations we learned about in the first half of the class to approximate the motion. If we choose to partition the motion into smaller time intervals, the approximation becomes much better since the differences between adjacent intervals become smaller. In this lab, this numerical analysis method will be applied to the motion of a falling object under the influence of gravity and drag force. If you are adventurous, you can even try to extend this to 2D and compute the realistic trajectory of a baseball. You might even try some other problems, such as a mass attached to a spring.

Euler’s Method Foundations This method is well suited for problems where the acceleration is a function of the velocity, as in the case of a falling object under the influence of gravity and drag:

a = f(v) (1)

Falling object with drag force The model for drag fits the prescription for using Euler’s method since the net force on a falling object with drag is given by:

ma = −mg −FD (2)

ma = −mg − 1

2 ρairACDv

2 · sgn(v) (3)

a = −g (

1 + ρairACDv

2 · sgn(v) 2mg

) (4)

a = f(v) ← Equation of Motion (5)

where m is the mass of the falling object, a is the acceleration of the object (which is positive when pointed up), ρair is the density of air (about 1.29 ·10−3kg/m3), A is the cross-sectional area, CD is the drag coefficient, v is the object’s velocity, and sgn(v) is the signum function which returns the sign of the argument. The second signum function is there to guarantee that the direction of the drag force is always in the opposite direction of the velocity function. Note that we see that the acceleration is an explicit function of v, which sort of makes this a chicken or egg problem. This is because we need a to get v, but to get a we need v, so which one do we compute first? Hold that thought. We’ll talk more on how to program this in EXCEL or Google Sheets later.

1

Figure 1: FBD for an object falling under the pull of gravity and resistance by drag force

Terminal Velocity In lecture, we talked about how after waiting for some time, if the object started at rest the speed will increase and the drag force will also become larger and eventually balance out with the gravitational force. When this happens, we have reached terminal velocity vterm = −v. This can be solved by setting a = 0:

0 = −mg − 1

2 ρairACDv

2 · sgn(v) (6)

2mg = ρairACDv 2 term (7)

→ vterm = √

2mg

ρairACD (8)

Using this definition for the terminal velocity, you can actually simplify Eq. 4 after a little bit of algebra to give a simpler version of the equation of motion.

a = −g (

1 + v2

v2term sgn(v)

) (9)

The Algorithm What we didn’t talk about or was able to do was to compute the exact evolution of the velocity of the system from t = 0 to see how velocity went form 0 to −vterm. We can achieve this numerically without knowing how to solve a differential equation if we divide the problem into many smaller time intervals ∆t, say starting from ti = 0 to tf = N∆t, where N is the number of time intervals from the start to the end time. When you increase the number of intervals N or reduce the interval ∆t the approximation becomes better. (It should, or you will have an unstable solution.) Thus, we have N + 1 values of time tn = n∆t for N intervals from t0 to tN.

The Euler’s method will solve the equation of motion problem in each time interval by com- puting new velocities and new positions using the information from the previous interval. Let’s define our initial values for a0 and v0 at t0 = 0:

a(t = 0) = a0 = −g (10) v(t = 0) = v0 = 0 at rest initially (11)

We compute the new values at the next time step t1 from the old values by direction applica- tion of kinematics equations. The new velocity v1 is computed from the old velocity v0 and old acceleration a0. The new acceleration a1 is computed from the new velocity v1.

v(t = t1) = v1 = v0 + a0∆t (12)

a(t = t1) = a1 = −g (

1 + v21

v2term sgn(v1)

) (13)

2

Since the acceleration is negative, the new velocity will be negative and the term inside the parentheses will be smaller than 1. This means the drag force will begin to counter the effects of gravity. To compute the next set of acceleration and velocity, we repeat the same process, but the a1 and v1 values become old values and a2 and v2 become new values.

v(t = t2) = v2 = v1 + a1∆t (14)

a(t = t2) = a2 = −g (

1 + v22

v2term sgn(v1)

) (15)

And we continue to churn out new values by repeating this process. The technical term is recursion and the recursive relation is given by:

v(t = tn+1) = vn+1 = vn + an∆t (16)

a(t = tn+1) = an+1 = −g (

1 + v2n v2term

sgn(v1) )

(17)

Thus, you can now compute to arbitrary precision the velocity as a function of time for an object falling under the influence of gravity and drag. Solving for the position is as easy as solving another set of kinetmatics equations once the velocities are know. This is best done by using the average velocity in each time interval.

Procedures You’re welcome to look at the example Excel worksheet to get some ideas as to how to do this, but here is a quick tutorial on getting things started:

1. Open Excel

2. Enter the initial values on the spreadsheet. This will allow you to quickly change things later. Your table should look something like this in excel: (note that mass is set to 1)

parameter value units g 9.8 m/s2

vterm 0.7 m/s v0 0 m/s dt 0.01 s t0 0 s

Table 1: Initial values for modeling drag

3. Next, find an empty section in your spreadsheet and enter the column headings and the first row of initial values:

time(s) velocity(m/s) accel(m/s2) 0 0 -9.8

Table 2: Excel Worksheet Example

4. Next, enter for appropriate formulas in the next set of cells, using the appropriate algorithm. See EXCEL file for hints on the equations to enter.

5. In the box for time, we simply add the value stored at t0 to the previous, so 0 + 0.1 = 0.1.

6. In the box for velocity, we update this using the old velocity and the old acceleration using the formula v1 = v0 + a0∆t . Thus, we get 0 + (−9.8) ∗ 0.01 = −0.098.

3

time(s) velocity(m/s) accel(m/s2) 0 0 -9.8

0.01 -0.098 -9.60792

Table 3: Excel Worksheet Example

7. In the box for accel, we calculate the acceleration using Eq. 17. Thus, we get −9.8 · (1 + (0.0982/0.72)(−1)) = −9.60792...

8. This is for one iteration.

9. Repeat Steps 4 to 7 for as many as you need to get to terminal velocity.

10. Make a graph showing the velocity as a function of time. Can you see when the velocity becomes 90% of the terminal velocity?

11. Make another graph showing the acceleration as a function of time. Can you see when the net force become zero?

12. Try changing the ∆t to 0.1 or 0.001. What do you find? Explain what you found.

13. Try changing the vterm to 2.1 m/s. How does this affect your graphs?

14. Try a different starting velocity, say vi = −2.

15. Now, remember your results to the coffee filter drag force lab. Does it make sense?

4