8.13: Mathematics Lesson Extensions
Calculate the Center of Pressure
Overview: This exploration can be used as an enrichment exercise for a geometry class. Students will first need to know what a centroid is and how to find it for triangles and rectangles. They will need to be able to calculate areas of those shapes.
Learning Objective: Students will accurately calculate the centroid of a complex model.
Prerequisite: Geometry; can be extended to include calculus concepts.
Materials: Graph paper, ruler
Time Required: 80 minutes
Procedure: A profile of a rocket (like a shadow cutout) is used to calculate the center of pressure of the rocket. The centroid (the geometric center) of this profile can be computed, and it correlates with the center of pressure of the rocket. In flight dynamics, it is essential to know where the center of pressure is relative to the center of mass in order to assure stability. You should find through your calculation that the center of mass for the rocket will be "forward," meaning closer to the top of the rocket, than the center of pressure. Most of the mathematics required for this process is relatively easy, except for the nose cone if it is nonconical (which may require calculus).
The centroid of an object, like a rocket, is a calculation of the product of areas and distances from a reference point divided by the total area. Examine the rocket profile below and the mathematics that follow.
The entire shape can be divided into basic geometric shapes whose areas are easily calculated: the body is a rectangle and the fins are divided, as shown, into triangles and rectangles. Students should calculate the area of each region — all of which are easy to calculate, perhaps with the exception of the nose cone, which will be addressed later.
The centroid (shown in the figure above as a black dot with a cross) of each geometric region should then be located as accurately as possible. Centroids of rectangles are easy and intuitive. Methods for finding the centroid of a triangle can be researched in a geometry text or on the Internet.
For this twodimensional representation, we need to define two axes that are going to be our references for centroid calculations. Since the profile of the entire rocket is symmetric, it makes sense to have one of the orthogonal axes run along the axis of symmetry. Without any computation we can claim that the centroid of the entire rocket profile will be on this axis of symmetry. The other axis is chosen to be at some distance below the rocket passing through the reference point O (the origin, if you will). Students should draw a line from each centroid to the symmetry axis.
Since there is symmetry, the centroids of the fin sections will be directly across the axis of symmetry from each other, and the line that connects them will be perpendicular to that axis. The points of intersection (shown in the figure below as A, B, and C) are measurable distances from the origin, O. Points D (centroid of the rectangle) and E (centroid of the nose) are already on the axis.
Using the reference point O, measure the distances OA, OB, OC, OD and OE (recall that the position of point E will be discussed later).
The position of the centroid of the entire figure on the axis measured from O on Figure 2 is calculated by the following formula:
\begin{align*}d=\frac{\sum (\text{Area of object} * \text{distance})}{\text{Total Area}}\end{align*}
Letting the areas of the trianglerectangletriangle portion of the fins be represented by Area A, Area B, and Area C, respectively; the area of the body rectangle by Area D; and the nose by Area E; then the formula above would be:
\begin{align*}d=\frac{(\text{Area A} * \text{OA} * 2 + \text{Area B} * \text{OB} * 2+ \text{Area C} * \text{OC} * 2 + \text{Area D} * \text{OD} + \text{Area E} * \text{OE})}{(\text{Area A} * 2 + \text{Area B} * 2 + \text{Area C} * 2+ \text{Area D} + \text{Area E})}\end{align*}
The area and centroid of the nose must be calculated using the following guidelines. First, a curve match must be made. The nonconical nose can probably be matched by using a quadratic function. The vertex may have to be offset if the point on the nose is not blunt. Let’s look at a simple parabola first.
This area is calculated by a simple integration
\begin{align*}Area=2 \int\limits_{0}^{4} \sqrt{y}dy=\frac{32}{3}\end{align*}
The centroid is the point on the axis of the symmetry such that a horizontal line drawn through it would divide the area in half, so
\begin{align*}& 2 \int\limits_{0}^{h} \sqrt{y}dy=\int\limits_{0}^{4} \sqrt{y}dy\\ & \frac{4}{3}h^{3/2} = \frac{16}{3} \rightarrow h = 2.5\end{align*}
So, the centroid would be approximately 2.5 units from \begin{align*}y = 0\end{align*} and 1.5 units from \begin{align*}y = 4\end{align*}, which would be the measurement required for the example shown.
If a more pointed nose were required, then an offset quadratic function might be tried. Integrate what is needed (see the figure below). Although, if the edges of the curve are close to a straight line, as is shown in this example, it may be just as easy to approximate the centroid by using that of an isosceles triangle, which is just a third the length of the base altitude.
In fact, in this particular model, the value of h is 2.05, and, if it was calculated using an isosceles approximation, it would be exactly 2. The difference is within the error of measurement (assuming units are centimeters or inches). So it would seem that as nose sections become more linear, the isosceles approximation would be easier.
A student worksheet follows.
Alignment with National Math Standards: NCTM, Grades 9–12
Problem Solving, Communication, Connection and Representation:
 Problem Solving
 Communication
 Connection
 Representation
