Mathematical Numerical methods

mathematical modeling numerical methods and complexes of programs and mathematical analysis and numerical methods for science and technology
GregDeamons Profile Pic
GregDeamons,New Zealand,Professional
Published Date:03-08-2017
Your Website URL(Optional)
cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 4 1 Mathematical Modeling, Numerical Methods, and Problem Solving CHAPTER OBJECTIVES The primary objective of this chapter is to provide you with a concrete idea of what numerical methods are and how they relate to engineering and scientific problem solving. Specific objectives and topics covered are • Learning how mathematical models can be formulated on the basis of scientific principles to simulate the behavior of a simple physical system. • Understanding how numerical methods afford a means to generate solutions in a manner that can be implemented on a digital computer. • Understanding the different types of conservation laws that lie beneath the models used in the various engineering disciplines and appreciating the difference between steady-state and dynamic solutions of these models. • Learning about the different types of numerical methods we will cover in this book. YOU’VE GOT A PROBLEM uppose that a bungee-jumping company hires you. You’re given the task of predict- ing the velocity of a jumper (Fig. 1.1) as a function of time during the free-fall part S of the jump. This information will be used as part of a larger analysis to determine the length and required strength of the bungee cord for jumpers of different mass. You know from your studies of physics that the acceleration should be equal to the ratio of the force to the mass (Newton’s second law). Based on this insight and your knowledge 4cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 5 1.1 A SIMPLE MATHEMATICAL MODEL 5 of physics and fluid mechanics, you develop the following mathematical model for the rate Upward force due to air of change of velocity with respect to time, resistance dv c d 2 = g − v dt m where v = downward vertical velocity (m/s), t = time (s), g = the acceleration due to 2 ∼ gravity (= 9.81 m/s ), c = a lumped drag coefficient (kg/m), and m = the jumper’s d mass (kg). The drag coefficient is called “lumped” because its magnitude depends on fac- tors such as the jumper’s area and the fluid density (see Sec. 1.4). Because this is a differential equation, you know that calculus might be used to obtain an analytical or exact solution for v as a function of t. However, in the following pages, we will illustrate an alternative solution approach. This will involve developing a computer- oriented numerical or approximate solution. Aside from showing you how the computer can be used to solve this particular prob- lem, our more general objective will be to illustrate (a) what numerical methods are and Downward force due (b) how they figure in engineering and scientific problem solving. In so doing, we will also to gravity show how mathematical models figure prominently in the way engineers and scientists use numerical methods in their work. FIGURE 1.1 Forces acting on a free-falling bungee jumper. 1.1 A SIMPLE MATHEMATICAL MODEL A mathematical model can be broadly defined as a formulation or equation that expresses the essential features of a physical system or process in mathematical terms. In a very gen- eral sense, it can be represented as a functional relationship of the form   Dependent independent forcing = f , parameters, (1.1) variable variables functions where the dependent variable is a characteristic that typically reflects the behavior or state of the system; the independent variables are usually dimensions, such as time and space, along which the system’s behavior is being determined; the parameters are reflective of the system’s properties or composition; and the forcing functions are external influences acting upon it. The actual mathematical expression of Eq. (1.1) can range from a simple algebraic relationship to large complicated sets of differential equations. For example, on the basis of his observations, Newton formulated his second law of motion, which states that the time rate of change of momentum of a body is equal to the resultant force acting on it. The math- ematical expression, or model, of the second law is the well-known equation F = ma (1.2) 2 where F is the net force acting on the body (N, or kg m/s ), m is the mass of the object (kg), 2 and a is its acceleration (m/s ).cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 6 6 MATHEMATICAL MODELING, NUMERICAL METHODS, AND PROBLEM SOLVING The second law can be recast in the format of Eq. (1.1) by merely dividing both sides by m to give F a = (1.3) m where a is the dependent variable reflecting the system’s behavior, F is the forcing func- tion, and m is a parameter. Note that for this simple case there is no independent variable because we are not yet predicting how acceleration varies in time or space. Equation (1.3) has a number of characteristics that are typical of mathematical models of the physical world. • It describes a natural process or system in mathematical terms. • It represents an idealization and simplification of reality. That is, the model ignores neg- ligible details of the natural process and focuses on its essential manifestations. Thus, the second law does not include the effects of relativity that are of minimal importance when applied to objects and forces that interact on or about the earth’s surface at veloc- ities and on scales visible to humans. • Finally, it yields reproducible results and, consequently, can be used for predictive pur- poses. For example, if the force on an object and its mass are known, Eq. (1.3) can be used to compute acceleration. Because of its simple algebraic form, the solution of Eq. (1.2) was obtained easily. However, other mathematical models of physical phenomena may be much more complex, and either cannot be solved exactly or require more sophisticated mathematical techniques than simple algebra for their solution. To illustrate a more complex model of this kind, Newton’s second law can be used to determine the terminal velocity of a free-falling body near the earth’s surface. Our falling body will be a bungee jumper (Fig. 1.1). For this case, a model can be derived by expressing the acceleration as the time rate of change of the velocity (dv/dt) and substituting it into Eq. (1.3) to yield dv F = (1.4) dt m where v is velocity (in meters per second). Thus, the rate of change of the velocity is equal to the net force acting on the body normalized to its mass. If the net force is positive, the object will accelerate. If it is negative, the object will decelerate. If the net force is zero, the object’s velocity will remain at a constant level. Next, we will express the net force in terms of measurable variables and parameters. For a body falling within the vicinity of the earth, the net force is composed of two oppos- ing forces: the downward pull of gravity F and the upward force of air resistance F D U (Fig. 1.1): F = F + F (1.5) D U If force in the downward direction is assigned a positive sign, the second law can be used to formulate the force due to gravity as F = mg (1.6) D 2 where g is the acceleration due to gravity (9.81 m/s ).cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 7 1.1 A SIMPLE MATHEMATICAL MODEL 7 Air resistance can be formulated in a variety of ways. Knowledge from the science of fluid mechanics suggests that a good first approximation would be to assume that it is pro- portional to the square of the velocity, 2 F =−c v (1.7) U d where c is a proportionality constant called the lumped drag coefficient (kg/m). Thus, the d greater the fall velocity, the greater the upward force due to air resistance. The parameter c accounts for properties of the falling object, such as shape or surface roughness, that af- d fect air resistance. For the present case, c might be a function of the type of clothing or the d orientation used by the jumper during free fall. The net force is the difference between the downward and upward force. Therefore, Eqs. (1.4) through (1.7) can be combined to yield dv c d 2 = g − v (1.8) dt m Equation (1.8) is a model that relates the acceleration of a falling object to the forces acting on it. It is a differential equation because it is written in terms of the differential rate of change (dv/dt) of the variable that we are interested in predicting. However, in contrast to the solution of Newton’s second law in Eq. (1.3), the exact solution of Eq. (1.8) for the velocity of the jumper cannot be obtained using simple algebraic manipulation. Rather, more advanced techniques such as those of calculus must be applied to obtain an exact or analytical solution. For example, if the jumper is initially at rest (v = 0 at t = 0), calculus can be used to solve Eq. (1.8) for    gm gc d v(t) = tanh t (1.9) c m d 1 where tanh is the hyperbolic tangent that can be either computed directly or via the more elementary exponential function as in x −x e − e tanh x = (1.10) x −x e + e Note that Eq. (1.9) is cast in the general form of Eq. (1.1) where v(t) is the dependent variable, t is the independent variable, c and m are parameters, and g is the forcing function. d EXAMPLE 1.1 Analytical Solution to the Bungee Jumper Problem Problem Statement. A bungee jumper with a mass of 68.1 kg leaps from a stationary hot air balloon. Use Eq. (1.9) to compute velocity for the first 12 s of free fall. Also determine the terminal velocity that will be attained for an infinitely long cord (or alternatively, the jumpmaster is having a particularly bad day). Use a drag coefficient of 0.25 kg/m. 1 MATLAB allows direct calculation of the hyperbolic tangent via the built-in function tanh(x).cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 8 8 MATHEMATICAL MODELING, NUMERICAL METHODS, AND PROBLEM SOLVING Solution. Inserting the parameters into Eq. (1.9) yields     9.81(68.1) 9.81(0.25) v(t) = tanh t = 51.6938 tanh(0.18977t) 0.25 68.1 which can be used to compute t, s v, m/s 00 2 18.7292 4 33.1118 6 42.0762 8 46.9575 10 49.4214 12 50.6175 ∞ 51.6938 According to the model, the jumper accelerates rapidly (Fig. 1.2). A velocity of 49.4214 m/s (about 110 mi/hr) is attained after 10 s. Note also that after a sufficiently long FIGURE 1.2 The analytical solution for the bungee jumper problem as computed in Example 1.1. Velocity increases with time and asymptotically approaches a terminal velocity. 60 Terminal velocity 40 20 0 048 12 t, s u, m/scha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 9 1.1 A SIMPLE MATHEMATICAL MODEL 9 time, a constant velocity, called the terminal velocity, of 51.6983 m/s (115.6 mi/hr) is reached. This velocity is constant because, eventually, the force of gravity will be in bal- ance with the air resistance. Thus, the net force is zero and acceleration has ceased. Equation (1.9) is called an analytical or closed-form solution because it exactly satis- fies the original differential equation. Unfortunately, there are many mathematical models that cannot be solved exactly. In many of these cases, the only alternative is to develop a numerical solution that approximates the exact solution. Numerical methods are those in which the mathematical problem is reformulated so it can be solved by arithmetic operations. This can be illustrated for Eq. (1.8) by realizing that the time rate of change of velocity can be approximated by (Fig. 1.3): dv v v(t ) − v(t ) i+1 i ∼ = = (1.11) dt t t − t i+1 i where v and t are differences in velocity and time computed over finite intervals, v(t ) i is velocity at an initial time t , and v(t ) is velocity at some later time t . Note that i i+1 i+1 ∼ dv/dt = v/t is approximate because t is finite. Remember from calculus that dv v = lim dt t→0 t Equation (1.11) represents the reverse process. FIGURE 1.3 The use of a finite difference to approximate the first derivative of v with respect to t. v(t ) i1 True slope dv/dt v Approximate slope v(t )  v(t ) v(t ) v i1 i i  t t  t i1 i t t t i i1 tcha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 10 10 MATHEMATICAL MODELING, NUMERICAL METHODS, AND PROBLEM SOLVING Equation (1.11) is called a finite-difference approximation of the derivative at time t . i It can be substituted into Eq. (1.8) to give v(t ) − v(t ) c i+1 i d 2 = g − v(t ) i t − t m i+1 i This equation can then be rearranged to yield   c d 2 v(t ) = v(t ) + g − v(t ) (t − t ) (1.12) i+1 i i i+1 i m Notice that the term in brackets is the right-hand side of the differential equation itself Eq. (1.8). That is, it provides a means to compute the rate of change or slope of v. Thus, the equation can be rewritten more concisely as dv i v = v + t (1.13) i+1 i dt where the nomenclature v designates velocity at time t and t = t − t . i i i+1 i We can now see that the differential equation has been transformed into an equation that can be used to determine the velocity algebraically at t using the slope and previous val- i+1 ues ofv and t. If you are given an initial value for velocity at some time t , you can easily com- i pute velocity at a later time t . This new value of velocity at t can in turn be employed to i+1 i+1 extend the computation to velocity at t and so on. Thus at any time along the way, i+2 New value = old value + slope × step size This approach is formally called Euler’s method. We’ll discuss it in more detail when we turn to differential equations later in this book. EXAMPLE 1.2 Numerical Solution to the Bungee Jumper Problem Problem Statement. Perform the same computation as in Example 1.1 but use Eq. (1.12) to compute velocity with Euler’s method. Employ a step size of 2 s for the calculation. Solution. At the start of the computation (t = 0), the velocity of the jumper is zero. 0 Using this information and the parameter values from Example 1.1, Eq. (1.12) can be used to compute velocity at t = 2 s: 1   0.25 2 v = 0 + 9.81 − (0) × 2 = 19.62 m/s 68.1 For the next interval (from t = 2 to 4 s), the computation is repeated, with the result   0.25 2 v = 19.62 + 9.81 − (19.62) × 2 = 36.4137 m/s 68.1cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 11 1.1 A SIMPLE MATHEMATICAL MODEL 11 60 Terminal velocity Approximate, numerical solution 40 Exact, analytical solution 20 0 48 12 0 t, s FIGURE 1.4 Comparison of the numerical and analytical solutions for the bungee jumper problem. The calculation is continued in a similar fashion to obtain additional values: t, s v, m/s 00 2 19.6200 4 36.4137 6 46.2983 8 50.1802 10 51.3123 12 51.6008 ∞ 51.6938 The results are plotted in Fig. 1.4 along with the exact solution. We can see that the nu- merical method captures the essential features of the exact solution. However, because we have employed straight-line segments to approximate a continuously curving function, there is some discrepancy between the two results. One way to minimize such discrepan- cies is to use a smaller step size. For example, applying Eq. (1.12) at 1-s intervals results in a smaller error, as the straight-line segments track closer to the true solution. Using hand calculations, the effort associated with using smaller and smaller step sizes would make such numerical solutions impractical. However, with the aid of the computer, large num- bers of calculations can be performed easily. Thus, you can accurately model the velocity of the jumper without having to solve the differential equation exactly. v, m/scha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 12 12 MATHEMATICAL MODELING, NUMERICAL METHODS, AND PROBLEM SOLVING As in Example 1.2, a computational price must be paid for a more accurate numerical result. Each halving of the step size to attain more accuracy leads to a doubling of the num- ber of computations. Thus, we see that there is a trade-off between accuracy and computa- tional effort. Such trade-offs figure prominently in numerical methods and constitute an important theme of this book. 1.2 CONSERVATION LAWS IN ENGINEERING AND SCIENCE Aside from Newton’s second law, there are other major organizing principles in science and engineering. Among the most important of these are the conservation laws. Although they form the basis for a variety of complicated and powerful mathematical models, the great conservation laws of science and engineering are conceptually easy to understand. They all boil down to Change = increases − decreases (1.14) This is precisely the format that we employed when using Newton’s law to develop a force balance for the bungee jumper Eq. (1.8). Although simple, Eq. (1.14) embodies one of the most fundamental ways in which conservation laws are used in engineering and science—that is, to predict changes with respect to time. We will give it a special name—the time-variable (or transient) computation. Aside from predicting changes, another way in which conservation laws are applied is for cases where change is nonexistent. If change is zero, Eq. (1.14) becomes Change = 0 = increases − decreases or Increases = decreases (1.15) Thus, if no change occurs, the increases and decreases must be in balance. This case, which is also given a special name—the steady-state calculation—has many applications in engi- neering and science. For example, for steady-state incompressible fluid flow in pipes, the flow into a junction must be balanced by flow going out, as in Flow in = o w out For the junction in Fig. 1.5, the balance can be used to compute that the flow out of the fourth pipe must be 60. For the bungee jumper, the steady-state condition would correspond to the case where the net force was zero or Eq. (1.8) with dv/dt = 0 2 mg = c v (1.16) d Thus, at steady state, the downward and upward forces are in balance and Eq. (1.16) can be solved for the terminal velocity  gm v = c d Although Eqs. (1.14) and (1.15) might appear trivially simple, they embody the two funda- mental ways that conservation laws are employed in engineering and science. As such, they will form an important part of our efforts in subsequent chapters to illustrate the connection between numerical methods and engineering and science.cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 13 1.3 NUMERICAL METHODS COVERED IN THIS BOOK 13 Pipe 2 Flow in  80 Pipe 1 Pipe 4 Flow in  100 Flow out  ? Pipe 3 Flow out  120 FIGURE 1.5 A flow balance for steady incompressible fluid flow at the junction of pipes. Table 1.1 summarizes some models and associated conservation laws that figure promi- nently in engineering. Many chemical engineering problems involve mass balances for reactors. The mass balance is derived from the conservation of mass. It specifies that the change of mass of a chemical in the reactor depends on the amount of mass flowing in minus the mass flowing out. Civil and mechanical engineers often focus on models developed from the conserva- tion of momentum. For civil engineering, force balances are utilized to analyze structures such as the simple truss in Table 1.1. The same principles are employed for the mechanical engineering case studies to analyze the transient up-and-down motion or vibrations of an automobile. Finally, electrical engineering studies employ both current and energy balances to model electric circuits. The current balance, which results from the conservation of charge, is simi- lar in spirit to the flow balance depicted in Fig. 1.5. Just as flow must balance at the junction of pipes, electric current must balance at the junction of electric wires. The energy balance specifies that the changes of voltage around any loop of the circuit must add up to zero. It should be noted that there are many other branches of engineering beyond chemical, civil, electrical, and mechanical. Many of these are related to the Big Four. For example, chem- icalengineeringskillsareusedextensivelyinareassuchasenvironmental,petroleum,andbio- medical engineering. Similarly, aerospace engineering has much in common with mechanical engineering. I will endeavor to include examples from these areas in the coming pages. 1.3 NUMERICAL METHODS COVERED IN THIS BOOK Euler’s method was chosen for this introductory chapter because it is typical of many other classes of numerical methods. In essence, most consist of recasting mathematical opera- tions into the simple kind of algebraic and logical operations compatible with digital com- puters. Figure 1.6 summarizes the major areas covered in this text. cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 14 14 MATHEMATICAL MODELING, NUMERICAL METHODS, AND PROBLEM SOLVING TABLE 1.1 Devices and types of balances that are commonly used in the four major areas of engineering. For each case, the conservation law on which the balance is based is specified. Field Device Organizing Principle Mathematical Expression Chemical Conservation Mass balance: engineering of mass Input Output Reactors Over a unit of time period mass  inputs  outputs F Civil Conservation Force balance: V engineering of momentum Structure F F H H F V At each node  horizontal forces (F )  0 H  vertical forces (F )  0 V Mechanical Conservation Force balance: Machine Upward force engineering of momentum x  0 Downward force 2 d x m  downward force  upward force 2 dt Electrical Conservation Current balance: i i 1 3 engineering of charge  For each node  i 2  current (i)  0 Circuit i R 1 1 Conservation Voltage balance: of energy i R  2 2 i R 3 3 Around each loop  emf’s   voltage drops for resistors  0     iR  0cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 15 1.3 NUMERICAL METHODS COVERED IN THIS BOOK 15 (a) Part 2 : Roots and optimization f(x) Roots: Solve for x so that f(x)  0 Roots Optimization: Solve for x so that f'(x)  0 x Optima (b) Part 3 : Linear algebraic equations x 2 Given the a’s and the b’s, solve for the x’s Solution a x  a x  b 11 1 12 2 1 a x  a x  b 21 1 22 2 2 x 1 (c) Part 4 : Curve fitting f(x) f(x) Interpolation Regression x x (d) Part 5 : Integration and differentiation dy/dx y Integration: Find the area under the curve Differentiation: Find the slope of the curve I x y (e) Part 6 : Differential equations Slope  f(t , y ) Given i i dy y   f(t, y) dt t solve for y as a function of t t y  y  f(t , y )t i1 i i i t FIGURE 1.6 Summary of the numerical methods covered in this book.cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 16 16 MATHEMATICAL MODELING, NUMERICAL METHODS, AND PROBLEM SOLVING Part Two deals with two related topics: root finding and optimization. As depicted in Fig. 1.6a, root location involves searching for the zeros of a function. In contrast, optimiza- tion involves determining a value or values of an independent variable that correspond to a “best” or optimal value of a function. Thus, as in Fig. 1.6a, optimization involves identify- ing maxima and minima. Although somewhat different approaches are used, root location and optimization both typically arise in design contexts. Part Three is devoted to solving systems of simultaneous linear algebraic equations (Fig. 1.6b). Such systems are similar in spirit to roots of equations in the sense that they are concerned with values that satisfy equations. However, in contrast to satisfying a single equation, a set of values is sought that simultaneously satisfies a set of linear algebraic equations. Such equations arise in a variety of problem contexts and in all disciplines of en- gineering and science. In particular, they originate in the mathematical modeling of large systems of interconnected elements such as structures, electric circuits, and fluid networks. However, they are also encountered in other areas of numerical methods such as curve fit- ting and differential equations. As an engineer or scientist, you will often have occasion to fit curves to data points. The techniques developed for this purpose can be divided into two general categories: regression and interpolation. As described in Part Four (Fig. 1.6c), regression is employed where there is a significant degree of error associated with the data. Experi- mental results are often of this kind. For these situations, the strategy is to derive a sin- gle curve that represents the general trend of the data without necessarily matching any individual points. In contrast, interpolation is used where the objective is to determine intermediate val- ues between relatively error-free data points. Such is usually the case for tabulated infor- mation. The strategy in such cases is to fit a curve directly through the data points and use the curve to predict the intermediate values. As depicted in Fig. 1.6d, Part Five is devoted to integration and differentiation. A physical interpretation of numerical integration is the determination of the area under a curve. Integration has many applications in engineering and science, ranging from the determination of the centroids of oddly shaped objects to the calculation of total quan- tities based on sets of discrete measurements. In addition, numerical integration formu- las play an important role in the solution of differential equations. Part Five also covers methods for numerical differentiation. As you know from your study of calculus, this involves the determination of a function’s slope or its rate of change. Finally, Part Six focuses on the solution of ordinary differential equations (Fig. 1.6e). Such equations are of great significance in all areas of engineering and science. This is because many physical laws are couched in terms of the rate of change of a quantity rather than the magnitude of the quantity itself. Examples range from population-forecasting models (rate of change of population) to the acceleration of a falling body (rate of change of velocity). Two types of problems are addressed: initial-value and boundary-value problems.cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 17 1.4 CASE STUDY 17 1.4 CASE STUDY IT’S A REAL DRAG Background. In our model of the free-falling bungee jumper, we assumed that drag depends on the square of velocity (Eq. 1.7). A more detailed representation, which was originally formulated by Lord Rayleigh, can be written as 1 2 F =− ρv AC v  (1.17) d d 2 3 where F  the drag force (N),   fluid density (kg/m ), A  the frontal area of the object d 2 on a plane perpendicular to the direction of motion (m ), C  a dimensionless drag coef- d ficient, and v   a unit vector indicating the direction of velocity. This relationship, which assumes turbulent conditions (i.e., a high Reynolds number), allows us to express the lumped drag coefficient from Eq. (1.7) in more fundamental terms as 1 C = ρ AC (1.18) d d 2 Thus, the lumped drag coefficient depends on the object’s area, the fluid’s density, and a dimensionless drag coefficient. The latter accounts for all the other factors that contribute to air resistance such as the object’s “roughness”. For example, a jumper wearing a baggy outfit will have a higher C than one wearing a sleek jumpsuit. d Note that for cases where velocity is very low, the flow regime around the object will be laminar and the relationship between the drag force and velocity becomes linear. This is referred to as Stokes drag. In developing our bungee jumper model, we assumed that the downward direction was positive. Thus, Eq. (1.7) is an accurate representation of Eq. (1.17), because v  1 and the drag force is negative. Hence, drag reduces velocity. But what happens if the jumper has an upward (i.e., negative) velocity? In this case, v –1 and Eq. (1.17) yields a positive drag force. Again, this is physically correct as the pos- itive drag force acts downward against the upward negative velocity. Unfortunately, for this case, Eq. (1.7) yields a negative drag force because it does not include the unit directional vector. In other words, by squaring the velocity, its sign and hence its direction is lost. Consequently, the model yields the physically unrealistic result that air resistance acts to accelerate an upward velocity In this case study, we will modify our model so that it works properly for both downward and upward velocities. We will test the modified model for the same case as Example 1.2, but with an initial value of v(0) 40 m/s. In addition, we will also illustrate how we can extend the numerical analysis to determine the jumper’s position. Solution. The following simple modification allows the sign to be incorporated into the drag force: 1 F =− ρvvAC (1.19) d d 2cha01102_ch01_001-023.qxd 12/18/10 1:56 PM Page 18 18 MATHEMATICAL MODELING, NUMERICAL METHODS, AND PROBLEM SOLVING 1.4 CASE STUDY continued or in terms of the lumped drag: F =−c vv (1.20) d d Thus, the differential equation to be solved is dv c d = g − vv (1.21) dt m In order to determine the jumper’s position, we recognize that distance travelled, x (m), is related to velocity by dx =−v (1.22) dt In contrast to velocity, this formulation assumes that upward distance is positive. In the same fashion as Eq. (1.12), this equation can be integrated numerically with Euler’s method: x = x − v(t )t (1.23) i+1 i i Assuming that the jumper’s initial position is defined as x(0)  0, and using the parame- ter values from Examples 1.1 and 1.2, the velocity and distance at t  2 s can be computed as   0.25 v(2)=−40 + 9.81 − (−40)(40) 2=−8.6326 m/s 68.1 x(2) = 0 − (−40)2 = 80 m Note that if we had used the incorrect drag formulation, the results would be 32.1274 m/s and 80 m. The computation can be repeated for the next interval (t  2 to 4 s):   0.25 v(4)=−8.6326 + 9.81 − (−8.6326)(8.6326) 2 = 11.5346 m/s 68.1 x(4) = 80 − (−8.6326)2 = 97.2651 m The incorrect drag formulation gives –20.0858 m/s and 144.2549 m. The calculation is continued and the results shown in Fig. 1.7 along with those obtained with the incorrect drag model. Notice that the correct formulation decelerates more rapidly because drag always diminishes the velocity. With time, both velocity solutions converge on the same terminal velocity because eventually both are directed downward in which case, Eq. (1.7) is correct. However, the impact on the height prediction is quite dramatic with the incorrect drag case resulting in a much higher trajectory. This case study demonstrates how important it is to have the correct physical model. In some cases, the solution will yield results that are clearly unrealistic. The current exam- ple is more insidious as there is no visual evidence that the incorrect solution is wrong. That is, the incorrect solution “looks” reasonable. cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 19 1.4 CASE STUDY 19 1.4 CASE STUDY continued (a) Velocity, m/s 60 40 Correct drag 20 0 48 12 t, s Incorrect drag –20 –40 (b) Height, m 200 Incorrect drag 100 0 48 12 t, s –100 Correct drag –200 FIGURE 1.7 Plots of (a) velocity and (b) height for the free-falling bungee jumper with an upward (negative) initial velocity generated with Euler’s method. Results for both the correct (Eq. 1 20) and incorrect (Eq. 1.7) drag formulations are displayed. x, m v, m/scha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 20 20 MATHEMATICAL MODELING, NUMERICAL METHODS, AND PROBLEM SOLVING PROBLEMS 1.1 Use calculus to verify that Eq. (1.9) is a solution of 1.5 Rather than the nonlinear relationship of Eq. (1.7), you Eq. (1.8) for the initial condition v(0)  0. might choose to model the upward force on the bungee 1.2 Use calculus to solve Eq. (1.21) for the case where the ini- jumper as a linear relationship: tial velocity is (a) positive and (b) negative. (c) Based on your  F =−c v U results for (a) and (b), perform the same computation as in Ex-  where c = a first-order drag coefficient (kg/s). ample 1.1 but with an initial velocity of 40 m/s. Compute (a) Using calculus, obtain the closed-form solution for the values of the velocity from t  0 to 12 s at intervals of 2 s. Note case where the jumper is initially at rest (v = 0 at t = 0). that for this case, the zero velocity occurs at t  3.470239 s. (b) Repeat the numerical calculation in Example 1.2 with 1.3 The following information is available for a bank account: the same initial condition and parameter values. Use a  value of 11.5 kg/s for c . Date Deposits Withdrawals Balance 1.6 For the free-falling bungee jumper with linear drag (Prob. 1.5), assume a first jumper is 70 kg and has a drag co- 5/1 1512.33 efficient of 12 kg/s. If a second jumper has a drag coefficient 220.13 327.26 6/1 of 15 kg/s and a mass of 80 kg, how long will it take her to 216.80 378.61 reach the same velocity jumper 1 reached in 9 s? 7/1 1.7 For the second-order drag model (Eq. 1.8), compute the 450.25 106.80 velocity of a free-falling parachutist using Euler’s method 8/1 for the case where m = 80 kg and c = 0.25 kg/m. Perform d 127.31 350.61 the calculation from t = 0 to 20 s with a step size of 1 s. Use 9/1 an initial condition that the parachutist has an upward veloc- ity of 20 m/s at t = 0. At t = 10 s, assume that the chute is instantaneously deployed so that the drag coefficient jumps Note that the money earns interest which is computed as to 1.5 kg/m. Interest = iB i 1.8 The amount of a uniformly distributed radioactive con- taminant contained in a closed reactor is measured by its where i  the interest rate expressed as a fraction per month, concentration c (becquerel/liter or Bq/L). The contaminant and B the initial balance at the beginning of the month. i decreases at a decay rate proportional to its concentration; (a) Use the conservation of cash to compute the balance on that is 6/1, 7/1, 8/1, and 9/1 if the interest rate is 1% per month Decay rate=−kc (i  0.01/month). Show each step in the computation. (b) Write a differential equation for the cash balance in the −1 where k is a constant with units of day . Therefore, accord- form ing to Eq. (1.14), a mass balance for the reactor can be written as dB = f D(t), W(t), i dt dc =−kc where t  time (months), D(t)  deposits as a function dt     of time (/month), W(t)  withdrawals as a function of change decrease = time (/month). For this case, assume that interest is in mass by decay compounded continuously; that is, interest  iB. (c) Use Euler’s method with a time step of 0.5 month to (a) Use Euler’s method to solve this equation from t = 0 to –1 simulate the balance. Assume that the deposits and with- 1 d with k = 0.175 d . Employ a step size of t = 0.1 d. drawals are applied uniformly over the month. The concentration at t = 0 is 100 Bq/L. (d) Develop a plot of balance versus time for (a) and (c). (b) Plot the solution on a semilog graph (i.e., ln c versus t) 1.4 Repeat Example 1.2. Compute the velocity to t = 12 s, and determine the slope. Interpret your results. with a step size of (a) 1 and (b) 0.5 s. Can you make any 1.9 A storage tank (Fig. P1.9) contains a liquid at depth y statement regarding the errors of the calculation based on the where y = 0 when the tank is half full. Liquid is withdrawn results? at a constant flow rate Q to meet demands. The contents arecha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 21 PROBLEMS 21 y y r top top Q in y 0 y out s Q out 1 0 FIGURE P1.11 FIGURE P1.9 2 2 resupplied at a sinusoidal rate 3Q sin (t). Equation (1.14) The liquid flows in at a sinusoidal rate of Q = 3 sin (t) and in can be written for this system as flows out according to 1.5 d(Ay) Q = 3(y − y ) y y 2 out out out = 3Q sin (t) − Q Q = 0 y ≤ y dt out out   change in = (in o w) − (out o w) 3 where flow has units of m /d and y = the elevation of the volume water surface above the bottom of the tank (m). Use Euler’s method to solve for the depth y from t = 0 to 10 d with a step or, since the surface area A is constant size of 0.5 d. The parameter values are r  2.5 m, y  4 m, top top and y  1 m. Assume that the level is initially below the out dy Q Q 2 outlet pipe with y(0)  0.8 m. = 3 sin (t) − dt A A 1.12 A group of 35 students attend a class in an insulated room which measures 11 m by 8 m by 3 m. Each student 3 takes up about 0.075 m and gives out about 80 W of heat Use Euler’s method to solve for the depth y from t = 0 to (1 W = 1 J/s). Calculate the air temperature rise during the first 10 d with a step size of 0.5 d. The parameter values are A = 2 3 20 minutes of the class if the room is completely sealed and in- 1250 m and Q = 450 m /d. Assume that the initial condition sulated. Assume the heat capacity C for air is 0.718 kJ/(kg K). is y = 0. v Assume air is an ideal gas at 20°C and 101.325 kPa. Note 1.10 For the same storage tank described in Prob. 1.9, sup- that the heat absorbed by the air Q is related to the mass of pose that the outflow is not constant but rather depends on the air m the heat capacity, and the change in temperature by the depth. For this case, the differential equation for depth the following relationship: can be written as T 2 1.5 dy Q α(1 + y) Q = m C dT = mC (T − T ) v v 2 1 2 = 3 sin (t) − T 1 dt A A The mass of air can be obtained from the ideal gas law: Use Euler’s method to solve for the depth y from t = 0 to m PV = RT 10 d with a step size of 0.5 d. The parameter values are A = Mwt 2 3 1250 m , Q = 450 m /d, and α = 150. Assume that the ini- tial condition is y = 0. where P is the gas pressure, V is the volume of the gas, Mwt 1.11 Apply the conservation of volume (see Prob. 1.9) to sim- is the molecular weight of the gas (for air, 28.97 kg/kmol), 3 ulate the level of liquid in a conical storage tank (Fig. P1.11). and R is the ideal gas constant 8.314 kPa m /(kmol K).cha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 22 22 MATHEMATICAL MODELING, NUMERICAL METHODS, AND PROBLEM SOLVING (c) Use calculus to obtain the closed form solution where Skin v = v at x = 0. 0 (d) Use Euler’s method to obtain a numerical solution from Urine Feces x = 0 to 100,000 m using a step of 10,000 m where the initial velocity is 1500 m/s upward. Compare your result Food Air with the analytical solution. BODY 1.15 Suppose that a spherical droplet of liquid evaporates at a rate that is proportional to its surface area. Drink Sweat dV =−kA dt Metabolism 3 where V = volume (mm ), t = time (min), k = the evapora- 2 tion rate (mm/min), and A = surface area (mm ). Use FIGURE P1.13 Euler’s method to compute the volume of the droplet from t = 0 to 10 min using a step size of 0.25 min. Assume that k = 0.08 mm/min and that the droplet initially has a radius of 2.5 mm. Assess the validity of your results by determining 1.13 Figure P1.13 depicts the various ways in which an aver- the radius of your final computed volume and verifying that age man gains and loses water in one day. One liter is ingested it is consistent with the evaporation rate. as food, and the body metabolically produces 0.3 liters. In 1.16 A fluid is pumped into the network shown in Fig. P1.16. breathing air, the exchange is 0.05 liters while inhaling, and 3 If Q = 0.7, Q = 0.5, Q = 0.1, and Q = 0.3 m /s, determine 0.4 liters while exhaling over a one-day period. The body will 2 3 7 8 the other flows. also lose 0.3, 1.4, 0.2, and 0.35 liters through sweat, urine, 1.17 Newton’s law of cooling says that the temperature of a feces, and through the skin, respectively. To maintain steady body changes at a rate proportional to the difference between state, how much water must be drunk per day? its temperature and that of the surrounding medium (the am- 1.14 In our example of the free-falling bungee jumper, we bient temperature), assumed that the acceleration due to gravity was a constant 2 value of 9.81 m/s . Although this is a decent approximation dT =−k(T − T ) when we are examining falling objects near the surface of a dt the earth, the gravitational force decreases as we move where T = the temperature of the body (°C), t = time (min), above sea level. A more general representation based on k = the proportionality constant (per minute), and T = the Newton’s inverse square law of gravitational attraction can a ambient temperature (°C). Suppose that a cup of coffee orig- be written as inally has a temperature of 70°C. Use Euler’s method to 2 compute the temperature from t = 0 to 20 min using a step R g(x) = g(0) size of 2 min if T = 20 °C and k = 0.019/min. 2 a (R + x) where g(x) = gravitational acceleration at altitude x (in m) 2 measured upward from the earth’s surface (m/s ), g(0) = 2 ∼ gravitational acceleration at the earth’s surface (= 9.81 m/s ), Q Q Q 1 3 5 6 ∼ and R = the earth’s radius (= 6.37 × 10 m). (a) In a fashion similar to the derivation of Eq. (1.8), use a force balance to derive a differential equation for veloc- Q Q Q Q 2 4 6 7 ity as a function of time that utilizes this more complete representation of gravitation. However, for this deriva- tion, assume that upward velocity is positive. (b) For the case where drag is negligible, use the chain rule Q Q Q 10 9 8 to express the differential equation as a function of alti- tude rather than time. Recall that the chain rule is dv dv dx FIGURE P1.16 = dt dx dtcha01102_ch01_001-023.qxd 12/17/10 7:58 AM Page 23 PROBLEMS 23 using the same parameters and conditions as in Example 1.2. 1.18 You are working as a crime scene investigator and Develop a plot of your results. must predict the temperature of a homicide victim over a 1.20 In addition to the downward force of gravity (weight) 5-hour period. You know that the room where the victim was and drag, an object falling through a fluid is also subject to a found was at 10 °C when the body was discovered. buoyancy force which is proportional to the displaced vol- (a) Use Newton’s law of cooling (Prob. 1.17) and Euler’s ume. For example, for a sphere with diameter d (m), the method to compute the victim’s body temperature for 3 sphere’s volume is V  π d /6, and its projected area is the 5-hr period using values of k  0.12/hr and t  2 A  πd /4. The buoyancy force can then be computed as 0.5 hr. Assume that the victim’s body temperature at F Vg. We neglected buoyancy in our derivation of the time of death was 37 °C, and that the room tempera- b Eq. (1.8) because it is relatively small for an object like a ture was at a constant value of 10°C over the 5-hr bungee jumper moving through air. However, for a more period. dense fluid like water, it becomes more prominent. (b) Further investigation reveals that the room temperature (a) Derive a differential equation in the same fashion as had actually dropped linearly from 20 to 10 °C over the Eq. (1.8), but include the buoyancy force and represent 5-hr period. Repeat the same calculation as in (a) but in- the drag force as described in Sec. 1.4. corporate this new information. (b) Rewrite the differential equation from (a) for the special (c) Compare the results from (a) and (b) by plotting them case of a sphere. on the same graph. (c) Use the equation developed in (b) to compute the terminal 1.19 The velocity is equal to the rate of change of distance, velocity (i.e., for the steady-state case). Use the following x (m): parameter values for a sphere falling through water: dx 3 sphere diameter  1 cm, sphere density  2,700 kg/m , = v(t) (P1.19) dt 3 water density  1,000 kg/m , and C  0.47. d (d) Use Euler’s method with a step size of t  0.03125 s to Use Euler’s method to numerically integrate Eq. (P1.19) and numerically solve for the velocity from t  0 to 0.25 s (1.8) in order to determine both the velocity and distance with an initial velocity of zero. fallen as a function of time for the first 10 seconds of freefall

Advise: Why You Wasting Money in Costly SEO Tools, Use World's Best Free SEO Tool Ubersuggest.