Handson Mathematics:
 Geometry
 Measurement.
Calculate Center of Pressure: Student Worksheet
Locate the centroid of the rocket profiled below. Each grid is 3 cm square.
Research Questions: The centroid of this profile coincides with the center of pressure of a rocket.
 Why is the location of the center of pressure of a rocket an important point to know?
 How does the center of pressure differ from the center of mass of a rocket?
 Both points, CoP and CoM, share a relationship to the rocket’s stability in flight. Can you determine what that is?
Calculate Center of Pressure: Student Worksheet Key
Area of a fin: 32.76 sq cm (approx)
Area of body: 180 sq cm
Area of nose: 27 sq cm
Sum of distances * Areas = 2*32.76*7.3 + 21*180 + 39*27 = 5311.3
Total area = 2*33 + 180 + 27 = 273
Quotient} = 5311.3/273 = 19.5 (approx)
So, the centroid of the profile is 19.5 cm from O (shown on the model as ).
Figuring out the centroid of the rocket profile as the center of pressure is a good approximation for locating the center of pressure of the model rocket. Accurate measurement of the center of pressure is a difficult engineering problem and takes a lot of resources. In the real rocket engineering world, the center of pressure is found by simulating air flow past the rocket using complex computational methods or by using wind tunnel testing and measuring aerodynamic forces. Both of these methods are expensive and time consuming, but locating the center of pressure on a real rocket is a very important parameter and needs to be located accurately. The center of mass of the rocket is the balance point and can be determined for model rockets by balancing them while suspended on a string. In order for a rocket to have stability in flight, the CoP and CoM should be a distance apart of at least 1.5 times the rocket body diameter. The CoM should always be higher (closer to the nose) for the rocket to be stable during flight.
Additional Resource: This NASA link explains this process for an airfoil, but the concept is the same as is done in this exploration — different variables, but the same concept.
http://www.grc.nasa.gov/WWW/K12/airplane/cp.html
Determine Maximum Height
Overview: This lesson requires knowledge of right triangle trigonometry and law of sines. In the original project, this math extension was used to determine the maximum height a rocket achieves while in flight, as shown in Figure 1, below. However, it is a nice standalone math lesson in using trigonometry along with two independent observers to measure heights of static objects like buildings or trees.
Learning Objective: Students will use triangle trigonometry to calculate the length of a side of a three dimensional object.
Prerequisite: Advanced algebra or precalculus.
Materials Needed: Measuring tape, inclinometer
Time Required: 80 minutes
Procedure: Within the basic lesson plan, students assume that their rocket flies directly vertical, in which case simple right triangle trigonometry can be used to determine the rocket’s maximum height. However, many rockets will fly at least slightly off a vertical path. When this happens, the triangle is no longer a right triangle, so the maximum altitude calculated with the tangent function would be inaccurate.
Distribute the student exploration, "Calculate the Height Your Rocket Achieves." Begin by outlining the one fundamental assumption: all measurements are going to be made from a level plane (ground level). Then, allow the students to brainstorm, with no other constraints, the problem of determining a method to accurately measure the maximum height of a model rocket. Students’ first inclination might be to use a simple right triangle with a single observer who will measure the angle of elevation as the rocket flies to its maximum (as shown above). If students agree that this is the best way, the teacher might ask if two independent observers might be better. It is essential that students understand that in all likelihood, the rocket’s path will not be exactly vertical; thus, the method used to calculate maximum height must be "independent" of the position from which the rocket is launched.
Pass out the student exploration "How high is that tree?"
Model 1:
Paul and Elise would like to estimate the height of a tree. Assuming they are on level ground, Paul measures an angle of elevation of 30 degrees to the top of the tree from point \begin{align*}A\end{align*}. Elise, standing 10 m from Paul and in a straight line to the tree, measures an angle of elevation to the top of the tree of 40 degrees. Using this information, find the height of the tree. Both Paul and Elise do not need to know how far they are from the base of the tree.
Model 1 Answer: approximately 18.5 meters
(Define height as H. Create tangent function equations for triangles ADC and BCD. Simultaneously solve these two equations.)
Model 2:
Now, suppose Paul and Elise do this. Paul stands at point A and Elise stands at point B. It is known that distance AB is 54 meters. Again, the assumption is that everything is on level ground, so point points A, B, C are in a plane that is perpendicular to line segment CD. Paul and Elise each measure the angle of elevation to the top of the tree. Paul’s measurement is 30° and Elise’s is 37°. Then Paul measures angle CAB as 15° and Elise measures angle CBA as 20°. They each use these measurements to calculate independent measurements on the height and compare. What do they find?
Note: To measure angles CAB and CBA, after recording the angle of max height, a vertical angle, students will move their inclinometer from the max height to a point on the ground directly underneath the max height (correlating to segment AC in the diagram). Students will then measure the angle parallel to the ground from that line to the line connecting the two observers (segment AB in the diagram).
Students should be guided through the mathematical derivation shown below, which results in a single equation appropriate for determining the height of an object. To aid students in following the derivation, pass out the worksheet "Deriving the Equation."
Measure values of length L, and angles a, b, c, d.
Author’s sidebar note:
It is a nice derivation to show that \begin{align*}\sin \angle ACB=\sin(b+c)\end{align*} by showing that the sine of an angle is the same as the sine of its supplement.
Model 2 Answer: Paul calculates the height as approximately 18.6 m; Elise calculates it as 18.4 m. Their average is 18.5 meters.
Ideally, two independent values of H (the maximum height of the rocket) are calculated and compared. If careful measurements are taken, they should be close. They may be averaged to obtain the final value.
Please note that this derivation does not reference the launch position at any time. Therefore, the maximum height can be accurately determined with this equation, regardless of whether any xdirection motion has occurred between liftoff and attaining max height.
Have students examine the similarities and difference between the two models: both models require multiple measurements to calculate a height that cannot be directly measured. Model 1 requires two observers, who stand a known distance apart and along the straight line segment AC, to each measure the angle of elevation to the top of the tree. Alternately, Method 2 has no requirements as to where the observers must stand, only that the distance between them be a known or measured value.
Have students identify which model would more accurately determine the maximum height if the object of interest were a rocket instead of a tree; have them explain their reasoning. If necessary, ask students if they think the path of the rocket will be exactly vertical and which model would fit better if it were not. Guide students towards understanding that Model 1, which requires the object to be measured and the two observers to be located along a straight line, would require the observers to know before the launch exactly where the rocket would be located when it is at max height. Since that is not possible, Model 2 would yield a more accurate result.
Note: There are several methods to solve analytically for the height in Model 2; only one of those methods is derived above. For more, access the following NASA site:
http://www.grc.nasa.gov/WWW/K12/rocket/rkthowhi.html
Student worksheets follow.
Calculate The Height Your Rocket Achieves
Student Exploration
In this exploration, you will determine a method to calculate the maximum height your rocket achieves. You must begin with the assumptions that:
 The plane on which your rocket is launched and from which you are measuring is an approximately level plane.
 You have a device that you can use to measure angles vertically (this would be referred to as an angle of elevation in a mathematics setting) and horizontally (the acute angle between two fixed points on the plane relative to your position).
Brainstorm various methods of determining this height measurement and describe and/or sketch them below. Select the one you think is best and explain your reasons.
How high is that tree?
Student Exploration
Model 1:
Paul would like to estimate the height of a tree. Assuming he is on level ground, he measures an angle of elevation of 30 degrees to the top of the tree from point A. He then moves 10 meters to point B and measures an angle of elevation to the top of the tree of 40 degrees. Using this information, find the height of the tree.
Model 2:
Now, suppose Paul does this. He enlists the help of his companion Elise. Paul stands at point A and Elise stands at point B. It is known that distance AB is 54 meters. Again, the assumption is that everything is on level ground so point points A, B, C are in a plane that is perpendicular to line segment CD. Paul and Elise each measure the angle of elevation to the top of the tree. Paul’s measurement is 30 degrees and Elise’s is 37 degrees. Then Paul measures angle CAB as 15 degrees, and Elise measures angle CBA as 20 degrees. They each use these measurements to calculate independent measurements on the height and compare. What do they find?
Note: To measure angles CAB and CBA, after recording the angle of max height – a vertical angle – students will move their inclinometer from the max height to a point on the ground directly underneath the max height and then measure a horizontal angle to the position of the other observer.
Deriving the Equation
Student Exploration
Lesson Three
Activity Title: Analyze a Thrust Curve
Overview: In the original project, thrust curves for rocket engines were modeled using piecewise linear functions. In the spreadsheet developed in the physics part of the project, time increments are being tracked in column B. A rocket engine thrust curve is a graphic depiction of how the engine thrust changes as a function of time. Modeling the curve using piecewise linear functions allows evaluation of thrust at any desired time, by breaking the curve into smaller time increments. In this lesson, the student will use this process to divide the thrust curve – after first defining it as a piecewise linear function – into finer time increments. The second use of this process will be to gain greater accuracy in the degradation of propellant mass over time. This is addressed in Lesson Four.
Learning Objective: Students will use piecewise linear functions and sequences to estimate instantaneous values of a construct that varies over time.
Prerequisite: Advanced algebra or precalculus. This lesson requires knowledge of piecewise functions.
Time Required: 80 minutes
Procedure: One of the first steps to using the numerical simulator (spreadsheet) within the basic lesson plan is to estimate the instantaneous value for thrust at each defined point in time. This is accomplished by analyzing the curve through piecewise linear functions. It is important to note that the provided thrust curves are not exact, so small errors students incur in defining piecewise functions will not significantly alter their results.
Begin by guiding students through the following example.
A thrust curve for an Estes rocket engine model B4.
The piecewise linear concept would seem to be the simplest way to estimate instantaneous values of thrust. The entire curve above could be approximated by breaking it into several linear regions between the following time intervals: [0, 0.12), [0.12, 0.20), [0.20, 0.24), [0.24, 0.30), [0.30, 1.0), [1.0, 1.02]. The regions where the slope of the thrust curve is changing rapidly should be broken down into smaller time increments to more accurately approximate it with a series of straight lines. The Table below summarizes the time and thrust values in these intervals
Time  Thrust 

0.0  0 
0.12  13 
0.20  6.0 
0.24  4.0 
0.30  3.5 
1.0  3.5 
1.02  0.0 
Plotting these points and connecting them with straight lines produces the following graph, which is very similar to the Estes B4 engine curve.
The most effective way to estimate thrust values at specific points in time would be to use the concept of a linear sequence within each of the identified linear regions. Two examples are shown below for the first two linear regions [0, 0.12) and [0.12, 0.20) in the curve above. The spreadsheet uses a time increment of 0.01 seconds. So we have to be able to evaluate thrust during the entire burn at every 0.01 second time increment.
The number of common differences (or the number of time increments) within each interval is found by dividing the change in time across the entire interval divided by the time increment
\begin{align*}N=\frac{t_ft_i}{\Delta t}\end{align*}
Common difference (change in thrust in each time increment) is found by dividing the change in thrust across the entire interval by the number of common differences (number of time increments in the interval).
\begin{align*}D=\frac{T_fT_i}{N}\end{align*}
Where
\begin{align*}N\end{align*} number of common differences (number of time increments in the interval)
\begin{align*}t_f\end{align*} final time in the interval
\begin{align*}t_i\end{align*} initial time in the interval
\begin{align*}\Delta t\end{align*} time increment
\begin{align*}T_f\end{align*} thrust at the end of the interval
\begin{align*}T_i\end{align*} thrust at the start of the interval
\begin{align*}D\end{align*} common difference (thrust change in each time increment)
At each time increment, thrust can be found by adding \begin{align*}D\end{align*} to the value of thrust from the last time increment.
\begin{align*}T_2 &= T_1 + D\\ T_3 &= T_2 + D\end{align*}
... and so on.
For region 1, time interval is [0, 0.12), and thrust values are [0, 13). On the spreadsheet, the time values begin at zero and advance by increments of 0.01. The total time difference in this interval (0.12 minus 0) divided by the time increment (0.01 seconds) will give the number of time intervals (common differences) in this region needed to advance a sequence from the first thrust value to the last on the interval.
\begin{align*}N &= \frac{0.120}{0.01}=12\\ D &= \frac{130}{12}1.08 \overline{3}\\ T_1 &= 0.0\\ T_2 &= 1.083\\ T_3 &= 2.167\\ T_4 &= 3.250\end{align*}
.. and so on.
For region 2, time interval is [0.12, 0.20) and thrust values are [13, 6).
\begin{align*}N &= \frac{0.200.12}{0.01}=8\\ D &= \frac{613}{8}=0.875\\ T_1 &= 13\\ T_2 &= 12.125\\ T_3 &= 11.250\\ T_4 &= 10.375\end{align*}
.. and so on.
The process in repeated for the remaining intervals. Note that the thrust values on time interval [0.30, 1.0) seem to be constant at an approximate value of 3.5. The thrust drops to zero at 1.02 seconds. The red dots below show all the points that were computed using the process described.
Distribute the student exploration and ask students to work through it.
Solution for the C11 Engine (which students will develop in the "Student Exploration"): Solutions may vary, depending on how accurately students wish to define piecewise functions.
For the equations in Table below,
 independent variable time = t
 dependent variable thurst = Th
Time Interval  Thrust Interval  Slope  Known Point x & y Coordinates  Equation 

0  .26  0  22  22/.26 
x: 0 y: 0 
Th = (22/.26)t 
.26  .36  22 – 10  12/.1 
x: .26 y: 22 
Th = (12/.1)t + 22 
.36  .90  10 – 5  5/.54 
x: .36 y: 10 
Th = (5/.54)t + 10 
.90 – .92  5 – 0  5/.02 
x: .90 y: 5 
Th = (5/.02)t + 5 
Using the concept of sequences with a timestep of 0.02:
For time interval 0 – .26,

 Number of common differences = .26/.02 = 13; D = 22/13

 Thrust values: 0, 1.69, 3.38, 5.08, 6.77, 8.46, 10.15, 11.85, 13.53, 15.23, 16.92, 18.62, 20.31
For time interval .26 – .36

 Number of common differences = .10/.02 = 5; D = 12/5 = 2.4

 Thrust values: 22, 19.6, 17.2, 14.8, 12.4
For time interval .36 – .90

 Number of common differences = .54/.02 = 27; D = 5/27

 Thrust values: 10, 9.81, 9.62, 9.44, 9.26, 9.07, 8.89, 8.70, 8.52, 8.33, 8.15, 7.96, 7.77, 7.59, 7.41, 7.22, 7.04, 6.85, 6.67, 6.48, 6.30, 6.11, 5.93, 5.74, 5.56, 5.37, 5.19
For time interval .90 – .92

 Number of common differences = .02/.02 = 1; D = 5/1 = 5

 Thrust values: 5, 0
This numerical processing would be easier to do in a spreadsheet. In the application developed at NASA, the author used a computer program (Octave, which is a hybrid of Matlab) to create a matrix that was entered into the spreadsheet.
Determine The Thrust Values For A Specific Engine
Student Exploration
Use the thrust curve of a C11 engine (shown below) to estimate the instantaneous thrust values for a sequence of time values.
Thrust profile for Estes CII Engine (Source: http://www2.estesrockets.com/pdf/Estes_TimeThrust_Curves.pdf)
Use the concepts of a piecewise linear function and sequences to come up with the values for thrust at a time interval of your choice.
Lesson Four
Activity Title: Calculate a More Accurate Mass: An Error Analysis
Overview: The original project was conceived as a capstone project for AP Physics and AP Calculus students. In Lesson Three, students created a spreadsheet that recorded thrust values (for engine propellant) for small time increments. In this lesson they extend that spreadsheet to account for mass degradation over the time interval.
Learning Objective: Students will explore methods to increase the accuracy of estimations, including how calculus is used by engineers to solve differential equations by numeric means.
Prerequisite: Calculus. This lesson requires knowledge of piecewise functions and integral calculus as it is used to determine area.
Time Required: 80 minutes
Procedure: Within the basic lesson plan, students assume the mass of the propellant degrades linearly. That is actually not the case, though; in fact, the mass depletion is proportional to the thrust. As the engine thrust changes, the rate at which the propellant is being ejected changes proportionally as well. In the world of rocketry, the constant of proportionality is called specific impulse, represented by the symbol \begin{align*}I_{sp}\end{align*}. The procedure below shows the derivation of \begin{align*}I_{sp}\end{align*} and how it is calculated for each specific engine. Quantities needed to calculate \begin{align*}I_{sp}\end{align*} are total mass depleted from the engine, and the thrust profile. Each engine will have a unique value for \begin{align*}I_{sp}\end{align*}.
\begin{align*}m_i\end{align*} initial mass
\begin{align*}m_f\end{align*} final mass
\begin{align*}T\end{align*} thrust at each instance of time
\begin{align*}t\end{align*} time
\begin{align*}g\end{align*} acceleration of gravity at Earth sea level, 9.81 m/s^{2}
\begin{align*}I_{sp}\end{align*} specific impulse
By definition, the rate of depletion of propellant is proportional to thrust. Note that thrust is a function of time as it changes with time.
\begin{align*}T(t) \propto \dot m\end{align*}
The right side is multiplied by standard Earth’s gravity field constant at sea level \begin{align*}g\end{align*} (9.81 m/s^{2}). The proportionality relationship should still hold.
\begin{align*}T(t) \propto \dot m g\end{align*}
As mentioned above, \begin{align*}I_{sp}\end{align*} is the constant of proportionality. The minus sign is included so that the value of \begin{align*}I_{sp}\end{align*} will be positive.
\begin{align*}T(t)=I_{sp} \dot m g\end{align*}
Let’s integrate both side of the equation from start of the engine burn to the end.
\begin{align*}\int\limits_{0}^{t_f} T(t)dt=\int\limits_{0}^{t_f} I_{sp} \frac{dm}{dt}g dt\end{align*}
\begin{align*}I_{sp}\end{align*} and \begin{align*}g\end{align*} are both constants and are moved out of the integral, and \begin{align*}dt\end{align*} terms cancel out of the right hand side. So the integral on the right hand side is now in terms of \begin{align*}dm\end{align*}.
\begin{align*}\int\limits_{0}^{t_f} T(t)dt &=I_{sp}g \int\limits_{0}^{t_f} dm\\ \int\limits_{0}^{t_f} T(t)dt &= I_{sp}g(m_ff_i)\\ \int\limits_{0}^{t_f} T(t)dt &= I_{sp}g(m_if_f)\end{align*}
Finally the specific impulse \begin{align*}I_{sp}\end{align*} is found using the equation below. The \begin{align*}m_im_f\end{align*} quantity is the amount of propellant mass expelled during the engine burn, and is computed by weighting the engine before and after the burn, and subtracting the two values. The propellant mass is given in specifications for rocket engines by the engine manufacturer. The integral \begin{align*}\int\limits_{0}^{t_f}T(t)dt\end{align*} is the area under the thrust curve and can easily be calculated.
\begin{align*}I_{sp}=\frac{1}{g(m_if_f)} \int\limits_{0}^{t_f}T(t)dt\end{align*}
Author Note: the notation \begin{align*}\dot m\end{align*}, which is more commonly seen by high school teachers of calculus as \begin{align*}m^{\prime}\end{align*}, is the notation the scientists we worked with at NASA used. I thought it interesting enough to use here. I’d never seen it before.
Also, it is important to note that the mass flow rate will be negative and that is why the calculation for \begin{align*}I_{sp}\end{align*} must include a negative sign, in order to make the value of \begin{align*}I_{sp}\end{align*} a positive quantity.
\begin{align*}I_{sp}\end{align*} can now be used to come up with a mass depletion model for the engine. Negative sign means the mass is decreasing.
\begin{align*}\dot m = \frac{T(t)}{I_{sp}g}\end{align*}
In each time step, it is assumed that mass flow rate is constant, so the amount of mass lost is the mass depletion rate above times the integration time step.
\begin{align*}\Delta m=\frac{T(t)}{I_{sp}g} \Delta t\end{align*}
Note: Physics students should recognize that the final term \begin{align*}\int\limits Tdt\end{align*} as the area under the thrust curve or total impulse of the engine.
Using the idea of Euler’s method (taking a continuous process and breaking it into small incremental pieces, thus changing it into a discrete process) is employed to get:
\begin{align*}new \ mass &= prior \ mass \left(\frac{thrust * (time  time \ prior)}{gI_{sp}}\right)\\ new \ mass &= prior \ mass  \left(\frac{thrust * \Delta t}{gI_{sp}}\right)\end{align*}
The "thrust" in the equation above is the incremental value that was calculated in the "Analyzing a Thrust Curve" extension.
In the numerical simulator being used in the basic lesson, students should already have columns for time and thrust in small, incremental pieces. A good approximation of the total thrust is the sum of the product of these incremental pieces. An excellent discussion can be generated about discrete approximations for continuous processes (see note at the end of this explanation). Once this idea is imparted, it is also important for the students to see that a Riemann sum is being applied as well
\begin{align*}\frac{1}{I_{sp}g} \int\limits_{0}^{t_f} T(t)dt=\frac{1}{I_{sp}g} \sum^{n}_{i=1} Thrust_n * \Delta t\end{align*}
and that a continuous process is being modeled by a discrete process. The result of this calculation can become another column in the spreadsheet. Then, through an iterative process, the change in propellant mass can be calculated. From there, the physics takes over.
Once students have corrected the mass column in their spreadsheet, they should compare the calculated maximum height using this new series of mass values with the calculated maximum height using the linear degradation of mass model. These two values should be different and a percent error should be calculated.
Adding yet another dimension to the error analysis, students can then account for the mass of the delay and tracking layer. This can be approximated by looking at the engine specifications: for example, an Estes B42 rocket engine has an initial mass of 19.8 grams, and a B44 has an initial mass of 21.0 grams. The propellant masses for both engines are the same, so the 1.2 gram difference must be the approximate mass of the delay and tracking layer. Assume the degradation is linear over the time period of the delay and account for this change in mass on the spreadsheet. Compare the maximum height calculated after this change with the two maximum height values compared in the previous paragraph; students should find the percent error is smaller. They should then do another error analysis.
Note about discrete approximations for continuous processes: In the author’s discussions with engineers and scientists at the NASA Langley Research Center, this topic arose. As it turns out, this method is very useful in the modeling and simulation process. The mathematics of many of these engineering models is based on differential equations that either cannot be solved explicitly or have complicated explicit solutions. Thus, they can often only be solved using numerical methods. The concept of using small, discrete increments provides solutions to problems that, while not error free, could not easily be done by other means. Euler method is the simplest, most crude method to perform numerical integration. However, it was chosen for this lesson plan because it provides an easytounderstand introduction into the concept of numerical integration. There are many other numerical integration schemes that are more complex but produce less error. One such numerical integration method used often is the RungeKutta method. The fourth order RungeKutta method is commonly used as a good trade between accuracy and computer processing time. Writeups are available in Wikipedia among other places.^{*} Whenever possible, solutions are always compared to empirical realities in order to refine the process.
^{*} Thanks to Jim Batterson, NASA Educational Consultant, for the information about RungeKutta. This method was also mentioned in several conversations with mentor Behzad Raiszadeh.
Calculate A More Accurate Mass: An Error Analysis
Student Exploration
In order to estimate the incremental change in propellant, the first quantity that must be calculated is specific impulse, or \begin{align*}I_{sp}\end{align*}, which is the thrust per unit mass expelled.
The variable set:
\begin{align*}m_i\end{align*} initial mass
\begin{align*}m_f\end{align*} final mass
\begin{align*}T\end{align*} thrust
\begin{align*}t\end{align*} time
\begin{align*}\dot m=\frac{dm}{dt}\end{align*} mass flow rate
\begin{align*}\dot w =\frac{dw}{dt}=\dot m g\end{align*} where \begin{align*}w\end{align*} is weight
\begin{align*}I_{sp}\end{align*} specific impulse
Questions:
1. If you are unfamiliar with the notation \begin{align*}\dot m\end{align*}, what do you think it means?
\begin{align*} \\ \\ \\ \end{align*}
2. Using the notation referenced in question 1, what is \begin{align*}\frac{d}{dt} \dot m\end{align*}?
\begin{align*} \\ \\ \\ \end{align*}
3. What is the sign of \begin{align*}\dot m\end{align*} (for rocket engines)?
\begin{align*} \\ \\ \\ \end{align*}
Now to the derivation of a formula for \begin{align*}I_{sp}\end{align*}.
4. Begin with:
\begin{align*}I_{sp} &=\frac{T}{\dot w}=\frac{T}{\dot m g}\\ I_{sp} \dot m g &= T\end{align*}
Integrate both sides (remember, \begin{align*}\dot m =\frac{dm}{dt}\end{align*} and consider \begin{align*}g I_{sp}\end{align*} a constant)
\begin{align*} \\ \\ \\ \end{align*}
5. Now show that:
\begin{align*}I_{sp}=\frac{1}{g(m_im_f)} \int\limits_{0}^{t_f} Tdt\end{align*}
and give an explanation of how you arrived at this.
The constant of proportionality for \begin{align*}I_{sp}\end{align*} is the term
\begin{align*}\frac{1}{g(m_im_f)}\end{align*}
which has a calculable value for any given engine, where \begin{align*}(m_im_f)\end{align*} is the change in propellant mass (in specifications for rocket engines, propellant mass is given).
\begin{align*} \\ \\ \\ \end{align*}
6. Calculate the value for \begin{align*}I_{sp}\end{align*} for the engine you are using.
\begin{align*} \\ \\ \\ \end{align*}
7. Once you have your \begin{align*}I_{sp}\end{align*} value, examine what we do next (the next three lines):
\begin{align*}I_{sp} &= \frac{1}{g(m_im_f)} \int\limits_{0}^{t_f} Tdt\\ (m_im_f) &= \frac{1}{gI_{sp}} \int\limits_{0}^{t_f} Tdt\\ new \ mass &= prior \ mass\left(\frac{thrust * (time  time \ prior)}{gI_{sp}}\right)\\ new \ mass &= prior \ mass\left(\frac{thrust * \Delta t}{gI_{sp}}\right)\end{align*}
The "thrust" in the equation above is the incremental value that was calculated in the thrust curve exploration.
Now, using your spreadsheet, the incremental change in mass can be calculated for each timestep. If your time increment is constant, which it should be, and you know that \begin{align*}g I_{sp}\end{align*} is constant, you can see how the change in mass is proportional to the change in thrust. The method shown here, using a numerical calculation with small, incremental steps to solve a differential equation, is known as Euler’s method.
8. Error analysis. There are two levels of analysis that can be performed.
Part 1: Without the concept of variably changing mass, most physics classes will assume the propellant depletes linearly. So, your task here is to find the maximum height when the mass is changing variably (the analysis above) and then find the maximum height when it is changing linearly. Calculate the percent error between these two values.
Part 2: The variably changing mass analysis does not account for the mass of the delay and tracking layer. If it is possible to determine this, do so and add this to the propellant mass. Remember that the mass of the delay and tracking layer does not begin to deplete until the propellant is exhausted. It will be important to account for this in the proper place in your spreadsheet. Once this is done, you can compare the maximum heights again.
Lesson Five
Activity Title: Computer Programming through Octave
Overview:
In the base lesson plan, students created a numerical simulator within a spreadsheet. A numerical simulator such as this can be created using computer programming languages as well. In this extension, a computer program written in a language called Octave is used to interpolate values for thrust. The concept is the same as what students did in the thrust curve exploration, but the computer is doing the work for them. Octave is a free download but has little support. However, it is a hybrid of another language called Matlab. Matlab is not free, but it does have online support. Matlab is used by engineers in industry and can be very expensive per seat license. However, Matlab has special deals for K–12 educators, and, while still not cheap, it is considerably less expensive than the standard industry price. We found that code written in Octave using Matlab support was relatively easy.
The use of a computer program to generate a matrix of time / thrust values may resonate with some teachers and students. The code for the Octave program that does this is given below. The output of the program below is a text file, "thrust_profile.txt" on the desktop (or in the folder where the program resides). "thrust_profile.txt" can be loaded into the spreadsheet for use in the simulation.
X = [0.00 0.22 0.28 0.66 0.72 100.0]; Y = [0.00 10.00 5.00 3.00 0.00 0.0]; time = 0:0.01:1.1; thrust = interp1(X,Y,time); fp = fopen('thrust_profile.txt','w'); for n=1:length(time) fprintf(fp,'%8.2f %8.4f\n',time(n),thrust(n)); end fclose(fp);
Students would have to change the \begin{align*}X\end{align*} and \begin{align*}Y\end{align*} values (these are just endpoints of the piecewise linear functions, there can be as many as needed) and the endtime (currently set at 1.1 in the third line of code).
The script above can be written in a regular text editor and saved with a ".m" extension. Let’s say we named the file containing the script above thrust_curve.m.
Starting Octave brings up the screen below. Octave accepts user input, line by line, at the command prompt. In the screen shot below, the command prompt is the last line, octave3.4.0:1>. User input is entered immediately after the > sign. In order to run the program above, the user needs to type "thrust_curve" at the command prompt and press "Enter." Octave then processes the contents of the file sequentially and generates the output file, thrust_profile.txt. Octave is only able to see files that are at the current working directory. thrust_curve.m will not run if it is not in the current working directory of Octave. The current working directory can be changed by issuing the cd directory_name command at the command prompt. directory_name is the name of the folder (or directory) where thrust_curve.m resides.
Octave allows users to copy and paste commands directly on the Octave command prompt. What the author did was to start Octave and type the command "cd Desktop" (no quotes) on the command prompt of Octave to change the current working directory to "Desktop". On Macs, "Desktop" has to be typed with the first character in capital letter. The author then copied the code above into the clipboard and pasted it on the Octave window. Octave ran the contents of the program line by line and produced the output file "thrust_profile.txt" on the "Desktop".
A screenshot of the computer follows as the example of what happened. Note that now there is a file that has been generated called "thrust_profile.txt."