Matlab intro

Solving Dynamics Problems in MATLAB Brian D. Harper Mechanical Engineering The Ohio State University A supplement to accompany Engineering Mechanics: Dynamics, 6th Edition by J.L. Meriam and L.G. Kraige JOHN WILEY & SONS, INC. New York • Chichester • Brisbane • Toronto • Singapore CONTENTS Introduction Chapter 1 5 An Introduction to MATLAB 7 Numerical Calculations Writing Scripts (m-files) Defining Functions Graphics Symbolic Calculations Differentiation and Integration 7 10 12 13 21 24 Solving Equations 26 Chapter 2 37 Kinematics of Particles 2.1 Sample Problem 2/4 (Rectilinear Motion) 2.2 Problem 2/87 (Rectangular Coordinates) 2.3 Problem 2/126 (n-t Coordinates) 2.4 Sample Problem 2/9 (Polar Coordinates) 2.5 Sample Problem 2/10 (Polar Coordinates) 2.6 Problem 2/183 (Space Curvilinear Motion) 2.7 Sample Problem 2/16 (Constrained Motion of Connected Particles) 38 41 46 48 52 55 Chapter 3 61 Kinetics of Particles 3.1 Sample Problem 3/3 (Rectilinear Motion) 3.2 Problem 3/98 (Curvilinear Motion) 3.3 Sample Problem 3/17 (Potential Energy) 3.4 Problem 3/218 (Linear Impulse/Momentum) 3.5 Problem 3/250 (Angular Impulse/Momentum) 3.6 Problem 3/365 (Curvilinear Motion) 57 62 65 67 70 72 73 Chapter 4 Kinetics of Systems of Particles 77 4.1 Problem 4/26 (Conservation of Momentum) 4.2 Problem 4/62 (Steady Mass Flow) 4.3 Problem 4/86 (Variable Mass) 78 80 83 Chapter 5 87 Plane Kinematics of Rigid Bodies 5.1 Problem 5/3 (Rotation) 5.2 Problem 5/44 (Absolute Motion) 5.3 Sample Problem 5/9 (Relative Velocity) 5.4 Problem 5/123 (Relative Acceleration) 5.5 Sample Problem 5/15 (Absolute Motion) 88 93 95 99 100 Chapter 6 107 Plane Kinetics of Rigid Bodies 6.1 Sample Problem 6/2 (Translation) 6.2 Sample Problem 6/4 (Fixed-Axis Rotation) 6.3 Problem 6/98 (General Plane Motion) 6.4 Problem 6/104 (General Plane Motion) 6.5 Sample Problem 6/10 (Work and Energy) 6.6 Problem 6/206 (Impulse/Momentum) Chapter 7 Introduction to Three-Dimensional Dynamics of Rigid Bodies 108 113 115 118 120 125 129 7.1 Sample Problem 7/3 (General Motion) 7.2 Sample Problem 7/6 (Kinetic Energy) 130 132 Chapter 8 137 Vibration and Time Response 8.1 Sample Problem 8/2 (Free Vibration of Particles) 8.2 Problem 8/139 (Damped Free Vibrations) 8.3 Sample Problem 8/6 (Forced Vibration of Particles) 138 140 143 INTRODUCTION Computers and software have had a tremendous impact upon engineering education over the past several years and most engineering schools now incorporate computational software such as MATLAB in their curriculum. Since you have this supplement the chances are pretty good that you are already aware of this and will have to learn to use MATLAB as part of a Dynamics course. The purpose of this supplement is to help you do just that. There seems to be some disagreement among engineering educators regarding how computers should be used in an engineering course such as Dynamics. I will use this as an opportunity to give my own philosophy along with a little advice. In trying to master the fundamentals of Dynamics there is no substitute for hard work. The old fashioned taking of pencil to paper, drawing free body and mass acceleration diagrams, struggling with equations of motion and kinematic relations, etc. is still essential to grasping the fundamentals of Dynamics. A sophisticated computational program is not going to help you to understand the fundamentals. For this reason, my advice is to use the computer only when required to do so. Most of your homework can and should be done without a computer. The problems in this booklet are based upon problems taken from your text. The problems are slightly modified since most of the problems in your book do not require a computer for the reasons discussed in the last paragraph. One of the most important uses of the computer in studying Mechanics is the convenience and relative simplicity of conducting parametric studies. A parametric study seeks to understand the effect of one or more variables (parameters) upon a general solution. This is in contrast to a typical homework problem where you generally want to find one solution to a problem under some specified conditions. For example, in a typical homework problem you might be asked something about the trajectory of a particle launched at an angle of 30 degrees from the horizontal with an initial speed of 30 ft/sec. In a parametric study of the same problem you might typically find the trajectory as a function of two parameters, the launch angle θ and initial speed v. You might then be asked to plot the trajectory for different launch angles and speeds. A plot of this type is very beneficial in visualizing the general solution to a problem over a broad range of variables as opposed to a single case. 6 INTRODUCTION As you will see, it is not uncommon to find Mechanics problems that yield equations that cannot be solved exactly. These problems require a numerical approach that is greatly simplified by computational software such as MATLAB. Although numerical solutions are extremely easy to obtain in MATLAB this is still the method of last resort. Chapter 1 will illustrate several methods for obtaining symbolic (exact) solutions to problems. These methods should always be tried first. Only when these fail should you generate a numerical approximation. Many students encounter some difficulties the first time they try to use a computer as an aid to solving a problem. In many cases they are expecting that they have to do something fundamentally different. It is very important to understand that there is no fundamental difference in the way that you would formulate computer problems as opposed to a regular homework problem. Each problem in this booklet has a problem formulation section prior to the solution. As you work through the problems be sure to note that there is nothing peculiar about the way the problems are formulated. You will see free-body and mass acceleration diagrams, kinematic equations etc. just like you would normally write. The main difference is that most of the problems will be parametric studies as discussed above. In a parametric study you will have at least one and possibly more parameters or variables that are left undefined during the formulation. For example, you might have a general angle θ as opposed to a specific angle of 20°. If it helps, you can “pretend” that the variable is some specific number while you are formulating a problem. This supplement has eight chapters. The first chapter contains a brief introduction to MATLAB. If you already have some familiarity with MATLAB you can skip this chapter. Although the first chapter is relatively brief it does introduce all the methods that will be used later in the book and assumes no prior knowledge of MATLAB. Chapters 2 through 8 contain computer problems taken from chapters 2 through 8 of your textbook. Thus, if you would like to see some computer problems involving the kinetics of particles you can look at the problems in chapter 3 of this supplement. Each chapter will have a short introduction that summarizes the types of problems and computational methods used. This would be the ideal place to look if you are interested in finding examples of how to use specific functions, operations etc. This supplement uses the student edition of MATLAB version 7.1. MATLAB is a registered trademark of The Mathworks, Inc., 24 Prime Park Way, Natick, Massachusetts, 01760. AN INTRODUCTION TO MATLAB This chapter provides an introduction to the MATLAB programming language. Although the chapter is introductory in nature it will cover everything needed to solve the computer problems in this booklet. 1.1 Numerical Calculations When you open MATLAB you should see a prompt something like the following. EDU» To the right of this prompt you will write some sort of command. The following example assigns a value of 200 to the variable a. Simply type “a = 200” at the prompt and then press enter. EDU» a=200 a= 200 There are many situations where you may not want MATLAB to print the result of a command. To suppress the output simply type a semi-colon “;” after the command. For example, to suppress the output in the above case type “a = 200;” and the press enter. The basic mathematical operations of addition, subtraction, multiplication, division and raising to a power are accomplished exactly as you would expect if you are familiar with other programming languages, i.e. with the keys +, -, *, \, and ^. Here are a couple of examples. EDU» c=2*5-120/4 c= -20 1 8 CH. 1 AN INTRODUCTION TO MATLAB EDU» d=2^4 d= 16 MATLAB has many built in functions that you can use in calculations. If you already know the name of the function you can simply type it in. If not, you can find a list of functions in the “Help Desk (HTML)” file available under Help in the main menu bar. Here are a few examples. EDU» e=sin(.5)+tan(.2) e= 0.6821 As with most mathematical software packages, the default unit for angles will be radians. EDU» f=sqrt(16) f= 4 EDU» g=log(5) g= 1.6094 EDU» h=log10(5) h= 0.6990 In the last two examples be careful to note that log is the natural logarithm. For a base 10 logarithm (commonly used in engineering) you need to use log10. Range Variables In the above examples we have seen many cases where a variable (or name) has been assigned a numerical value. There are many instances where we would like a single variable to take on a range of values rather than just a single value. Variables of this type are often referred to as range variables. For example, to have the variable x take on values between 0 and 3 with an increment of 0.25 we would type “x=0:0.25:3”. EDU» x=0:0.25:3 x= Columns 1 through 7 0 0.2500 0.5000 0.7500 1.0000 Columns 8 through 13 1.2500 1.5000 INTRODUCTION TO MATLAB 9 1.7500 2.0000 2.2500 2.5000 2.7500 3.0000 Range variables are very convenient in that they allow us to evaluate an expression over a range with a single command. This is essentially equivalent to performing a loop without actually having to set up a loop structure. Here’s an example. EDU» x=0:0.25:3; EDU» f=3*sin(x)-2*cos(x) f= Columns 1 through 7 -2.0000 -1.1956 -0.3169 0.5815 1.4438 2.2163 Columns 8 through 13 3.3084 3.5602 3.5906 3.3977 2.9936 2.4033 2.8510 Although range variables are very convenient and almost indispensable when making plots, they can also lead to considerable confusion for those new to MATLAB. To illustrate, suppose we wanted to compute a value y = 2*x-3*x^2 over a range of values of x. EDU» x=0:0.5:3; EDU» y = 2*x-3*x^2 ??? Error using ==> mpower Matrix must be square. What went wrong? We have referred to x as a range variable because this term is most descriptive of how we will actually use variables of this type in this manual. As the name implies, MATLAB is a program for doing matrix calculations. Thus, internally, x is really a row matrix and operations such as multiplication or division are taken, by default, to be matrix operations. Thus, MATLAB takes the operation x^2 to be matrix multiplication x*x. Since matrix multiplication of two row matrices doesn’t make sense, we get an error message. In a certain sense, getting this error message is very fortunate since we are not actually interested in a matrix calculation here. What we want is to calculate a certain value y for each x in the range specified. This type of operation is referred to as term-by-term. For term-by-term operations we need to place a period “.” in front of the operator. Thus, for term by term multiplication, division, and raising to a power we would write “.*”, “./” and “.^” respectively. Since matrix and term-by-term addition or subtraction are equivalent, you do not need a period before + or -. The reason that our earlier calculations did not produce an error is that they all involved scalar quantities. To correct the error above we write. EDU» x=0:0.5:3; EDU» y = 2*x-3*x.^2 y= 10 CH. 1 AN INTRODUCTION TO MATLAB 0 0.2500 -1.0000 -3.7500 -8.0000 -13.7500 -21.0000 Note that we wrote “2*x” instead of “2.*x” since 2 is a scalar. If in doubt, it never hurts to include the period. 1.2 Writing Scripts (m-files) If you have been working through the examples above you have probably noticed by now that you cannot modify a line once it has been entered into the worksheet. This makes it rather difficult to debug a program or run it again with a different set of parameters. For this reason, most of your work in MATLAB will involve writing short scripts or m-files. The term m-file derives from the fact that the programs are saved with a “.m” extension. An m-file is essentially just a list of commands that you would normally enter at a prompt in the worksheet. After saving the file you just enter the name of the file at the prompt after which MATLAB will execute all the lines in the file as if they had been entered. The chief advantage is that you can alter the script and run it as many times as you like. To get started, select “File…New…M-File” from the main menu. A screen like that below will then pop up. INTRODUCTION TO MATLAB 11 Now enter the program line by line into the editor. Basically, you just enter the commands you would normally enter at a prompt (do not, however, enter >> at the beginning of a line). If you want to enter comments, begin the line with “%”. These lines will not be executed. Once you have finished entering the script into the editor you can save the file (be sure to remember where you saved it). The file name can contain numbers but cannot begin with a number. It also cannot contain any mathematical operators such as +, -, *, / etc. Also, the name of the file should not be the same as a variable it computes. To run the script all you have to do is type the name of the program at a prompt (>>), but without the “.m” extension. A common error at this point is that MATLAB will give some response indicating that it hasn’t a clue what you are doing. The reason this happens is that the file you saved is not in the current path. You can check (or change) the path by selecting File…Set Path… from the main menu. Following is a simple example of a script file that sets up a range variable x and then calculates two functions over the range. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This script calculates two % functions over a specified range x=0:0.2:1; f=sin(x)./(1+x) g=sin(x).*exp(x) %%%%%%%%% end of script %%%%%%%%%%%%%%%%%%%% The program was saved as example1, so here is the output you should see when the program is entered at a prompt. EDU» example1 f= 0 0.1656 0.2782 0.3529 0.3985 g= 0 0.2427 0.5809 1.0288 1.5965 0.4207 2.2874 To save space we will, in the future, show the output of a script immediately after the script is given without actually showing the file name being entered at a prompt. 12 CH. 1 AN INTRODUCTION TO MATLAB 1.3 Defining Functions Although MATLAB has many built in functions (see the “Help Desk (HTML)” file available under Help in the main menu bar), it is sometimes advantageous to define your own functions. While user defined functions are actually special types of m-files (they will be saved with the .m extension) their use is really quite different from the m-files discussed above. For example, suppose you were to create a function that you call f1. You would then use this function as you would any built in function like sin, cos, log etc. Since functions are really special types of m-files you start by opening an editor window with “File…New…M-File”. The structure of the file is, however, quite different from a script file. Following is an m-file for a function called load1. Note that you should save the file under the same name as the function. The file would thus show up as load1.m in whatever directory you save it. As with script files, you need to be sure the current path contains the directory in which the file was stored. %%%%%% a function m-file %%%%%%%%%%%%%%%% function y = load1(x) % this file defines the function load1(x) y = 2*(x-1)+3*(x-2).^2; %%%%%% end of function m-file %%%%%%%%%%% It is important to understand the structure of a function definition before you try to create some on your own. The name of this function is “load1” and will be used as you would any other function. y and x are local variables which MATLAB uses to calculate the function. There is nothing special about y and x, other variable names will work as well. Now, suppose you were to type load1(3) at a prompt or within a script file. MATLAB will calculate the expression for y substituting x = 3. This value of y will then be returned as the value for load1(3). Here are a few examples using our function load1. EDU» load1(3) ans = 7 EDU» f=load1(5)-load1(3) f= 28 INTRODUCTION TO MATLAB 13 EDU» load1(load1(3)) ans = 87 EDU» x=1:0.2:2; EDU» g=load1(x) g= 3.0000 2.3200 1.8800 1.6800 1.7200 2.0000 1.4 Graphics One of the most useful things about a computational software package such as MATLAB is the ability to easily create graphs of functions. As we will see, these graphs allow one to gain a lot of insight into a problem by observing how a solution changes as some parameter (the magnitude of a load, an angle, a dimension etc.) is varied. This is so important that practically every problem in this supplement will contain at least one plot. By the time you have finished reading this supplement you should be very proficient at plotting in MATLAB. This section will introduce you to the basics of plotting in MATLAB. MATLAB has the capability of creating a number of different types of graphs. Here we will consider only the X-Y plot. The most common and easiest way to generate a plot of a function is to use range variables. The following example will guide you through the basic procedure. EDU» x=-3:0.1:3; EDU» f=x.*exp(-x.^2); EDU» plot(x,f) The procedure is very simple. First define a range variable covering the range of the plot, then define the function to be plotted and issue the plot command. After typing in the above you should see a graph window pop up which looks something like the following figure. 14 CH. 1 AN INTRODUCTION TO MATLAB Things like titles and axis labels can be added in one of two ways. The first is by line commands such as those shown below. EDU» xlabel('x') EDU» ylabel('x*exp(-x^2)') EDU» title('A Simple Plot') Take a look at the graph window after you type each of the lines above. At this point, the graph window should look something like that shown below. INTRODUCTION TO MATLAB 15 There are many other things that you can add or change on the plot. The easiest way to do this is with the second of the two approaches mentioned above. In the graph menu select “Tools…Edit Plot”. Now you can edit the plot in a number of different ways. To see some of the options try either double clicking or right clicking on the graph. You should spend some time experimenting with a graph window to see the range of possibilities. In particular, be sure to try clicking the show plot tools button . Once you have things like you want them you can save the graph to a file. You can also export the graph to some format suitable for inserting in, say, a word processor. Select “File…Export Setup” in the graph window to see the range of possible formats. 16 CH. 1 AN INTRODUCTION TO MATLAB Parametric Studies One of the most important uses of the computer in studying Mechanics is the convenience and relative simplicity of conducting parametric studies (not to be confused with parametric plotting discussed below). A parametric study seeks to understand the effect of one or more variables (parameters) upon a general solution. This is in contrast to a typical homework problem where you generally want to find one solution to a problem under some specified conditions. For example, in a typical homework problem you might be asked to find the reactions at the supports of a structure with a concentrated force of magnitude 200 lb that is oriented at an angle of 30 degrees from the horizontal. In a parametric study of the same problem you might typically find the reactions as a function of two parameters, the magnitude of the force and its orientation. You might then be asked to plot the reactions as a function of the magnitude of the force for several different orientations. A plot of this type is very beneficial in visualizing the general solution to a problem over a broad range of variables as opposed to a single case. Parametric studies generally require making multiple plots of the same function with different values of a particular parameter in the function. As a simple example, consider the following function. f = 5 + x − 5x 2 + ax3 What we would like to do is gain some understanding of how f varies with both x and a. It might be tempting to make a three dimensional plot in a case like this. Such a plot can, in some cases, be very useful. Usually, however, it is too difficult to interpret. This is illustrated by the following three dimensional plot of f versus x and a. % This script produces a 3-d plot of % f versus x and a. x=linspace(-10,10,30); a=linspace(-3,3,30); [X,A]=meshgrid(x,a); F=5+X-5*X.^2+A.*X.^3; B=0.*A+2; mesh(X,A,F) colormap gray %%%%%%% end of script %%%%%%%%%%%%%%%%%%%% INTRODUCTION TO MATLAB 17 The plot above certainly is interesting but, as mentioned above, not very easy to interpret. In most cases it is much better to plot the function several times (with different values of the parameter of interest) on a single two-dimensional graph. We will illustrate this by plotting f as a function of x for a = -1, 0, and 1. This is accomplished in the following script. % script for plotting f versus x for several % values of a. x=-5:0.1:5; g=5+x-5*x.^2; f1=g-x.^3; % a = -1 f2=g; % a = 0 f3=g+x.^3; % a = 1 plot(x,f1,x,f2,x,f3) xlabel('x') ylabel('f') %%%%%%%% end of script %%%%%%%%%%%%%%%%%%%%%% 18 CH. 1 AN INTRODUCTION TO MATLAB Running this script results in the following plot. Parametric Plots It often happens that one needs to plot some function y versus x but y is not known explicitly as a function of x. For example, suppose you know the x and y coordinates of a particle as a function of time but want to plot the trajectory of the particle, i.e. you want to plot the y coordinate of the particle versus the x coordinate. A plot of this type is generally called a parametric plot. Parametric plots are easy to obtain in MATLAB. You start by defining the two functions in terms of the common parameter and then define the common parameter as a range variable. Then issue a plot command using the two functions as arguments. The following script illustrates this procedure. INTRODUCTION TO MATLAB 19 % This script file illustrates parametric plotting. % A function f is plotted versus another function g. % The two functions are related by the % common parameter a. a=-1:0.05:3.5; f=10*a.*(2-a); g=sin(3*a); plot(g,f) xlabel('g') ylabel('f') %%%%%%%% end of script %%%%%%%%%%%%%%%%%%%%%%%%%%% The range selected for the parameter can have a big, and sometimes surprising effect on the resulting graph. To illustrate, try increasing the upper limit on the range on a a few times and see how the graph changes. You can, of course, also plot g as a function of f. 20 CH. 1 AN INTRODUCTION TO MATLAB % Another example of parametric plotting a=-1:0.05:6; f=10*a.*(2-a); g=sin(3*a); plot(f,g) xlabel('f') ylabel('g') %%%%%%%%% end of script %%%%%%%%%%%%%%%%% INTRODUCTION TO MATLAB 21 1.5 Symbolic Calculations One of the most useful features of many modern mathematical software packages is the capability of doing symbolic math. By symbolic math we mean mathematical operations on expressions containing undefined variables rather than numbers. Maple is, perhaps, the best known computer algebra program. Many other programs (including MATLAB) have a Maple engine that is called whenever symbolic operations are performed. The Maple engine is part of the Symbolic Math Toolbox so, if you do not have this toolbox, you will not be able to do symbolic math. If you have a student edition of MATLAB you should have the Symbolic Math Toolbox. If you do not have the Symbolic Math Toolbox, you will need to perform whatever symbolic calculations are done in this manual by hand. How does MATLAB decide whether to call Maple and attempt a symbolic calculation? Perhaps the best way to learn this is by studying the examples below and by practicing. Maple will always be used to manipulate any expression represented as a character string (i.e. an expression enclosed in ‘quotation’ marks). In the remainder of this manual we will refer to character strings of this type as Maple expressions. Here are a few examples of Maple expressions. EDU» 'x+a*x^2-b*x^3' ans = x+a*x^2-b*x^3 EDU» f='A*sin(x)*exp(-x^2)' f= A*sin(x)*exp(-x^2) Instead of placing the entire expression in quotations, you can also first declare certain variables to be symbols and then write the expressions as usual. This is more convenient in some situations but can get confusing. Here are some examples of this type. EDU» syms x a b A % Sets x, a, b and A as symbols EDU» g=x+a*x^2-b*x^3 g= x+a*x^2-b*x^3 EDU» f=A*sin(x)*exp(-x^2) f= A*sin(x)*exp(-x^2) The above are just examples of Maple expressions and do not involve any symbolic math. Two of the most important applications of symbolic math will be discussed in the next two sections, namely symbolic calculus (integration and 22 CH. 1 AN INTRODUCTION TO MATLAB differentiation) and symbolic solution of one or more equations. The purpose of the present section is to introduce you to the basic procedures of symbolic math as well as to give a few other useful applications. One particularly useful symbolic operation is substitution and is illustrated by the following examples. EDU» g='x+a*x^2-b*x^3' g= x+a*x^2-b*x^3 EDU» subs(g,'a',3) ans = x+3*x^2-b*x^3 % substitutes 3 for a in g. Note the quotation marks about a. EDU» g1=subs(g,'x','(y-c)') g1 = ((y-c))+a*((y-c))^2-b*((y-c))^3 EDU» syms x y A B EDU» f=A*sin(x)*exp(-B*x^2) f= A*sin(x)*exp(-B*x^2) EDU» subs(f,B,4) % Note that B doesn't have to be in quotes since it was declared to be a symbol. ans = A*sin(x)*exp(-4*x^2) Another very convenient feature available in the Symbolic Toolbox provides the capability of defining “inline” functions thus avoiding the necessity for creating special m-files as described above. Here is an example. EDU» f=inline('x^2*exp(-x)') f= Inline function: f(x) = x^2*exp(-x) EDU» f(3) ans = 0.4481 In the previous section we plotted the function f = 5 + x − 5x 2 + ax3 for several values of a. It turns out to be somewhat more convenient to do this with an inline function. Let’s see how it works. INTRODUCTION TO MATLAB 23 EDU» f=inline('5+x-5*x^2+a*x^3') f= Inline function: f(a,x) = 5+x-5*x^2+a*x^3 At this point we might anticipate a minor problem. When we plot this as a function of x we will have a range variable as one of the arguments for f. This creates a problem as MATLAB will be performing the numerical calculations and it will take the operations in f as matrix operations rather than term by term. Fortunately we can use the “vectorize” operator to insert periods in the appropriate places. EDU» f=vectorize(f) f= Inline function: f(a,x) = 5+x-5.*x.^2+a.*x.^3 EDU» x=-5:0.1:5; EDU» plot(x,f(-1,x),x,f(0,x),x,f(1,x)) 24 CH. 1 AN INTRODUCTION TO MATLAB 1.6 Differentiation and Integration Here we will consider only symbolic differentiation and integration since these are most useful in the mechanics applications in this manual. For differentiation we will use the Maple diff command. Care should be exercised here since diff is also a MATLAB command. To guarantee that you are using the Maple command, be sure that you use diff only on Maple expressions. EDU» f='a*sec(b*t)'; EDU» diff(f, 't') % Note that t must be placed in quotation marks since it is a symbol. ans = a*sec(b*t)*tan(b*t)*b We can avoid using the quotation marks by making a, b and t symbols. EDU» syms a b t EDU» f=a*sec(b*t); EDU» diff(f, t) ans = a*sec(b*t)*tan(b*t)*b To find higher order derivatives use the command diff(f, x, n) where n is the order of the derivative. For example, to find the third derivative of alog(b+x) we would write EDU» diff('a*log(b+x)', 'x', 3) ans = 2*a/(b+x)^3 Symbolic integration will be performed with the int command. The general format of this command is int(f, x, a, b) where f is the integrand (a Maple expression), x is the integration variable and a and b are the limits of integration. If x, a and b have not been declared symbols with the syms command, they must be placed in quotation marks. If the integration limits are omitted, the indefinite integral will be evaluated. Here are several examples of definite and indefinite integrals. EDU» int('sin(b*x)','x') ans = -1/b*cos(b*x) INTRODUCTION TO MATLAB 25 EDU» int('sin(b*x)','x','c','d') ans = -(cos(b*d)-cos(b*c))/b EDU» syms x b c d EDU» f=b*log(x); EDU» int(f,x) ans = b*x*log(x)-b*x EDU» int(f,x,c,d) ans = b*d*log(d)-b*d-b*c*log(c)+b*c If a definite integral contains no unknown parameters either in the integrand or the integration limits, the int command will provide numerical answers. Here are a few examples. EDU» int('x+3*x^3','x',0,3) ans = 261/4 EDU» int('log(x)','x',2,5) % Don’t forget that log is the natural logarithm. ans = 5*log(5)-3-2*log(2) Note that Maple will always try to return an exact answer. This usually results in answers containing fractions or functions as in the above examples. This is very useful in some situations; however, one often wants to know the numerical answer without having to evaluate a result such as the above with a calculator. To obtain numerical answers use Maple’s eval function as in the following examples. EDU» eval(int('x+3*x^3','x',0,3)) ans = 65.2500 EDU» eval(int('log(x)','x',2,5)) ans = 3.6609 26 CH. 1 AN INTRODUCTION TO MATLAB 1.7 Solving Equations Our preferred approach to solving equations is to use Maple’s solve command available in the Symbolic Toolbox. We will also discuss briefly MATLAB’s built in function fzero which can be used to solve single equations numerically. Solving single equations Here are a couple of examples illustrating how to use Maple’s solve command to solve a single equation symbolically. The first example is the (hopefully) familiar quadratic equation. EDU» f='a*x^2+b*x+c=0'; EDU» solve(f, 'x') ans = 1/2/a*(-b+(b^2-4*a*c)^(1/2)) 1/2/a*(-b-(b^2-4*a*c)^(1/2)) EDU» g='a*sin(x)-cos(x)=1'; EDU» solve(g, 'x') ans = pi atan2(2*a/(1+a^2),(a^2-1)/(1+a^2)) Note that both examples produced two solutions with each solution being placed on separate lines. One of the solutions to the second example is the number π. The second solution illustrates something you should always watch out for. In many cases you may get a solution in terms of a function that you are not familiar with, i.e. atan2 that appears in the second solution to the second example above. Whenever this happens you should type “help atan2” at a prompt in order to learn more about the function. If the variable being solved for is the only unknown in the equation, solve will return a number as the result. Here are a couple of examples. EDU» f='2*x^2+4*x-12'; EDU» solve(f, 'x') ans = -1+7^(1/2) -1-7^(1/2) EDU» g='5*sin(x)-cos(x)=1'; EDU» solve(g, 'x') ans = INTRODUCTION TO MATLAB 27 pi atan(5/12) If you want the answer as a floating point number instead of an exact result, use the Maple command eval. For example, with g defined as above we would write EDU» eval(solve(g, 'x')) ans = 3.1416 0.3948 Although Maple will always attempt to return an exact (symbolic) result, you may occasionally have a case where a floating point number is returned without using eval. Here is an example. EDU» f='10*sin(z)=2*z+1'; EDU» solve(f, 'z') ans = .12541060097666150276238831440340 Whenever this occurs, Maple was unable to solve the equation exactly and switched automatically to a numerical solution. In a certain sense this is very convenient as the same command is being used for both symbolic and numerical calculations. The problem is that nonlinear equations quite often have more than one solution, as illustrated by the examples above. If you have the full version of Maple then you will have various options available for finding multiple numerical solutions. But since we are dealing now with numerical rather than symbolic solutions it makes more sense just to use MATLAB. The appropriate MATLAB command for the numerical solution of an equation is fzero. This command numerically determines the point at which some defined function has a value of zero. This means that we need to rewrite our equation in terms of an expression whose zeros (roots) are solutions to the original expression. This is easily accomplished by rearranging the original equation into the form expression = 0. We will call this expression g in this example. Rewriting the above equation f we have g = 2 z + 1 − 10 sin( z ) . The general format for fzero is fzero(‘function_name’, x0) where ‘function_name’ is the name of the m-file function that you have created and x0 is an initial guess at the root (zero) of the function. The initial guess is very important since, if there is more than one solution, MATLAB will find the solution closest to the initial guess. Well, you may be wondering how we can determine an approximate solution if we do not yet know the solution. Actually, this is very easy to do. First we define a function g(z) whose roots will be the solution to the equation of interest. Next, we plot this function in order to estimate the location of points where g(z) = 0. Here’s how it works for the present example. 28 CH. 1 AN INTRODUCTION TO MATLAB %%%%% function m-file %%%%%% function y = gz(z) y = 2*z + 1 - 10*sin(z); %%%%% end of m-file %%%%%%%% Above is our user defined function (m-file) which we have named gz. Now we can plot this function over some range as illustrated below. EDU» z=-1:0.01:3; EDU» plot(z,gz(z),z,0) % we also plot a line at g=0 to help find the roots. EDU» xlabel('z') EDU» ylabel('g(z)') The plot shows that the function is zero at about z = 0 and 2. The first root is clearly that found by Maple’s solve command above. Now we can easily find both roots with the fzero command. EDU» fzero('gz', 0) ans = 0.1254 EDU» fzero('gz', 2) ans = 2.4985 INTRODUCTION TO MATLAB 29 Thus, the two solutions are z = 0.1254 and 2.4985. Finding Maxima and Minima of Functions Here we will illustrate two methods for finding maxima and minima of a function. Method 1. Using the symbolic diff and solve functions. The usual method for finding maxima or minima of a function f(x) is to first determine the location(s) x at which maxima or minima occur by solving the df equation dx = 0 for x. One then substitutes the value(s) of x thus determined into f(x) to find the maximum or minimum. Consider finding the maximum and minimum values of the following function: f = 1 + 2x − x3 Before proceeding, it is a good idea to make a plot of the function f. This will give us a rough idea where the minima and maxima are as a check on the results we obtain below. The following script will plot f as a function of x. % This script plots f as a function of x x = -2:0.05:2; f = 1 + 2*x - x.^3; plot(x, f) xlabel('x') ylabel('f') %%%%%%%%%% end of script %%%%%%%%%%%%% 30 CH. 1 AN INTRODUCTION TO MATLAB This plot shows we have a minimum of about 0 at x of about −0.9 and a maximum of about 2 at x of about 0.9. Now we can proceed to find the results more precisely. EDU» f = '1+2*x-x^3' f= 1+2*x-x^3 EDU» dfdx = diff(f,'x') dfdx = 2-3*x^2 EDU» solve(dfdx,'x') ans = 1/3*6^(1/2) -1/3*6^(1/2) If an equals sign is omitted (as above) solve will find the root of the supplied Maple expression. In other words, it will solve the equation dfdx = 0. It is more convenient in the present case to have the solutions as floating point numbers. Thus, we should have used eval as in the following. INTRODUCTION TO MATLAB 31 EDU» eval(solve(dfdx,'x')) ans = 0.8165 -0.8165 The above results are the locations (x) where the minima and maxima occur. To find the maximum and minimum values of the function f we use subs to substitute these results back into f. EDU» subs(f,'x',0.8165) ans = 2.0887 EDU» subs(f,'x',-0.8165) ans = -0.0887 Thus, fmax = 2.0887 at x = 0.8165 and fmin = − 0.0887 at x = − 0.8165. Method 2. Using MATLAB’s min and max functions. EDU» x = -2:0.05:2; EDU» f = 1 + 2*x - x.^3; These two lines are identical to those in the script used to plot f above. If you have just run that script you do not need to execute these lines. Now we can find the minimum and maximum values of f by typing min(f) and max(f). EDU» min(f) ans = -3 EDU» max(f) ans = 5 These are clearly not the answers we obtained above. What went wrong? It is important to understand that min and max do not find true minima and maxima in the mathematical sense, i.e. they do not find locations where df/dx = 0. Instead, they find the minimum or maximum value of all those calculated in the vector f. To see this, go back and look at the plot of f above and you will see that the minimum and maximum values of f are indeed − 3 and 5. This is another good reason to plot the function first. To avoid the above difficulty we can change the range of x to insure we pick up the true minima and maxima. EDU» x = -1:0.05:1; 32 CH. 1 AN INTRODUCTION TO MATLAB EDU» f = 1 + 2*x - x.^3; In most cases, we want not only a minimum or maximum but also the location where they occur. This being the case we should use the following command. EDU» [f_min,i] = min(f) f_min = -0.0880 i=5 The command [f_min, i] = min(f) finds the minimum value of f and assigns it to the variable f_min. The index of this value is assigned to the variable i. Now we can find the location of the minimum by printing the 5th value of the vector x. EDU» x_min=x(i) x_min = -0.8000 Now we can repeat the above for the maximum. EDU» [f_max,i] = max(f) f_max = 2.0880 i = 37 EDU» x_max=x(i) x_max = 0.8000 You may have noticed at this point that the answers are somewhat different from those found by the first method above. The reason is, once again, that min and max do not find true minima and maxima. To obtain more accurate results you can use a closer spacing for the range variable (x). To see this, repeat the above for a spacing of 0.01 instead of 0.05. INTRODUCTION TO MATLAB 33 Solving several equations simultaneously EDU» eqn1='x^2+y^2=12' eqn1 = x^2+y^2=12 EDU» eqn2='x*y=4' eqn2 = x*y=4 In the above we have defined two equations which we will now solve for the two unknowns, x and y. EDU» [x,y] = solve(eqn1,eqn2) x= 5^(1/2)-1 -1-5^(1/2) 5^(1/2)+1 1-5^(1/2) y= 5^(1/2)+1 1-5^(1/2) 5^(1/2)-1 -1-5^(1/2) Note that four solutions have been found, the first being x = 5 − 1 , y = 5 + 1 . The following example illustrates a case where there are more unknowns than there are equations. In this case, the specified variables will be solved for in terms of the others. EDU» f='x^2 + a*y^2 = 0' f= x^2 + a*y^2 = 0 EDU» g='x-y = b' g= x-y = b 34 CH. 1 AN INTRODUCTION TO MATLAB EDU» [x,y] = solve(f,g) x= 1/2/(a+1)*(-2+2*(-a)^(1/2))*b+b 1/2/(a+1)*(-2-2*(-a)^(1/2))*b+b y= 1/2/(a+1)*(-2+2*(-a)^(1/2))*b 1/2/(a+1)*(-2-2*(-a)^(1/2))*b In this case, there are two solutions. Now let’s consider an example with three equations. EDU» eqn1='x^2+y^2=12' eqn1 = x^2+y^2=12 EDU» eqn2='x*y=4' eqn2 = x*y=4 EDU» eqn3='x-y=z' eqn3 = x-y=z EDU» [x,y,z] = solve(eqn1,eqn2,eqn3) x= 5^(1/2)+1 1-5^(1/2) 5^(1/2)-1 -1-5^(1/2) INTRODUCTION TO MATLAB 35 y= 5^(1/2)-1 -1-5^(1/2) 5^(1/2)+1 1-5^(1/2) z= 2 2 -2 -2 KINEMATICS OF PARTICLES Kinematics involves the study of the motion of bodies irrespective of the forces that may produce that motion. MATLAB can be very useful in solving particle kinematics problems. Problem 2.1 is a rectilinear motion problem illustrating integration with the int command. The formulation of this problem results in an equation that cannot be solved exactly except with some rather sophisticated mathematics. When this occurs it is generally easiest to obtain either a graphical or numerical solution. This problem illustrates both approaches. Problem 2.2 is a rectangular coordinates problem that illustrates diff and solve. Problem 2.3 is a relatively straightforward problem where MATLAB is used to generate a plot that might be useful in a parametric study. The r-θ components of velocity and acceleration are plotted in problem 2.4 as well as the trajectory of a particle. In problem 2.5, the r-θ components of the velocity are determined using symbolic differentiation (diff). The problem also illustrates how computer algebra can simplify what might normally be a rather tedious algebra problem. The diff command is further illustrated in problems 2.6 and 2.7. Problem 2.7 is particularly interesting in that it requires differentiation with respect to time of a function whose explicit time dependence is unknown. This happens rather frequently in Dynamics so it is useful to know how to accomplish this with MATLAB. 2 38 CH. 2 KINEMATICS OF PARTICLES 2.1 Sample Problem 2/4 (Rectilinear Motion) A freighter is moving at a speed of 8 knots when its engines are suddenly stopped. From this time forward, the deceleration of the ship is proportional to the square of its speed, so that a = − kv 2 . The sample problem in your text shows that it is rather easy to determine the constant k by measuring the speed of the boat at some specified time. Show how k could be found by (a) measuring the speed after some specified distance and (b) measuring the time required to travel some specified distance. In both cases let the initial speed be v0. Problem Formulation (a) Since time is not involved, the easiest approach is to integrate the equation vdv = ads. v ∫ vdv = ads = −kv ds 2 v0 s v  ks = ln 0  v dv = − k ds v 0 ∫ With this result it is easy to find k given v at some specified s. To illustrate, assume that v0 = 8 knots and that the speed of the boat is determined to be 3.9 knots after it has traveled one nautical mile.  8  k (1) = ln    3. 9  k = 0.718 mi-1 (b) Here we follow the general approach in the sample problem. Integrating a = dv/dt yields v dv ∫v 2 t ∫ − kt = = − k dt v0 0 v − v0 vv0 v= v0 1 + ktv 0 To obtain the distance s as a function of time we integrate v = ds/dt s ∫ 0 t ∫ ds = s = vdt = 0 t v0 ∫ 1 + ktv 0 dt 0 s= 1 ln(1 + ktv 0 ) k KINEMATICS OF PARTICLES 39 This equation turns out to be very difficult to solve for k. A good mathematician or someone familiar with symbolic algebra software might be able to find the general solution for k in terms of the so-called LambertW function (LambertW(x) is the solution of the equation yey = x). Even if this solution were found it would be of little use in most practical situations. For example, you would have to spend some time familiarizing yourself with the function. Once this is done you would still have to use a program like Maple or a mathematical handbook to evaluate the function. For these reasons it is probably easiest to find k either graphically or numerically. Obtaining a numerical solution with MATLAB is so easy that there is little reason not to use this approach. It is generally advisable though to use a graphical approach even when a numerical solution is being obtained. This is the best way to identify whether there are multiple solutions to the problem and also serves as a useful check on the numerical results. Thus, both approaches are illustrated below. The usual way to generate a graphical solution is to rearrange the equation so as to give a function that is zero at points that are solutions to the original equation. Rearranging the equation above in this manner yields, f = ks − ln(1 + ktv0 ) = 0 Given values of s, t, and v0, f can be plotted versus k. The value of k at which f = 0 provides the solution to the original equation. MATLAB Worksheet and Scripts Although the integrations are simple in this problem, we'll go ahead and evaluate them symbolically for purposes of illustration. EDU» s_a='-1/k'*int('1/x','x','v0','v') % part (a) s_a = -1/k*(log(v)-log(v0)) EDU» s_b=int('v0/(1+k*x*v0)','x',0,'t') % part (b) s_b = log(1+t*k*v0)/k The following script implements a graphical solution for part (b). To illustrate, we take v0 = 8 knots and assume that the boat is found to move 1.1 nautical miles after 10 minutes. 40 CH. 2 KINEMATICS OF PARTICLES %%%%%%%%%%%% Script %%%%%%%%%%%%%%%%%%%% % This script plots f as a function of k for % given values of v0, s, and t v0 = 8; % initial speed (knots) s = 1.1; % distance (nautical miles) t = 10/60; % time (hours) k = 0:0.01:0.5; f = s*k-log(1+k*v0*t); line = k*0; plot(k,f,k,line) xlabel('k') ylabel('f') The above graph shows that k is about 0.34 mi-1. Now let's try to find a solution with the symbolic solve. First we write a Maple expression for f given our assumed values for v0, s, and t. EDU» f = '1.1*k - log(1+k*8*10/60) = 0' f = 1.1*k - log(1+k*8*10/60)=0 EDU» solve(f, 'k') ans = 0. .33923053342470867736031513009887 Thus, k = 0.3392 mi-1. KINEMATICS OF PARTICLES 41 2.2 Problem 2/87 (Rectangular Coordinates) A long-range rifle is fired at A with the projectile hitting the mountain at B. (a) If the muzzle velocity is u ⵜ 400 m/s, determine the two angles of elevation θ which will permit the projectile to hit the mountain target B and plot the two trajectories. (b) Determine the smallest muzzle velocity that will allow the projectile to strike at B and the angle at which it must be fired. Repeat the plot for part (a) and include the trajectory of the projectile for this minimum initial velocity. Problem Formulation Place a coordinate system at A with x positive to the right and y positive up. The initial components of the velocity are, (v x )0 = u cosθ (v y )0 = u sin θ The acceleration is constant with components a x = 0 and a y = − g . Integrating these two accelerations twice and applying initial conditions yields (see page 44 of your text if you need additional details), x = ucosθ t y = usinθ t -1/2gt2 Plotting y in terms of x for different times t will yield the trajectory of the projectile. This type of plot is called a parametric plot since the items plotted (x and y) are each known in terms of another parameter (t). Anytime you have a projectile motion problem and you know the coordinates of a point on the trajectory (our point B) you should solve for x and y (as we have done above) and then obtain two equations by substituting the coordinates of the points. These two equations can then be solved for two unknowns. Note that in most cases one of the two unknowns will be the time of flight. Part (a) Substituting x = 5,000 m, y = 1,500 m and u = 400 m/s gives 5000 = 400cosθ t 1500 = 400sinθ t -1/2(9.81)t2 42 CH. 2 KINEMATICS OF PARTICLES We will let MATLAB solve these two equations simultaneously. MATLAB actually finds four solutions, however two of the four can be discarded since they involve negative times. The other two solutions correspond to the two solutions shown in the illustration accompanying the problem statement. The results are, θ1 = 26.6° and θ2 = 80.6° Part (b) It should be intuitively obvious why there must be a minimum initial velocity below which the projectile cannot reach B. How do we go about finding it? We still have the two equations for the coordinates of point B, 5000 = ucosθ t 1500 = usinθ t -1/2(9.81)t2 however there are now three unknowns (u, θ, t). Suppose for the moment that the launch angle θ were given and we were asked to calculate the required initial speed u so that the projectile strikes B. In this case we would have two equations and two unknowns. From this observation we see that u is a function of θ from which we get our general solution strategy: (a) (b) Eliminate t from the above two equations and solve for u as a function of θ. Differentiate this function with respect to θ to find the location of the minimum. Solving the first equation for u gives u= 5000 t cos θ 1 2 gt . This equation is 2 now solved for t = 2(5000 tan θ − 1500 ) / g which can be substituted back into u to give Substituting into the second yields 1500 = 5000 tan θ − u= 5000 cos θ 2(5000 tan θ − 1500 ) / g We will let MATLAB differentiate this equation and solve for the minimum speed and the associated launch angle. The result is umin = 256.8 m/s at θ = 53.3° KINEMATICS OF PARTICLES 43 MATLAB Worksheet and Scripts Part (a) EDU>> eqn1='5000 = 400*cos(x)*y' eqn1 = 5000 = 400*cos(x)*y EDU>> eqn2='1500 = 400*sin(x)*y-1/2*9.81*y^2' eqn2 = 1500 = 400*sin(x)*y-1/2*9.81*y^2 EDU>> [x,y] = solve(eqn1,eqn2) x= -1.7350349681568043127091034228443 -2.6858972177500184529892570306420 .45569543583977478547338635263748 1.4065576854329889257535399604352 y= -76.452007728583426926316546391465 -13.920516408172829238156414017816 13.920516408172829238156414017816 76.452007728583426926316546391465 MATLAB finds four solutions but two can be excluded because they have time (y) being negative. The two positive times and corresponding angles are and t B1 = 13.9205 secs θ1 = 0.4557 rads (26.6°) t B2 = 76.4520 secs θ2 = 1.4066 rads (80.6°) %=========== script #1 ====================== % Note that we need to set up two time scales in order to % insure that both trajectories end at point B. t1 = 0:0.01:13.92; t2 = 0:0.01:76.45; x1 = 400*cos(.455695)*t1; y1 = 400*sin(.455695)*t1-1/2*9.81*t1.^2; x2 = 400*cos(1.4066)*t2; y2 = 400*sin(1.4066)*t2-1/2*9.81*t2.^2; plot(x1/1000,y1/1000,x2/1000,y2/1000) xlabel('x (km)') ylabel('y (km)') title('plot for part (a)') %======= end script #1 ====================== 44 CH. 2 KINEMATICS OF PARTICLES Part (b) EDU>> u = '5000/cos(x)/sqrt(2*(5000*tan(x)-1500)/9.81)' u = 5000/cos(x)/sqrt(2*(5000*tan(x)-1500)/9.81) EDU>> dudx = diff(u,'x'); EDU>> solve(dudx,'x') ans = .93112656063638185561346315653632 -.63966976615851476361785853510343 Only one of the two solutions is between 0 and 90°. EDU>> .9311*180/pi ans = 53.3481 KINEMATICS OF PARTICLES 45 EDU>> subs(u,.9311) ans = 256.7581 Thus, umin = 256.8 m/s at θ = 53.3° Before plotting we need to find the time to reach point B. This can be done from the equation x = 5000 = ucosθ. EDU>> 5000/256.7581/cos(.9311) ans = 32.6217 %=========== script #2 ====================== % this script should be run after script 1 t3 = 0:0.01:32.62; x3 = 256.8*cos(.9311)*t3; y3 = 256.8*sin(.9311)*t3-1/2*9.81*t3.^2; plot(x1/1000,y1/1000,x2/1000,y2/1000,x3/1000,y3/1000) title('plot for part (b)') %======= end script #2 ====================== 46 CH. 2 KINEMATICS OF PARTICLES 2.3 Problem 2/126 (n-t Coordinates) A baseball player releases a ball with initial conditions shown in the figure. Plot the radius of curvature of the path just after release and at the apex as a function of the release angle θ. Explain the trends in both results as θ approaches 90°. Problem Formulation Just after release a n = g cos θ = v 02 ρ v 02 ρ= g cos θ At the apex At the apex, vy = 0 and v = vx = v0cosθ. Since v is horizontal, the normal direction is vertically downward so that an = g. an = g = (v0 cosθ )2 ρ ρ= (v0 cosθ )2 g MATLAB Script %%%%%%%%%%%%%%%%% Script %%%%%%%%%%%%%%%%%%%% v0 = 100; g = 32.2; theta = 0:0.01:pi/2; rho_i = v0^2/g./cos(theta); rho_a = (v0*cos(theta)).^2/g; plot(theta*180/pi,rho_i, theta*180/pi,rho_a) xlabel('theta (deg)') title('radius of curvature (ft)') axis([0 90 0 800]) % note that we need to limit the vertical axis since the % initial rho approaches infinity as theta approaches % 90 degrees. %%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%% KINEMATICS OF PARTICLES 47 Note that as θ approaches 90°, the initial ρ goes to infinity while ρ at the apex approaches zero. When θ = 90°, the ball travels along a straight (vertical) path. As you recall, straight paths have a radius of curvature of infinity. At the apex, the velocity will be zero giving a radius of curvature of zero. 48 CH. 2 KINEMATICS OF PARTICLES 2.4 Sample Problem 2/9 (Polar Coordinates) Rotation of the radially slotted arm is governed by θ = 0.2t + 0.02t3, where θ is in radians and t is in seconds. Simultaneously, the power screw in the arm engages the slider B and controls its distance from O according to r = 0.2 + 0.04t2, where r is in meters and t is in seconds. Calculate the magnitudes of the velocity and acceleration of the slider as a function of time t. (a) Plot v, vr and vθ for t between 0 and 5 sec. (b) Plot a, ar and aθ for t between 0 and 5 sec. (c) Plot the path of the slider B and compare with the result in your book. Problem Formulation The first part of this problem solution will be identical to that in the Sample Problem in your text except that everything will be left in terms of t. To summarize, r = 0.2 + 0.04t 2 r& = 0.08t &r& = 0.08 θ = 0.2t + 0.02t 3 θ& = 0.2 + 0.06t 2 θ&& = 0.12t Now all we have to do is substitute these expressions into the definitions for the velocity and acceleration. As usual, there is no need to make an explicit substitution when using the computer. v r = r& = 0.08t vθ = rθ& v = v r2 + vθ2 a r = &r& − rθ& 2 aθ = rθ&& + 2r&θ& a = a r2 + aθ2 The plot for part (c) can be obtained by using the suggestion in your book. First write x = r cosθ y = r sin θ and then plot y versus x. This type of plot is called a parametric plot since y is no known explicitly as a function of x. Instead, both x and y are known in terms of another parameter, namely t. KINEMATICS OF PARTICLES 49 MATLAB Scripts %======= script #1 ============================ % this script plots the velocities for part (a) t = 0:0.01:5; r = 0.2+0.04*t.^2; rd = 0.08*t; rdd = 0.08; theta = 0.2*t+0.02*t.^3; thetad = 0.2+0.06*t.^2; thetadd = 0.12*t; vr = rd; vtheta = r.*thetad; v=sqrt(vr.^2+vtheta.^2); plot(t,v,t,vr,t,vtheta) xlabel('time (seconds)') ylabel('velocity (m/s)') title('part (a) velocity') %====== end of script#1======================= 50 CH. 2 KINEMATICS OF PARTICLES %======= script #2 ============================ % this script plots the accelerations for part (b) % run script # 1 first ar = rdd - r.*thetad.^2; atheta = r.*thetadd+2*rd.*thetad; a=sqrt(ar.^2+atheta.^2); plot(t,a,t,ar,t,atheta) xlabel('time (seconds)') ylabel('acceleration (m/s^2)') title('part (b) acceleration') %=========end of script 2====================== %======= script #3 ==================== % This script plots the path of the % slider for part (b). % Run script # 1 first. x = r.*cos(theta); y = r.*sin(theta); plot(x, y) xlabel('x (m)') ylabel('y (m)') title('part (c) Path of the Slider') %=========end of script 3=============== KINEMATICS OF PARTICLES 51 52 CH. 2 KINEMATICS OF PARTICLES 2.5 Sample Problem 2/10 (Polar Coordinates) A tracking radar lies in the vertical plane of the path of a rocket which is coasting in unpowered flight above the atmosphere. For the instant when θ = 30°, the tracking data give r = 25(104) feet, r& = 4000 ft/s, and θ& = 0.8 deg/s. Let this instant define the initial conditions at time t = 0 and plot vr and vθ as a function of time for the next 150 seconds. You may assume that g remains constant at 31.4 ft/s2 during this time interval. Problem Formulation Place a Cartesian coordinate system at the radar with x positive to the right and y positive up. Since the rocket is coasting in unpowered flight we can use the equations for projectile motion. x = x 0 + v 0 cos(β )t y = y 0 + v 0 sin (β )t − 1 2 gt 2 Where x0 = 25(104)sin(30) ft, y0 = 25(104)cos(30) ft, v0 is the initial speed (5310 ft/sec, see the sample problem) and β is the angle that v0 makes with the horizontal. From the figure shown to the right we can find the angle between v0 and the r axis as φ = tan −1 (3490 / 4000) = 41.11o . Since the r axis is 60° from the horizontal, β = 60 – 41.11 = 18.89°. With r and θ defined as in the sample problem we have, at any time t r = x2 + y2 θ = tan −1 ( x / y ) Now we find vr and vθ from their definitions. v r = r& vθ = rθ& Substitution of x and y into the above equations and carrying out the derivatives with respect to time gives vr and vθ as functions of time. The results are very KINEMATICS OF PARTICLES 53 messy and will not be given here. Remember, though, that substitutions such as this can be made automatically when using computer software such as MATLAB. MATLAB Scripts %%%%%%%%%%%%%%%%%%%% Script # 1 %%%%%%%%%%%%%%%%%%%%%%%%% % This script obtains symbolic results for the r and % theta components of the velocity syms x y t x0 y0 v0 beta g x = x0+v0*cos(beta)*t y = y0+v0*sin(beta)*t-1/2*g*t^2 r = sqrt(x^2+y^2) theta = atan(x/y) vr = diff(r,t); vtheta = r*diff(theta,t); % The results from these operations will be copied and % pasted into script #2 for plotting. Before they can be % used they have to be "vectorized", which places periods % in front of the operators to insure term by term rather % than matrix operations. We'll go ahead and do that here % and then copy and paste the vectorized results. vr = vectorize(vr) vtheta = vectorize(vtheta) Output from Script #1 x = x0+v0*cos(beta)*t y = y0+v0*sin(beta)*t-1/2*g*t^2 r = ((x0+v0*cos(beta)*t)^2+(y0+v0*sin(beta)*t-1/2*g*t^2)^2)^(1/2) theta = atan((x0+v0*cos(beta)*t)/(y0+v0*sin(beta)*t-1/2*g*t^2)) vr = 1./2./((x0+v0.*cos(beta).*t).^2+(y0+v0.*sin(beta).*t1./2.*g.*t.^2).^2).^(1./2).*(2.*(x0+v0.*cos(beta).*t).*v0.*cos(beta)+2.*(y0+v0.* sin(beta).*t-1./2.*g.*t.^2).*(v0.*sin(beta)-g.*t)) vtheta = ((x0+v0.*cos(beta).*t).^2+(y0+v0.*sin(beta).*t1./2.*g.*t.^2).^2).^(1./2).*(v0.*cos(beta)./(y0+v0.*sin(beta).*t-1./2.*g.*t.^2)(x0+v0.*cos(beta).*t)./(y0+v0.*sin(beta).*t-1./2.*g.*t.^2).^2.*(v0.*sin(beta)g.*t))./(1+(x0+v0.*cos(beta).*t).^2./(y0+v0.*sin(beta).*t-1./2.*g.*t.^2).^2) 54 CH. 2 KINEMATICS OF PARTICLES %%%%%%%%%%%%%%%%%%%% Script # 2 %%%%%%%%%%%%%%%%%%%%% % This script plots the r and theta components of the % velocity as functions of time. x0 = 25*10^4*sin(30*pi/180); y0 = 25*10^4*cos(30*pi/180); v0 = 5310; beta = 18.89*pi/180; g = 31.4; t = 0:0.1:150; % The following expressions are copied and pasted from the % results obtained with script #1. vr = 1./2./((x0+v0.*cos(beta).*t).^2+(y0+v0.*sin(beta).*t1./2.*g.*t.^2).^2).^(1./2).*(2.*(x0+v0.*cos(beta).*t).*v0.* cos(beta)+2.*(y0+v0.*sin(beta).*t-1./2.*g.*t.^2).* (v0.*sin(beta)-g.*t)); vtheta = ((x0+v0.*cos(beta).*t).^2+(y0+v0.*sin(beta).*t1./2.*g.*t.^2).^2).^(1./2).*(v0.*cos(beta)./(y0+v0.*sin(bet a).*t-1./2.*g.*t.^2)-(x0+v0.*cos(beta).*t)./ (y0+v0.*sin(beta).*t-1./2.*g.*t.^2).^2.*(v0.*sin(beta)g.*t))./(1+(x0+v0.*cos(beta).*t).^2./(y0+v0.*sin(beta).*t1./2.*g.*t.^2).^2); plot(t, vr, t, vtheta) xlabel('t (sec)') ylabel('velocity (ft/sec)') KINEMATICS OF PARTICLES 55 2.6 Problem 2/183 (Space Curvilinear Motion) The base structure of the firetruck ladder rotates about a vertical axis through O with a constant angular velocity θ& = Ω. At the same time, the ladder unit OB elevates at a constant rate φ& = Ψ, and section AB of the ladder extends from within section OA at the constant rate R& = Λ. Find general expressions for the components of acceleration of point B in spherical coordinates if, at time t = 0, θ = 0, φ = 0, and AB = 0. Express your answers in terms of Ω, Ψ, Λ, R0 and t, where R0 = OA and is constant. Plot the components of acceleration of B as a function of time for the case Ω =10 deg/s, Ψ = 7 deg/s, Λ = 0.5 m/s, and R0 = 9 m. Let t vary between 0 and the time at which φ = 90°. Problem Formulation The components of acceleration in spherical coordinates are, && − Rφ& 2 − Rθ& 2 cos 2 φ aR = R ( ) aθ = cos φ d 2 & R θ − 2 Rθ&φ& sin φ R dt aφ = 1 d 2& R φ + Rθ& 2 sin φ cos φ R dt ( ) The components may be obtained as functions of time by substituting, R = R0 + Λt , θ = Ωt and φ = Ψt Differentiation and substitution will be performed in MATLAB. The results are, ( ) a R = (R0 + Λt ) Ψ 2 − Ω 2 cos 2 (Ψt ) aθ = 2ΩΛ cos(Ψt ) − 2ΩΨ (R0 + Λt ) sin(Ψt ) 56 CH. 2 KINEMATICS OF PARTICLES aφ = 2ΨΛ + (R0 + Λt )Ω 2 sin (Ψt ) cos(Ψt ) MATLAB Scripts %%%%%%%%%%%%%%%%%%%%% Script #1 %%%%%%%%%%%%%%%%%%%%%%%%%% % This script calculates the components of the % acceleration symbolically % O = Omega; P = Phi; L = Lambda syms O P L t R0 R = R0+L*t; theta = O*t; phi = P*t; % Even though there are some obvious simplifications % in this case, we still write the most general % expressions for the spherical components of the % acceleration. In this way we can consider other types % of time dependence without modifying the script. a_R = diff(R,t,2)-R*diff(phi,t)^2-R*diff(theta,t)^2 *cos(phi)^2 a_theta = cos(phi)/R*diff(R^2*diff(theta,t),t)2*R*diff(theta,t)*diff(phi,t)*sin(phi) a_phi = 1/R*diff(R^2*diff(phi,t),t)+R*diff(theta,t)^2 *sin(phi)*cos(phi) %%%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%%%%%%%%%%%%% Output of script #1 a_R = -(R0+L*t)*P^2-(R0+L*t)*O^2*cos(P*t)^2 a_theta = 2*cos(P*t)*O*L-(2*R0+2*L*t)*O*P*sin(P*t) a_phi = 2*P*L+(R0+L*t)*O^2*sin(P*t)*cos(P*t) %%%%%%%%%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%%%%%%%%%%%% % This script plots the components of the acceleration % as functions of time % O = Omega; P = Phi; L = Lambda L = 0.5; O = 10*pi/180; P = 7*pi/180; R0 = 9; tf = pi/2/P; % time at which phi=pi/2 t=0:0.05:tf; a_R = -(R0+L*t)*P^2-(R0+L*t).*O^2.*cos(P*t).^2; a_theta = 2*cos(P*t)*O*L-(2*R0+2*L*t)*O*P.*sin(P*t); a_phi = 2*P*L+(R0+L*t)*O^2.*sin(P*t).*cos(P*t); plot(t,a_R,t,a_theta,t,a_phi) KINEMATICS OF PARTICLES 57 xlabel('time (sec)') title('acceleration (m/s^2)') %%%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%%%%%%%%%%%%% 2.7 Sample Problem 2/16 (Constrained Motion of Connected Particles) The tractor A is used to hoist the bale B with the pulley arrangement shown. If A has a forward velocity vA, determine an expression for the upward velocity vB of the bale in terms of x. Put the result in nondimensional form by introducing the velocity ratio η = vB/vA and nondimensional position χ = x/h. Plot η versus χ for 0 ≤ χ ≤ 2. Problem Formulation The length L of the cable can be written L = 2(h − y ) + l + cons tan ts = 2(h − y ) + h 2 + x 2 + cons tan ts Now, L& = 0 will be used to obtain a relation between vA (= x& ) and vB (= y& ). L& = 0 = −2 y& + xx& h +x 2 2 vB = 1 2 xv A h2 + x2 58 CH. 2 KINEMATICS OF PARTICLES The nondimensional result is now obtained by substituting v B = ηv A and x = χh . η= 1 χ 2 1+ χ 2 Even though these operations are rather easily performed by hand, it is instructive to have MATLAB do them. In particular, it will be instructive to see how to evaluate L& even though x and y are not known explicitly as functions of time. MATLAB Worksheet and Script EDU» syms t h vA vB chi eta EDU» x = sym('x(t)'); y = sym('y(t)'); EDU» L = 2*(h-y)+sqrt(h^2+x^2) L= 2*h-2*y(t)+(h^2+x(t)^2)^(1/2) Note that we need to differentiate L with respect to time. Both x and y depend on time, however, exactly how they depend on time is not known. It turns out that this is not a problem. All we need to do is let MATLAB know x and y depend on time by writing x = sym('x(t)') and y = sym('y(t)'). To illustrate, consider the following. EDU» diff(x, t) ans = diff(x(t),t) This derivative would have been evaluated as zero had we not declared x to be a function of t. Now let’s proceed by differentiating L with respect to t. We give the result a name (eqn) to facilitate later substitutions. EDU» eqn = diff(L,t) eqn = -2*diff(y(t),t)+1/(h^2+x(t)^2)^(1/2)*x(t)*diff(x(t),t) EDU» eqn = subs(eqn,diff(y,t),eta*vA) % substitutes vB = ηvA for y& (diff(y,t)) eqn = -2*eta*vA+1/(h^2+x(t)^2)^(1/2)*x(t)*diff(x(t),t) EDU» eqn = subs(eqn,diff(x,t),vA) eqn = KINEMATICS OF PARTICLES 59 -2*eta*vA+1/(h^2+x(t)^2)^(1/2)*x(t)*vA EDU» eqn = subs(eqn,x,chi*h) eqn = -2*eta*vA+1/(h^2+chi^2*h^2)^(1/2)*chi*h*vA Now let’s remember that eqn is just a name for dL/dt which is zero. We now solve this equation for η (eta). Also recall that solve(eqn, eta) actually solves the equation eqn = 0 for eta. EDU» eta = solve(eqn, eta) eta = 1/2*chi*h/(h^2+chi^2*h^2)^(1/2) We note finally that the h cancels in the above expression yielding the result given in the problem formulation section above. Now we can produce the required plot. %%%%%%% Script for plotting chi versus eta %%%%% chi = 0:0.01:2; eta = 1/2*chi./sqrt(1+chi.^2); plot(chi,eta) xlabel('x/h') title('v_B/v_A') %%%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%% KINETICS OF PARTICLES The kinetics of particles is concerned with the motion produced by unbalanced forces acting on a particle. This chapter considers three approaches to the solution of particle kinetics problems: (1) direct application of Newton’s second law, (2) work and energy, and (3) impulse and momentum. Problem 3.1 is a rectilinear motion problem where solve is used solve three equations symbolically for three unknowns. In problem 3.2, dsolve is used to solve a second order differential equation with initial conditions. The absolute path of a particle is then plotted using polar. Problem 3.3 uses MATLAB to study the effect of initial spring stretch upon the velocity of a slider. A physical interpretation of the results is also required. Problem 3.4 is a typical ballistic pendulum problem requiring both work/energy and conservation of momentum to relate the velocity of a projectile to the maximum swing angle of a pendulum. Problem 3.5 is a relatively straightforward conservation of angular momentum problem where MATLAB is used to generate a plot that might be useful in a parametric study. In problem 3.6, two equations are solved symbolically for two unknowns using solve. The maximum value of a function is then determined using diff and solve. 3 62 CH. 3 KINETICS OF PARTICLES 3.1 Sample Problem 3/3 (Rectilinear Motion) The 250-lb concrete block A is released from rest in the position shown and pulls the 400-lb log up the 30° ramp. Plot the velocity of the block as it hits the ground at B as a function of the coefficient of kinetic friction µk between the log and the ramp. Let µk vary between 0 and 1. Why does the computer not plot results for the entire range specified? Problem Formulation The constant length of the cable is L = 2sC + sA (see figure). Differentiating this expression twice yields a relation between the acceleration of A and C (note that aC = aLOG). 0 = 2aC + a A (1) From the free-body diagram for the log [ΣF y = ma y = 0 [ΣFx = ma x ] ] N − 400 cos(30) = 0 µ k N − 2T + 400 sin(30) = 400 aC 32.2 Substituting N yields, 400µ k cos(30) − 2T + 400 sin(30) = 400 aC 32.2 (2) From the free-body diagram for block A [↓ ΣF = ma] 250 − T = 250 aA 32.2 (3) MATLAB will be used to solve the three equations above for aA, aC and T in terms of µk. Since the accelerations are constant, v 2A = 2a A d where d is the vertical distance through which block A has fallen. Thus, the velocity of A when it strikes the ground (d = 20 ft) is KINETICS OF PARTICLES 63 v Af = 40a A = v B MATLAB Scripts %%%%%%%%%%%%%%%%% Script #1 %%%%%%%%%%%%%%%%%%%%%%%%% % This script solves the three equations developed % in the problem formulation section for the two % accelerations (aC and aA) and the tension T. % When solving multiple equations, it is a good idea % to let the variables solved for be single characters. % Thus, we set x = aC, y = aA, z = T. eqn1 = '2*x+y = 0'; eqn2 = '400*muk*cos(theta)-2*z+400*sin(theta)=400/32.2*x'; eqn3 = '250-z = 250/32.2*y'; [x,y,z]=solve(eqn1,eqn2,eqn3) %%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%%%%%%%%%% Output from Script #1 x = 7.9674*muk - 6.9 y = -15.9349*muk + 13.8 z = 142.857+123.7179*muk The only thing we are interested in here is aA (y in our solution above). Thus, a A = 13.8 − 15.9349 µ k Note that the accelerations may be either positive or negative depending on the value of µk . The largest value of µk for which the block will move up can thus be found by solving the equation aA = 0 for µk . This yields µk = 13.8/15.935 = 0.866. %%%%%%%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%%%%%%%%%%%% % This script plots the vB (the velocity of the % block as it hits the ground) versus the friction % coeficient mu_k muk = 0:0.01:1; aA=-15.9349*muk+13.8; vB = sqrt(40*aA) plot(muk,vB) xlabel('mu_k') ylabel('v_B (ft/sec)') %%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%%%%%%%%%%%%% 64 CH. 3 KINETICS OF PARTICLES At first sight, it seems that no results are plotted beyond the limiting value for µk (0.866) that was determined above. If the plot were in color, you would see that MATLAB has actually plotted zeros beyond this point. From a numerical point of view this occurs because MATLAB will plot only the real part of complex numbers. If you have MATLAB print the values of vB you will find that the real parts of all the complex numbers generated are zero. Whenever imaginary or complex values result there is usually some physical explanation. In this problem, the physical explanation is that the log will not slide up the incline if the coefficient of friction is too large. KINETICS OF PARTICLES 65 3.2 Problem 3/98 (Curvilinear Motion) The particle P is released at time t = 0 from the position r = r0 inside the smooth tube with no velocity relative to the tube, which is driven at the constant angular velocity ω0 about the vertical axis. Determine the radial velocity vr, the radial position r, and the transverse velocity vθ as functions of time t. Plot the absolute path of the particle during the time that it is inside the tube for r0 = 0.1 m, l = 1 m, and ω0 = 1 rad/s. Problem Formulation From the free-body diagram to the right, ( ΣFr = 0 = ma r = m r&& − rθ& 2 ) &r& = rθ& 2 = rω 02 Any book on differential equations will have the solution to this equation in terms of the hyperbolic sine and cosine, r = A sinh(ω 0 t ) + B cosh(ω 0 t ) The constants A and B are found from the initial conditions. These conditions are that r = r0 and r& = 0 at t = 0. The second condition comes from the fact that the particle has no velocity (initially) relative to the tube. Before evaluating this condition we must first differentiate r with respect to time. r& = Aω 0 cosh(ω 0 t ) + Bω 0 sinh(ω 0 t ) r (t = 0) = r0 = A sinh(0) + B cosh(0) = B r&(t = 0) = 0 = Aω 0 cosh(0) + Bω 0 sinh(0) = Aω 0 From the above we have B = r0 and A = 0. Thus, r = r0 cosh(ω 0 t ) From this we can obtain the radial and transverse velocities, 66 CH. 3 KINETICS OF PARTICLES v r = r& = r0ω 0 sinh(ω 0 t ) vθ = rθ& = r0ω 0 cosh(ω 0 t ) The absolute path of the particle will be graphed using polar plotting. For this we need r as a function of θ. Since θ = ω0t we have, r = r0 cosh(θ ) We want to plot this function only up to the point where the particle leaves the tube. Substituting r = 1 we have 1 = 0.1cosh(θ), or θ = cosh-1(10) = 2.993 rads. Thus, the particle leaves the tube when θ = 2.993 rads (171.5°). MATLAB Worksheet EDU» theta = 0:0.01:2.993; EDU» r = 0.1*cosh(theta); EDU» polar(theta,r) KINETICS OF PARTICLES 67 3.3 Sample Problem 3/17 (Potential Energy) The 10-kg slider A moves with negligible friction up the inclined guide. The attached spring has a stiffness of 60 N/m and is stretched δ m at position A, where the slider is released from rest. The 250-N force is constant and the pulley offers negligible resistance to the motion of the cord. Plot the velocity of the slider as it passes C as a function of the initial spring stretch δ. Let δ vary between –0.4 and 0.8 m and explain the results when δ exceeds a value of about 0.65 m. Problem Formulation The change in the elastic potential energy is ∆Ve = ( ) ( 1 1 2 k x22 − x12 = k (1.2 + δ ) − δ 2 2 2 ) The other results in the sample problem are unchanged, ( ) 1 1 m v 2 − v 02 = (10)v 2 2 2 ∆V g = mg∆h = 10(9.81)(1.2 sin 30) = 58.9 J U 1′− 2 = 250(0.6) = 150 J U 1′− 2 = 150 = ∆T = ( 1 (10)v 2 + 58.9 + 1 (60 ) (1.2 + δ )2 − δ 2 2 2 ) This equation can be solved for v either by hand or by using MATLAB. The result is v= 1 958 − 1440δ 10 MATLAB Worksheet and Script First we will solve the work/energy equation symbolically for v. EDU» eqn='150=5*v^2+58.9+30*((1.2+del)^2-del^2)' 68 CH. 3 KINETICS OF PARTICLES eqn = 150=5*v^2+58.9+30*((1.2+del)^2-del^2) EDU» solve(eqn,'v') ans = .1*(958.-1440.*del)^(1/2) -.1*(958.-1440.*del)^(1/2) %%%%%%%%%%%%%%%% Script %%%%%%%%%%%%%%%%% % This script plots the velocity of the slider % at C versus the initial spring stretch del del = -0.4:0.01:0.8; v = 1/10*sqrt(958-1440*del); plot(del,v) axis([-0.4 0.8 0 4]) xlabel('initial spring stretch (m)') ylabel('velocity (m/s)') %%%%%%%%%%%%% end of script %%%%%%%%%%%%%%% KINETICS OF PARTICLES 69 At first sight, it seems that no results are plotted beyond δ ≅ 0.65 m. If the plot were in color, you would see that MATLAB has actually plotted zeros beyond this point. Why? One reason is to observe from the above equation that v becomes imaginary when δ > 958/1440 = 0.665 m. MATLAB interprets imaginary numbers as complex numbers with a real part equal to zero. When asked to plot a complex number it will plot only the real part. Thus, we have zeros plotted after δ = 0.665 m. But this is a numerical reason instead of a physical explanation. Usually, imaginary answers signify a situation that is physically impossible for some reason. One way of understanding this is as follows. If the spring is initially compressed it will, at least for some part of the motion, be pushing up and thus be aiding the 250 N force in overcoming the weight of the slider. If the spring is initially stretched, it will always be pulling back on the slider. Thus the 250 N force will have to overcome not only the weight but also the spring force. It stands to reason then that there will be some value for the initial spring stretch beyond which the 250 N force will not be able to pull the slider all the way to C. This value is found from the limiting case where v = 0. Thus, the block never reaches C if δ > 0.665 m. 70 CH. 3 KINETICS OF PARTICLES 3.4 Problem 3/218 (Linear Impulse/Momentum) The ballistic pendulum is a simple device to measure the projectile velocity v by observing the maximum angle θ to which the box of sand with embedded particle swings. As an aid for a laboratory technician, make a plot of the velocity v in terms of the maximum angle θ. Assume that the weight of the box is 50-lb while the weight of the projectile is 2-oz. Problem Formulation (1) Impulse/Momentum During impact, ∆G = 0 and G1 = G 2  2 / 16 + 50   2 / 16   50  v b   ( 0) =  v +   32.2   32.2   32.2  v = 401v b where v is the velocity of the projectile while vb is the velocity of the box of sand immediately after impact. (2) Work/Energy Now we use the work/energy equation with our initial position being the position where the pendulum is still vertical (θ = 0) and the final position is that where the pendulum has rotated through the maximum angle θ. U 1′− 2 = 0 = ∆T + ∆V g = ( ) 1 m 0 2 − v b2 + mg∆h 2 where m is the combined mass of the box and the projectile. v b = 2 g∆h = 2(32.2)(6)(1 − cos θ ) v = 401vb = 7882 1 − cos θ KINETICS OF PARTICLES 71 MATLAB Worksheet EDU» theta = 0:0.01:pi/2; EDU» v = 7882*sqrt(1-cos(theta)); EDU» plot(theta*180/pi,v) EDU» xlabel('theta (degrees)') EDU» title('velocity of projectile (ft/s)') 72 CH. 3 KINETICS OF PARTICLES 3.5 Problem 3/250 (Angular Impulse/Momentum) The assembly of two 5-kg spheres is rotating freely about the vertical axis at 40 rev/min with θ = 90°. The force F that maintains the given position is increased to raise the base collar and reduce the angle from 90° to an arbitrary angle θ. Determine the new angular velocity ω and plot ω as a function of θ for 0 ≤ θ ≤ 90°. Assume that the mass of the arms and collars is negligible. Problem Formulation Since the summation of moments about the vertical axis is zero we have conservation of angular momentum about that axis. The spheres are rotating in a circular path about the vertical axis. The angular momentum of a particle moving in a circular path of radius r with angular velocity ω is H = mr 2ω . Thus, from the conservation of angular momentum we have, 2mr02ω 0 = 2mr 2ω where ω= r0 = 0.1 + 0.6 cos( 45 0 ) r02 r2 ω0 and MATLAB Script %%%%%%%%%%% script %%%%%%%%%%%%%%% theta = 0:0.05:pi/2; r0 = 0.1+0.6*cos(pi/4); r = 0.1+0.6*cos(theta/2); w0 = 40*2*pi/60; w = r0^2./r.^2*w0; plot(theta*180/pi, w) xlabel('theta (degrees)') ylabel('omega (rad/s)') %%%%%%%% end of script %%%%%%%%%%% r = 0.1 + 0.6 cos(θ / 2) KINETICS OF PARTICLES 73 This diagram “begins” at θ = 90° where ω = ω0 = 40(2π)/60 = 4.19 rad/s. 3.6 Problem 3/365 (Curvilinear Motion) The 26-in. drum rotates about a horizontal axis with a constant angular velocity Ω = 7.5 rad/sec. The small block A has no motion relative to the drum surface as it passes the bottom position θ = 0. Determine the coefficient of static friction µs that would result in block slippage at an angular position θ; plot your expression for 0 ≤ θ ≤ 180°. Determine the minimum required coefficient value µmin that would allow the block to remain fixed relative to the drum throughout a full revolution. For a friction coefficient slightly less than µmin, at what angular position θ would slippage occur? Problem Formulation From the free body and mass acceleration diagrams, [ΣFn = man ] N − mg cos θ = mrΩ 2 [ΣFt = mat ] F − mg sin θ = 0 74 CH. 3 KINETICS OF PARTICLES For impending slip we have F = µ s N . Substituting F into the above and solving gives, µs = sin θ g sin θ = 2 1 . 8925 + cos θ g cos θ + rΩ The last two questions can be answered only after plotting µs as a function of θ. MATLAB Worksheet and Scripts %%%%%%%%%%%%%%%% Script #1 %%%%%%%%%%%%%%%%%%%% % This script solves our two equations symbolically % for mu_s and N % O = Omega % x = mu_s % y = N syms theta O g m r x y eqn1 = y-m*g*cos(theta)-m*r*O^2; eqn2 = x*y-m*g*sin(theta); % remember that we write our equations in the form % expression = 0 and then omit the "=0" [x,y] = solve(eqn1,eqn2) %%%%%%%%%%%%%%end of script %%%%%%%%%%%%%%%%%%%% Output of script #1 x= g*sin(theta)/(g*cos(theta)+r*O^2) y= m*g*cos(theta)+m*r*O^2 %%%%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%%%%% % This script plots mu_s as a function of theta theta = 0:0.02:pi; mu_s = sin(theta)./(1.8925+cos(theta)); plot(theta*180/pi, mu_s) xlabel('theta (degrees)') title('coefficient of static friction') %%%%%%%%%%%%%%end of script %%%%%%%%%%%%%%%%%%%% KINETICS OF PARTICLES 75 If the block is not to slip at any angle θ, the coefficient of friction must be greater than or equal to any value shown on the plot above. Thus, the minimum required coefficient value µmin that would allow the block to remain fixed relative to the drum throughout a full revolution is equal to the maximum value in the plot above. The location where this maximum occurs can be found by solving the equation dµ s / dθ = 0 for θ. This θ can then be substituted into µs to yield the required value for µmin. Here’s how we do this with MATLAB. EDU» syms theta EDU» mu_s = sin(theta)./(1.8925+cos(theta)); EDU» dmu = diff(mu_s,theta) dmu = cos(theta)/(757/400+cos(theta))+sin(theta)^2/(757/400+cos(theta))^2 EDU» theta_m = solve(dmu,theta) theta_m = -atan(1/302800*236697316401^(1/2))+pi atan(1/302800*236697316401^(1/2))-pi EDU» eval(theta_m) ans = 2.1275 -2.1275 EDU» mu_min = subs(mu_s,theta,2.1275) mu_min = 0.6224 From the above we see that µmin = 0.622. If µs is slightly less than this value, the block will slip when θ = 2.128 rads (121.9°). KINETICS OF SYSTEMS OF PARTICLES This chapter concerns the extension of principles covered in chapters two and three to the study of the motion of general systems of particles. The chapter first considers the three approaches introduced in chapter 3 (direct application of Newton’s second law, work/energy, and impulse/momentum) and then moves to other applications such as steady mass flow and variable mass. Problem 4.1 considers an application of the conservation of momentum to a system comprised of a small car and an attached rotating sphere. MATLAB is used to plot the velocity of the car as a function of the angular position of the sphere. The absolute position of the sphere is also plotted. Problem 4.2 uses the concept of steady mass flow to study the effects of geometry upon the design of a sprinkler system. One of the main purposes of this problem is to illustrate how a problem can be greatly simplified using non-dimensional analysis. In particular, an equation containing seven parameters is reduced to a non-dimensional equation with only three parameters. Problem 4.3 is a variable mass problem in which MATLAB is used to integrate the kinematic equation vdv = adx . 4 78 CH. 4 KINETICS OF SYSTEMS OF PARTICLES 4.1 Problem 4/26 (Conservation of Momentum) The small car, which has a mass of 20 kg, rolls freely on the horizontal track and carries the 5kg sphere mounted on the light rotating rod with r = 0.4 m. A geared motor drive maintains a constant angular speed θ& = 4 rad/s of the rod. If the car has a velocity v = 0.6 m/s when θ = 0, plot v as a function of θ for one revolution of the rod. Also plot the absolute position of the sphere for two revolutions of the rod. Neglect the mass of the wheels and any friction. Problem Formulation Since ΣFx = 0 we have conservation of momentum in the x direction. The diagram to the right shows the system at θ = 0 and at an arbitrary angle θ. From the relative velocity equation, the velocity of the sphere is the vector sum of the velocity of the car (v) and the velocity of the sphere relative to the car ( rθ& ). (G x )θ =0 = 20(0.6) + 5(0.6) = 15 (G x )θ ( N•s ) = 20v + 5 v − rθ& sin θ = 25v − 8 sin θ Setting (G x )θ =0 = (G x )θ and solving yields, v = 0.6 + 0.32 sin θ Now let time t = 0 be the time when θ = 0 and place an x-y coordinate system at the center of the car as shown in the diagram so that x(t) is the position of the center of the car. Since v = dx/dt and θ = 4t we have, t t ∫ ∫ (0.6 + 0.32 sin (4t ))dt = 0.6t + 0.08(1 − cos(4t )) 0 0 x = vdt = The x and y components of the sphere can now be determined as, x s = x + r cos θ = 0.08 + 0.6t + 0.32 cos(4t ) KINETICS OF SYSTEMS OF PARTICLES 79 y s = r sin θ = 0.4 sin(4t ) The absolute position of the sphere can be obtained by plotting ys versus xs. The time required for two revolutions of the arm is 4π/4 = π seconds. MATLAB Scripts %%%%%%%%%%%%%%%% Script #1 %%%%%%%%%%%%%%%%% % This script plots v as a function of theta theta = 0:0.02:2*pi; v = 0.6 + 0.32*sin(theta); plot(theta*180/pi, v) xlabel('theta (degrees)') ylabel('v (m/s)') %%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%% %%%%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%% % This script plots the position of the sphere t = 0:0.01:pi; xs = 0.08+0.6*t+0.32*cos(4*t); ys = 0.4*sin(4*t); plot(xs, ys) xlabel('x (m)') ylabel('y (m)') title('position of the sphere') %%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%% 80 CH. 4 KINETICS OF SYSTEMS OF PARTICLES 4.2 Problem 4/62 (Steady Mass Flow) The sprinkler is made to rotate at the constant angular velocity ω and distributes water at the volume rate Q. Each of the four nozzles has an exit area A. Write an expression for the torque M on the shaft of the sprinkler necessary to maintain the given motion. Here we would like to study the effects of the geometry of the sprinkler upon this torque. To this end, it is helpful to introduce the non-dimensional parameters M′ = M/4ρAru2, Ω = ωr/u, and β = b/r where u = Q/4A is the velocity of the water relative to the nozzle and ρ is the density of the water. Plot the non-dimensional torque M′ versus β for Ω = 0.5, 1, and 2. Let Ω0 be the non-dimensional velocity Ω at which the sprinkler will operate with no applied torque. Plot Ω0 versus β. For both plots let β range between 0 and 1. Problem Formulation The figure to the right shows the three components of the absolute velocity of the water at the exit. u (= Q/4A) is the velocity of the water relative to the nozzle. The mass flow rate m′ = ρQ. Taking clockwise as positive, the application of equation 4/19 of your text yields, ( ΣM 0 = − M = ρQ ωr 2 + ωb 2 − ur − 0 ) KINETICS OF SYSTEMS OF PARTICLES 81 ( ( M = ρQ ur − ω r 2 + b 2 )) Now we want to introduce the non-dimensional parameters defined in the problem statement. For many undergraduate students, non-dimensional analysis is a very confusing topic. It is important to realize that the difficulty is really that of determining which non-dimensional parameters are appropriate for a particular problem. If these parameters have already been defined, as in this problem, all you have to do is substitute. In this case we merely substitute M = 4ρAru2 M′, ω = Ωu/r, and b = rβ into the equation above. When this is done many terms will cancel yielding, ( M ′ =1− Ω 1+ β 2 ) Setting M′ = 0 we can solve for Ω0, Ω0 = 1 1+ β 2 MATLAB Scripts %%%%%%%%%%%%%%%% Script #1 %%%%%%%%%%%%%%%%% % This script plots the non-dimensional torque % M' as a function of beta = b/r for Omega = % 0.5, 1, and 2 beta = 0:0.01:1; Mp = inline('1-Omega*(1+beta.^2)') plot(beta, Mp(0.5,beta),beta, Mp(1,beta),beta, Mp(2,beta)) xlabel('beta = b/r') title('non-dimensional torque') %%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%% 82 CH. 4 KINETICS OF SYSTEMS OF PARTICLES %%%%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%% % This script plots the non-dimensional angular % velocity Omega0 as a function of beta = b/r beta = 0:0.01:1; Omega0 = 1./(1+beta.^2); plot(beta, Omega0) xlabel('beta = b/r') %%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%% KINETICS OF SYSTEMS OF PARTICLES 83 4.3 Problem 4/86 (Variable Mass) The open-link chain of length L and mass ρ per unit length is released from rest in the position shown, where the bottom link is almost touching the platform and the horizontal section is supported on a smooth surface. Friction at the corner guide is negligible. Determine (a) the velocity v1 of end A as it reaches the corner and (b) its velocity v2 as it strikes the platform. Plot v1 and v2 as functions of h for L = 5 m. Problem Formulation Let x be the displacement of the chain and T be the tension in the chain at the corner as shown in the diagram to the right. Note that the acceleration of the horizontal and vertical sections are both equal to &x& . For the horizontal section, [ΣFx = ma x ] T = ρ (L − h − x )&x& For the vertical section, [ ↓ ΣF y = ma y ] ρgh − T = ρh&x& Substituting T from the first equation into the second and simplifying gives, &x& = gh L−x Now we use the relation vdv = adx to write, L −h v1 ∫ 0 vdv = ∫ 0 gh dx L−x v1 = 2 gh ln (L / h ) L −h v12 = 2 ∫ 0 gh L dx = 2 gh ln   L−x h 84 CH. 4 KINETICS OF SYSTEMS OF PARTICLES After end A has passed the corner it will be in free-fall until it hits the platform. With y positive down we have vdv = gdy yielding, ( ) 1 2 v 2 − v12 = gh 2 Substituting for v1 and solving, v 2 = 2 gh(1 + ln( L / h) ) MATLAB Worksheet and Script EDU» syms v v1 v2 h g x L L −h v1 First we solve the equation ∫ 0 vdv − ∫ gh / (L − x )dx = 0 for v , omitting (as usual) 1 0 the “=0”. EDU» eqn1 = int(v,v,0,v1)-int(g*h/(L-x),x,0,L-h) eqn1 = 1/2*v1^2+log(h)*g*h-log(L)*g*h EDU» solve(eqn1,v1) ans = (-2*log(h)*g*h+2*log(L)*g*h)^(1/2) -(-2*log(h)*g*h+2*log(L)*g*h)^(1/2) MATLAB has found two solutions. The first is the one we want since it is positive. This solution can easily be simplified to the result given above in the problem formulation section. Once v1 is known it is rather easy to find v2. We’ll do it symbolically here for purposes of illustration. First we copy and paste the first solution above to define v1. Then we solve the equation ∫ v2 v1 ∫ h vdv − gdy = 0 for v2. Note how the result for v1 is automatically 0 substituted. EDU» v1 = (-2*log(h)*g*h+2*log(L)*g*h)^(1/2); EDU» eqn2 = int(v,v,v1,v2)-int(g,x,0,h) eqn2 = 1/2*v2^2+log(h)*g*h-log(L)*g*h-g*h KINETICS OF SYSTEMS OF PARTICLES 85 EDU» solve(eqn2,v2) ans = (-2*log(h)*g*h+2*log(L)*g*h+2*g*h)^(1/2) -(-2*log(h)*g*h+2*log(L)*g*h+2*g*h)^(1/2) Once again, the first solution will simplify to that given in the problem formulation section above. %%%%%%%%%%%%%%%%%%% Script %%%%%%%%%%%%%%%%%%%% % This script plots v1 and v2 as functions of % h for L = 5 m L = 5; g = 9.81; h = 0:0.01:L; v1 = sqrt(2*g*h.*log(L./h)); v2 = sqrt(2*g*h.*(1+log(L./h))); plot(h, v1, h, v2) xlabel('h (m)') ylabel('velocity (m/s)') %%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%% PLANE KINEMATICS OF RIGID BODIES This chapter extends the kinematic analysis of particles covered in chapter 2 to rigid bodies by taking into account the rotational motion of the body. Thus, the motion of rigid bodies involves both translation and rotation. Problem 5.1 is a straightforward problem that compares the angular position of a rotating disk to the total number of revolutions. Problem 5.2 is an interesting application of absolute motion analysis. The problem illustrates the usefulness of nondimensional analysis in conducting a parametric study of the effects of geometry upon the performance of a hydraulic lift. In problem 5.3 the velocity of the piston in a reciprocating engine is plotted versus the angular orientation of the crank. The maximum velocity and the corresponding orientation (of the crank) are obtained using max after first trying diff and solve. The instantaneous center of zero velocity is used in problem 5.4 to relate the velocity of a vertical control rod to the angular velocity of a rotating bar in a switching device. Problem 5.4 is a fairly straightforward relative acceleration problem. Problem 5.5 considers the reciprocating engine of problem 5.3 but uses absolute rather than relative motion analysis. The problem illustrates the use of computer algebra for carrying out some tedious calculations including differentiation and substitution. 5 88 CH. 5 PLANE KINEMATICS OF RIGID BODIES 5.1 Problem 5/3 (Rotation) The angular velocity of a gear is controlled according ω = 12 – 3t2 where ω, in radians per second, is positive in the clockwise sense and where t is the time in seconds. Find the net angular displacement ∆θ and the total number of revolutions N through which the gear turns in terms of the time t. Plot ∆θ and N for time t up to 4 seconds. Problem Formulation You may want to have a look at sample problem 5/1 in your text before continuing with this problem. In particular it is important to note the difference between the angular displacement ∆θ and total number of revolutions N. The angular displacement is the integral of ω over time and can be positive or negative depending on whether the rotation is clockwise or counterclockwise. Referring to the graph of ω to the right, ∆θ will be the total area under the curve up to a particular time. The area is negative when the curve dips below the time axis. As the name implies, the total number of revolutions simply counts the number of times the disk rotates and is not concerned with whether the rotation is clockwise or counterclockwise. Thus, N will be proportional to the magnitude of the area under the angular velocity curve. For example, suppose the gear rotates two revolutions clockwise followed by two revolutions counterclockwise. In this case ∆θ = 0 but N = 4. Based upon the above it is essential to know beforehand the time when ω = 0 (i.e. when the gear changes direction) in order to calculate N. Setting 12 – 3t2 = 0 gives t = 2 seconds. For time greater than 2 sec. we need to break the integral up into two parts when calculating N. For ∆θ we do not need to keep track of when the gear changes directions. Thus, dθ   ω =  dt  t t 0 0 ( ) ∆θ = ∫ ωdt = ∫ 12 − 3u 2 du = 12t − t 3 PLANE KINEMATICS OF RIGID BODIES 89 t t ( ) ( 1 1 ωdt = ∫ 12 − 3t 2 dt = 12t − t 3 ∫ 2π 0 π 2 0 ) for t < 2 sec N= for t > 2 sec N= t  1 2 1 t 1 2 2  − = − − ω ω 12 3 12 − 3t 2 dt  dt dt t dt ∫ ∫ ∫ ∫  2π 0 2π 2 2π  0 2  N= 1 32 − 12t + t 3 2π ( ( ) ( ) ) In summary, ∆θ = 12t − t 3 ( ( all t )  1 3  2π 12t − t N = 1  32 − 12t + t 3  2π )  t < 2 sec  t > 2 sec   Although the expression for N is rather simple, plotting the result can be a little tricky due to the change in the expression at time t = 2 sec. How do we tell MATLAB to stop plotting one expression and start plotting the other? Actually, this is considerably easier to do in MATLAB than in some other applications (e.g. Maple) due to the fact that you are not limited to a single range variable for the ordinate. Thus we can easily set up two time variables t1 and t2 and then write the two expressions in terms of those two variables. MATLAB Scripts %==== script #1 ========================== % This script plots the angular velocity as % a function of time. Note the way a horizontal % line at zero is plotted. This is useful for % visualizing where omega becomes negative t = 0:0.05:4; omega = 12 - 3*t.^2; zero = 0*t; plot(t,omega,t,zero) xlabel('time (sec)') title('angular velocity (rad/s)') %======= end of script =================== 90 CH. 5 PLANE KINEMATICS OF RIGID BODIES %==== script #2 ========================== % This script plots the angular position % (in degrees) as a function of time. t = 0:0.05:4; theta = (12*t - t.^3)*180/pi; zero = 0*t; plot(t,theta,t,zero) xlabel('time (sec)') title('angular position (degrees)') %======= end of script =================== PLANE KINEMATICS OF RIGID BODIES 91 %==== script #3 ========================== % This script plots the number of revolutions % as a function of time. Note that different % time variables are used for the two different % expressions for N t1 = 0:0.05:2; N1 = (12*t1 - t1.^3)/2/pi; t2 = 2:0.05:4; N2 = (32 - 12*t2 + t2.^3)/2/pi; plot(t1,N1,t2,N2) xlabel('time (sec)') title('total number of revolutions') %======= end of script =================== 92 CH. 5 PLANE KINEMATICS OF RIGID BODIES PLANE KINEMATICS OF RIGID BODIES 93 5.2 Problem 5/44 (Absolute Motion) Derive an expression for the upward velocity v of the car hoist system in terms of θ. The piston rod of the hydraulic cylinder is extending at the rate s& . Plot the non-dimensional velocity v / s& as a function of θ for b/L = 0.1, 0.5, 1, and 2. Problem Formulation From the diagram to the right, y = 2b sin θ y& = 2bθ& cos θ If the angular velocity θ& were known as a function of θ we would be finished. The motion of the car hoist system is controlled by the extension rate s& of the hydraulic cylinder rather than the angular velocity. Thus, to complete the problem we need to relate θ& and s& . s 2 = L2 + b 2 − 2 Lb cos θ 2ss& = 0 + 0 + 2 Lbθ& sin θ ss& Lb sin θ θ& = Substituting, v= 2bss& 2s& L2 + b 2 − 2 Lb cos θ cos θ = Lb sin θ L tan θ 2 1 + (b / L ) − 2(b / L ) cos θ 2 v / s& = where β = b/L. tan θ = 2 1 + β 2 − 2 β cos θ tan θ 94 CH. 5 PLANE KINEMATICS OF RIGID BODIES MATLAB Script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % v is the non-dimensional velocity v = inline('2/tan(theta)*sqrt(1^2+beta^2-2*beta*cos(theta))') v = vectorize(v); % the above adds "." to insure term by term rather than % matrix operations th = 10*pi/180:0.01:pi/2; deg = th*180/pi; plot(deg,v(.1,th),deg,v(.5,th),deg,v(1,th),deg,v(2,th)) axis([10 90 0 5]) xlabel('theta (degrees)') title('non-dimensional velocity') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PLANE KINEMATICS OF RIGID BODIES 95 5.3 Sample Problem 5/9 (Relative Velocity) The common configuration of a reciprocating engine is that of the slider crank mechanism shown. If crank OB has a clockwise rotational speed of 1500 rev/min; (a) Plot vA versus θ for one revolution of the crank. (b) Find the maximum speed of the piston A and the corresponding value of θ. Problem Formulation Let l be the length of connecting rod AB. Start with the relative velocity equation vB = vA + vB/A The crank pin velocity is vB = rω and is normal to OB. The velocity of A is horizontal while the velocity of B/A has magnitude lωAB and is directed perpendicular to AB. The angle β can be found in terms of θ by using the law of sines, sin β = r2 r sin θ Also, cos β = 1 − sin 2 β = 1 − 2 sin 2 θ l l From the vector diagram to the right, v B cos θ = v B / A cos β . Thus, v B / A = cos θ rω cos β Also from the diagram, v A = v B sin θ + v B / A sin β . Substitution yields    r cos θ v A = rω (sin θ + cos θ tan β ) = rω sin θ 1 + r2  − l 1 sin 2 θ  2 l         Note that vA has been expressed explicitly in terms of θ by substituting for cosβ and tanβ. This has been done only for sake of clarity. When working with a 96 CH. 5 PLANE KINEMATICS OF RIGID BODIES computer such substitutions will be automatic once β has been defined in terms of θ ( β = sin −1 (r sinθ / l ) ). In part (b) we first try to find the angle θ at which the maximum value of vA occurs by solving the equation dvA/dθ = 0 for θ. Although a solution is obtained it is extremely messy and occupies many screen widths. For this reason we switch to an approximate approach which utilizes Matlab’s max function. This method gives the maximum speed of the piston as vA = 69.6 ft/sec at θ = 72.2°. MATLAB Worksheet and Scripts %%%%%%%%%%%%%%% Script #1 %%%%%%%%%%%% % This script solves for vA using symbolic % algebra syms r L theta omega vB = r*omega; beta = asin(r/L*sin(theta)); vBA = cos(theta)/cos(beta)*r*omega; vA = vB*sin(theta)+vBA*sin(beta) %%%%%%%%%%%%% end of script %%%%%%%%%%% Output of script #1 vA = r*omega*sin(theta)+cos(theta)/(1-r^2/L^2*sin(theta)^2)^(1/2)* r^2*omega/L*sin(theta) %%%%%%%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%%%% % This script plots the velocity of piston A % as a function of theta r = 5/12; L = 14/12; o = 1500*2*pi/60; % o=omega x = 0:0.01:2*pi; % x=theta vA = r*o*sin(x)+cos(x)./(1-r^2/L^2*sin(x).^2).^(1/2)*r^2*o/L.*sin(x); plot(x*180/pi,vA) xlabel('theta (degrees)') title('velocity of piston A (ft/sec)') %%%%%%%%%%%%% end of script %%%%%%%%%%% PLANE KINEMATICS OF RIGID BODIES 97 Part (b) %%%%%%%%%%%%%%%%%%% Script #3 %%%%%%%%%%%%%% % This script solves for the values of theta which % make vA a maximum r = 5/12; L = 14/12; omega = 1500*2*pi/60; syms theta % The following is copied and pasted from the output of % script #1 vA = r*omega*sin(theta)+cos(theta)/(1r^2/L^2*sin(theta)^2)^(1/2)*r^2*omega/L*sin(theta) dvAdtheta = diff(vA,theta) solve(dvAdtheta,theta) %%%%%%%%%%%%%%%%%% end of script %%%%%%%%%%% Output of script #3 vA = 125/6*pi*sin(theta)+625/84*cos(theta)/(1-25/196*sin(theta)^2)^(1/2)* pi*sin(theta) dvAdtheta = 125/6*pi*cos(theta)-625/84*sin(theta)^2/(1-25/196*sin(theta)^2)^(1/2)*pi+ 15625/16464*cos(theta)^2/(1-25/196*sin(theta)^2)^(3/2)*pi*sin(theta)^2+ 625/84*cos(theta)^2/(1-25/196*sin(theta)^2)^(1/2)*pi 98 CH. 5 PLANE KINEMATICS OF RIGID BODIES ans = [output omitted] The output was omitted because it was very messy and extended over several screen widths making it practically unusable. For some reason, these problems do not occur in Matlab 5. Due to these difficulties we try a different approach, using Matlab’s max function. After running script #1 type the following into the command window, EDU>> [vA_max, i] = max(vA) vA_max = 69.5514 i= 127 EDU>> x(i) ans = 1.2600 EDU>> ans*180/pi ans = 72.1927 Thus, the maximum speed of the piston is vA = 69.6 ft/sec at θ = 72.2°. PLANE KINEMATICS OF RIGID BODIES 99 5.4 Problem 5/123 (Relative Acceleration) The two rotor blades of radius r = 800-mm rotate counterclockwise with a constant angular velocity ω about the shaft at O mounted in the sliding block. The acceleration of the block is aO. Determine the magnitude of the acceleration of the tip A of the blade in terms of r, ω, aO, and θ. Plot the acceleration of A as a function of θ for one revolution if aO = 3 m/s. Consider three cases: ω = 2, 4, and 6 rad/s. Problem Formulation The acceleration of A relative to O is r r r r a A = aO + (a A / O )n + (a A / O )t The acceleration of O is to the right while the normal relative acceleration must point from A towards O. Since ω is constant, the tangential relative acceleration will be zero. These considerations lead to the vector diagram shown to the right. Using the law of cosines, a A = aO2 + (a A / O )n2 − 2aO (a A / O )n cos θ ( ) a A = aO2 + rω 2 2 ( ) − 2aO rω 2 cosθ MATLAB Script %======== script ============================== aA = inline('sqrt(3^2+0.8^2*omega^4-2*3*0.8*omega^2*cos(theta))') th = 0:.05:2*pi; plot(th,aA(2,th),th,aA(4,th),th,aA(6,th)) xlabel('theta (radians)') ylabel('acceleration (m/s^2)') %========= end of script ======================= 100 CH. 5 PLANE KINEMATICS OF RIGID BODIES 5.5 Sample Problem 5/15 (Absolute Motion) The common configuration of a reciprocating engine is that of the slider crank mechanism shown. If crank OB has a clockwise rotational speed of 1500 rev/min; (a) Plot vA and vG versus time for two revolutions of the crank. (b) Plot aA and aG versus time for two revolutions of the crank. Problem Formulation This problem appears in sample problems 5/9 and 5/15 in your text. Sample problem 5/9 considers a relative velocity analysis while sample problem 5/15 uses a relative acceleration analysis. Generally speaking, the easiest approach to use with a computer is an absolute motion analysis, provided you have software capable of doing symbolic algebra and calculus such as MATLAB. We will use the present problem to illustrate this approach. PLANE KINEMATICS OF RIGID BODIES 101 We start by using the law of sines (sin(β)/r = sin(θ)/l) to express β as a function of θ r  β = sin −1  sin θ  l   where l is the length of connecting rod AB and β is the angle between AB and the horizontal. Now place an x-y coordinate system at O with x positive to the right and y positive up and write expressions for the coordinates of A and G in terms of θ and β x A = −r cos θ − l cos β x G = −r cos θ − r cos β y G = (l − r ) sin β where r is the distance from B to G (4 in. in the figure). All that is needed to find the velocities vA, vGx, and vGy is to differentiate these expressions with respect to time. The magnitude of the velocity of G is then found from v G = v Gx + vGy 2 2 The accelerations aA, aGx, and aGy are then found by differentiating vA, vGx, and vGy with the magnitude of the acceleration of G being obtained from, aG = aGx + a Gy 2 2 Since we will be differentiating with respect to time, the first thing we will do in the computer program is to define θ as a function of time. Then, when we write the above expressions for β, xA, xG, and yG, the computer will automatically substitute for θ rendering each of these as functions of time. Assuming that θ is initially zero, θ (t ) = ωt = 1500(2π ) t = 157.1t 60 The problem statement asks us to plot versus time for two revolutions (θ = 4π radians) of the crank. The time required for two revolutions is 4π/157.1 = 0.08 sec. 102 CH. 5 PLANE KINEMATICS OF RIGID BODIES MATLAB Worksheet and Scripts %%%%%%%%%%%%%% Script #1 %%%%%%%%%%%%%%%%%% % This script uses symbolic algebra to find % the velocity and acceleration of A and G % w = constant angular velocity % rb = r_bar (4 in. in the figure) syms r rb L w t theta = w*t; beta = asin(r/L*sin(theta)); xA = -r*cos(theta)-L*cos(beta) xG = -r*cos(theta)-rb*cos(beta) yG = (L-rb)*sin(beta) vA = diff(xA,t) vGx = diff(xG,t) vGy = diff(yG,t) vG = sqrt(vGx^2+vGy^2) aA = diff(vA,t) aGx = diff(vGx,t) aGy = diff(vGy,t) aG = sqrt(aGx^2+aGy^2) %%%%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%% Output of Script #1 xA = -r*cos(w*t)-L*(1-r^2/L^2*sin(w*t)^2)^(1/2) xG = -r*cos(w*t)-rb*(1-r^2/L^2*sin(w*t)^2)^(1/2) yG = (L-rb)*r/L*sin(w*t) vA = r*sin(w*t)*w+1/L/(1-r^2/L^2*sin(w*t)^2)^(1/2)*r^2*sin(w*t)*cos(w*t)*w vG = ((r*sin(w*t)*w+rb/(1-r^2/L^2*sin(w*t)^2)^(1/2)*r^2/L^2* sin(w*t)*cos(w*t)*w)^2+(L-rb)^2*r^2/L^2*cos(w*t)^2*w^2)^(1/2) aA = r*cos(w*t)*w^2+1/L^3/(1-r^2/L^2*sin(w*t)^2)^(3/2)*r^4*sin(w*t)^2* cos(w*t)^2*w^2+1/L/(1-r^2/L^2*sin(w*t)^2)^(1/2)*r^2*cos(w*t)^2*w^21/L/(1-r^2/L^2*sin(w*t)^2)^(1/2)*r^2*sin(w*t)^2*w^2 aG = ((r*cos(w*t)*w^2+rb/(1-r^2/L^2*sin(w*t)^2)^(3/2)*r^4/L^4*sin(w*t)^2* cos(w*t)^2*w^2+rb/(1-r^2/L^2*sin(w*t)^2)^(1/2)*r^2/L^2*cos(w*t)^2*w^2- PLANE KINEMATICS OF RIGID BODIES 103 rb/(1-r^2/L^2*sin(w*t)^2)^(1/2)*r^2/L^2*sin(w*t)^2*w^2)^2+(L-rb)^2* r^2/L^2*sin(w*t)^2*w^4)^(1/2) The results for the velocities and accelerations will be used to produce the required plots in the following two scripts. The easiest thing to do is to copy the results from the worksheet into the scripts. Before doing this you should use the vectorize command to place periods at appropriate places in order to insure term by term rather than matrix operations. After running the script type “vectorize(vA)” in the worksheet and then copy and paste the result into the script. Repeat this process for vG, aA, and aG. %%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%%% % This script plots vA and vG versus time r = 5/12; rb = 4/12; % r_bar L = 14/12; w = 1500*2*pi/60; % w=omega t = 0:0.001:0.08; vA = r.*sin(w.*t).*w+1./L./(1-r.^2./L.^2.*sin(w.*t).^2).^ (1./2).*r.^2.*sin(w.*t).*cos(w.*t).*w; vG = ((r.*sin(w.*t).*w+rb./(1-r.^2./L.^2.*sin(w.*t).^2).^ (1./2).*r.^2./L.^2.*sin(w.*t).*cos(w.*t).*w).^2+(Lrb).^2.*r.^2./L.^2.*cos(w.*t).^2.*w.^2).^(1./2); plot(t, vA, t, vG) %%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%% 104 CH. 5 PLANE KINEMATICS OF RIGID BODIES %%%%%%%%%%%%%%%% Script #3 %%%%%%%%%%%%%%%%% % This script plots aA and aG versus time r = 5/12; rb = 4/12; L = 14/12; w = 1500*2*pi/60; % w=omega t = 0:0.001:0.08; aA = r.*cos(w.*t).*w.^2+1./L.^3./(1-r.^2./L.^2.* sin(w.*t).^2).^(3./2).*r.^4.*sin(w.*t).^2.*cos(w.*t).^2.*w. ^2+1./L./(1-r.^2./L.^2.*sin(w.*t).^2).^(1./2).*r.^2.* cos(w.*t).^2.*w.^2-1./L./(1-r.^2./L.^2.*sin(w.*t).^2).^ (1./2).*r.^2.*sin(w.*t).^2.*w.^2; aG = ((r.*cos(w.*t).*w.^2+rb./(1-r.^2./L.^2.* sin(w.*t).^2).^(3./2).*r.^4./L.^4.*sin(w.*t).^2.*cos(w.*t). ^2.*w.^2+rb./(1-r.^2./L.^2.*sin(w.*t).^2).^(1./2). *r.^2./L.^2.*cos(w.*t).^2.*w.^2-rb./(1-r.^2./L.^2.* sin(w.*t).^2).^(1./2).*r.^2./L.^2.*sin(w.*t).^2.*w.^2).^2+( L-rb).^2.*r.^2./L.^2.*sin(w.*t).^2.*w.^4).^(1./2); plot(t, aA, t, aG) %%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%% PLANE KINEMATICS OF RIGID BODIES 105 For further study Suppose you wanted to solve this problem for the case where crank OB has a constant angular acceleration α = 60 rad/s2. It turns out that you can solve this problem using exactly the same approach as above, changing only one line defining the dependence of θ upon time. Assuming that the system starts from rest at θ = 0, the appropriate expression for θ is θ (t) = ½αt2 = 30t2. The accelerations for this case are shown below in case you want to give it a try. PLANE KINETICS OF RIGID BODIES This chapter concerns the motion (translation and rotation) of rigid bodies that results from the action of unbalanced external forces and moments. In problem 6.1, five equations are solved for five unknowns. The problem illustrates an alternative to the blunt (but straightforward) simultaneous solution of multiple equations. Instead, the equations are solved in such a way that there is never more than one unknown. In this way, the results are immediately obtained via automatic substitution, thus avoiding some tedious algebra. The force at the hinge of a pendulum is plotted versus the angular position of the pendulum in problem 6.2. The algebra is rather simple in this case and MATLAB is used primarily for purposes of plotting. Problems 6.3 and 6.4 consider rigid bodies in general plane motion. solve is used in problem 6.3 to solve two equations simultaneously for two unknowns. The maximum acceleration of a point on the rigid body is then obtained by using diff and solve. In this problem, solve finds five solutions to the equation and it is necessary to determine which of these is physically correct. This problem is also interesting since a very natural “guess” of the value for the maximum acceleration turns out to be incorrect. Problem 6.4 is an example of a kinetics problem that also requires some kinematics. The angular acceleration of a bar is determined by summing moments. The kinematic equation ωdω = αdθ is then integrated to obtain the angular velocity. Problem 6.5 is an interesting work and energy problem that is complicated considerably by the fact that a spring is engaged for only part of the motion of a rotating bar. Symbolic algebra simplifies this problem considerably, though it is still rather tedious. It is common in Dynamics to find problems that require a combination of methods for their solution. Problem 6.6 is a good example involving both conservation of momentum and work/energy. 6 108 CH. 6 PLANE KINETICS OF RIGID BODIES 6.1 Sample Problem 6/2 (Translation) The vertical bar AB has a mass of m =150 kg with center of mass G midway between the ends. The bar is elevated from rest at θ = 0 by means of the parallel links of negligible mass, with a couple M applied to the lower link at C. Plot the force at A and at B as functions of θ (between 0 and 60°) for two cases, (a) a constant couple M = 5 kN-m, and (b) a constant angular acceleration α = 5 rad/sec2. Problem Formulation The free-body diagram (FBD) and mass acceleration diagram (MAD) are shown to the right. Since the vertical bar undergoes curvilinear translation, the acceleration of all points on the bar will be identical. Thus, we can obtain the acceleration of G immediately from that of point A which moves in a circular path about C. Also note that BD is a two force member since the mass is negligible. Here we take a somewhat different approach than that in your textbook. The main reason for this is that the case of constant moment and constant angular acceleration require a slightly different approach. The difference between the two approaches is easiest seen by writing all the equations first. From the free-body diagram of the connecting link AC, [ΣM C = 0] M − r At = 0 (1) As pointed out in the sample problem, the force and moment equations are identical to the equilibrium equations whenever the mass is negligible. From the free-body diagram of the vertical bar, PLANE KINETICS OF RIGID BODIES 109 [ΣM A = ma d ] [ΣFt rAB B cos θ = mr ω 2 cos θ rAG + mr α sin θ rAG = mat ] [ΣFn = ma n ] (2) At − mg cosθ = mr α (3) B − An + mg sin θ = mr ω 2 (4) And, finally, the kinematics equation ω θ 0 0 ∫ ωdω = ∫ αdθ (5) At this point we have a total of 5 equations and 6 unknowns (M, α, ω, At, An, and B). In part (a), M is specified while in part (b) α is given. In both cases we will have 5 equations and 5 unknowns; however, it will be necessary to solve no more than one equation at a time provided they are done in the right order. This is what yields a different procedure for parts (a) and (b). Part (a) With M known we can find At from Equation (1) and then substitute the result into Equation (3) to get α as a function of θ. At = α= M r At g − cos θ mr r As usual, notice that there is no need to make an explicit substitution of At into α. The computer will make such substitutions automatically. Now we substitute α into Equation (5) and integrate, θ g 1 2  A ω = ∫  t − cos θ dθ mr r 2  0 g   At θ − sin θ  r   mr ω 2 = 2 Finally, we can substitute into Equations (2) and (4) to find B and then An. 110 CH. 6 PLANE KINETICS OF RIGID BODIES B= ( mr rAG ω 2 cos θ + α sin θ rAB cosθ ) An = B + mg sin θ − mr ω 2 A= An2 + At2 Part (b) With α known instead of M we have to take a different approach, starting with Equations (5) and (3), θ θ 0 0 ω = 2 ∫ αdθ = 2α ∫ dθ = 2αθ 2 At = mg cos θ + mr α Now we can substitute into Equations (1), (2) and (4) to find M, B and then An. M = r At B= ( mr rAG ω 2 cos θ + α sin θ rAB cosθ An = B + mg sin θ − mr ω 2 A= An2 + At2 ) PLANE KINETICS OF RIGID BODIES 111 MATLAB Scripts %========= script for part a =========== m = 0.15; g = 9.81; r = 1.5; rAB = 1.8; rAG = 1.2; rAC = r; M = 5; theta = 0:0.1:60; thr = theta*pi/180; At = M/r; alpha = At/m/r - g/r*cos(thr); omsq = 2*(At*thr/m/r - g/r*sin(thr)); B = m*r*rAG*(omsq.*cos(thr)+alpha.*sin(thr))/rAB./cos(thr); An = B + m*g*sin(thr)-m*r*omsq; A = sqrt(At.^2 + An.^2); plot(theta,B,theta,A) xlabel('theta (degrees)') ylabel('force (kN)') title('part (a): constant M') ======= end of script===================== 112 CH. 6 PLANE KINETICS OF RIGID BODIES %========= script for part b =========== m = 0.15; g = 9.81; r = 1.5; rAB = 1.8; rAG = 1.2; rAC = r; alpha = 5; theta = 0:0.1:60; thr = theta*pi/180; omsq = 2*alpha*thr; At = m*g*cos(thr) + m*r*alpha; M = r*At; B = m*r*rAG*(omsq.*cos(thr)+alpha.*sin(thr))/rAB./cos(thr); An = B + m*g*sin(thr)-m*r*omsq; A = sqrt(At.^2 + An.^2); plot(theta,B,theta,A) xlabel('theta (degrees)') ylabel('force (kN)') title('part (b): constant alpha') ======= end of script===================== PLANE KINETICS OF RIGID BODIES 113 6.2 Sample Problem 6/4 (Fixed-Axis Rotation) The pendulum has a mass of 7.5 kg with a mass center at G and a radius of gyration about the pivot O of 295 mm. If the pendulum is released from rest when θ = 0, plot the total force supported by the bearing at O along with its normal and tangential components as a function of θ. Let θ range between 0 and 180°. Problem Formulation The free body and mass acceleration diagrams are identical to those in the sample problem. The main difference in approach is that we will obtain results at an arbitrary angle θ rather than at 60°. [ΣM O = I Oα ] mgr cos θ = mk 02α [ωdω = αdθ ] ∫ ω 0 ωdω = ω2 = [ΣF n = mr ω 2 [ΣFt = mr α ] ] 2 gr k 02 gr k 02 α= gr k 02 θ gr 0 2 0 ∫ cos θdθ = k cos θ sin θ sin θ O n − mg sin θ = mr ω 2 − Ot + mg cosθ = mr α After substituting for α and ω we have,  r2 O n = mg 1 + 2 2  k0    sin θ   The magnitude of the force at O is, O= (On ) 2 + (Ot )2  r2 Ot = mg 1 − 2  k 0    cos θ   114 CH. 6 PLANE KINETICS OF RIGID BODIES After substituting r = 0.25 m, k 0 = 0.295 m, m = 7.5 kg, and g = 9.81 m/s2 all forces will be functions of θ only. MATLAB Script %%%%%%%%%%%%%%%% Script %%%%%%%%%%%%%%%%%%% % This script plots the force at O and its % components as functions of theta rb = 0.25; k0 = 0.295; mg = 7.5*9.81; theta = 0:0.01:pi; On = mg*(1+2*(rb/k0)^2)*sin(theta); Ot = mg*(1-(rb/k0)^2)*cos(theta); O=sqrt(On.^2+Ot.^2); td=180*theta/pi; % converts to degrees plot(td,O,td,On,td,Ot) xlabel('theta (degs)') title('Force at O (Newtons)') PLANE KINETICS OF RIGID BODIES 115 6.3 Problem 6/98 (General Plane Motion) The slender rod of mass m and length l is released from rest in the vertical position with the small roller at end A resting on the incline. (a) Determine the initial acceleration of A (aA) and plot aA versus θ for 0 ≤ θ ≤ 90°. (b) Determine the maximum value of aA over this range and the angle θ at which it occurs. Problem Formulation The free-body and mass acceleration diagrams are shown to the right. Shown on the mass acceleration diagram are the two components of the acceleration of the center of mass G obtained from the following kinematic relation, aG = aA + (aG/A)n + (aG/A)t where (aG/A)n = [ΣM A l 2 l ω =0, (aG/A)t = α . 2 2 = I α + ma d [ΣFx = ma x ] ] 0= l l l 1 ml 2α + m α − ma A cos θ 12 2 2 2 l   mg sin θ = m a A − α cos θ  2   These two equations can be solved simultaneously to give, α= 6( g / l ) sin θ cos θ 4 − 3 cos 2 θ aA = 4 g sin θ 4 − 3 cos 2 θ At first, the answer to part (b) seems obvious. Intuitively, we would like to say that the maximum acceleration is aA = g and occurs at θ = 90°. But this intuition neglects the effects of the bar’s rotation upon the acceleration. As we will see below, the maximum acceleration is somewhat larger than g. The maximum acceleration is obtained in the usual manner. The orientation where the maximum occurs is first found by solving the equation da A / dθ = 0 116 CH. 6 PLANE KINETICS OF RIGID BODIES for θ. This angle is then substituted back into the expression for aA to yield the maximum acceleration. MATLAB Worksheet and Script When using the solve command to solve simultaneous equations, be sure that the variables solved for are single characters. Thus, in the following we let x = aA and y = α. Also remember that the equations solved must be in the form where the right hand side is zero. The equals sign is omitted. EDU» syms x y L theta g m EDU» eqn1 = 1/12*m*L^2*y+m*L/2*y*L/2-m*x*L/2*cos(theta) eqn1 = 1/3*m*L^2*y-1/2*m*x*L*cos(theta) EDU» eqn2 = m*g*sin(theta)-m*(x-L/2*y*cos(theta)) eqn2 = m*g*sin(theta)-m*(x-1/2*L*y*cos(theta)) EDU» [x,y]=solve(eqn1,eqn2) x = -4*g*sin(theta)/(-4+3*cos(theta)^2) y = -6*cos(theta)*g*sin(theta)/L/(-4+3*cos(theta)^2) The result for x will be used in the script below to plot aA versus θ. It is convenient to go ahead and do part (b) first as we have everything already set up. EDU» aA = subs(x,g,9.81) aA = -981/25*sin(theta)/(-4+3*cos(theta)^2) EDU» daA = diff(aA,theta) daA = -981/25*cos(theta)/(-4+3*cos(theta)^2)-5886/25*sin(theta)^2/ (-4+3*cos(theta)^2)^2*cos(theta) EDU» solve(daA,theta) ans = 1/2*pi] atan(1/6*3^(1/2)*6^(1/2)) -atan(1/6*3^(1/2)*6^(1/2)) -atan(1/6*3^(1/2)*6^(1/2))+pi atan(1/6*3^(1/2)*6^(1/2))-pi MATLAB has found five solutions. Now we use the eval command to put these in a form easier to understand. PLANE KINETICS OF RIGID BODIES 117 EDU» eval(ans) ans = 1.5708 0.6155 -0.6155 2.5261 -2.5261 Only two of the five solutions are in the range from 0 to 90°. These two solutions are π/2 (90°) and 0.6155 rad (35.3°). Substitution will reveal which is the maximum. Of course, we could also look at the plot below to see that the second solution corresponds to a maximum. EDU» subs(aA,theta,.6155) ans = 11.3276 EDU» subs(aA,theta,pi/2) ans = 9.8100 % this result shouldn’t be surprising Thus, (aA)max = 11.33 m/s2 when θ = 35.3°. %%%%%%%%%%%%%%%%%% Script %%%%%%%%%%%%%%%%% % This script plots aA versus theta g = 9.81; theta = 0:0.01:pi/2; aA = -4*g*sin(theta)./(-4+3*cos(theta).^2); plot(theta*180/pi, aA) xlabel('theta (deg)') ylabel('a_A (m/s^2)') %%%%%%%%%%%%%%%% end of script %%%%%%%%%%%% 118 CH. 6 PLANE KINETICS OF RIGID BODIES 6.4 Problem 6/104 (General Plane Motion) The uniform 12-ft pole is hinged to the truck bed and released from the vertical position as the truck starts from rest with an acceleration a. If the acceleration remains constant during the motion of the pole, derive an expression for the angular velocity ω in terms of a, g, and L where L is the length of the pole. Plot ω versus θ for a = 3 ft/s2. Problem Formulation The free-body and mass acceleration diagrams are shown to the right. Shown on the mass acceleration diagram are the three components of the acceleration of the center of mass G obtained from the following kinematic relation, aG = aO + (aG/O)n + (aG/O)t where aO = a, (aG/O)n = r ω 2 , (aG/O)t = r α , and r = L / 2 = 6 feet. [ΣM O = I α + ma d ] mg L L L L 1 sin θ = mL2α + m α − ma cos θ 12 2 2 2 2 α= 3 (g sin θ + a cosθ ) 2L Now we integrate the relation ωdω = αdθ to obtain, 1 2 3 ω = 2 2L ω= θ ∫ (g sin θ + a cosθ )dθ = 2L (g (1 − cosθ ) + a sin θ ) 0 3 (g (1 − cosθ ) + a sin θ ) L 3 PLANE KINETICS OF RIGID BODIES 119 MATLAB Script %%%%%%%%%%%%%%%%%%%% Script %%%%%%%%%%%%%%%%%%%% L = 12; g = 32.2; a = 3; theta = 0:0.01:pi/2; omega = sqrt(3/L*(g*(1-cos(theta))+a*sin(theta))); plot(theta*180/pi,omega) xlabel('theta (deg)') ylabel('omega (rad/s)') %%%%%%%%%%%%%%%%%end of script %%%%%%%%%%%%%%%%% 120 CH. 6 PLANE KINETICS OF RIGID BODIES 6.5 Sample Problem 6/10 (Work and Energy) The 4-ft slender bar weighs 40 lb with a mass center at B and is released from rest in the position for which θ is essentially zero. Point B is confined to move in the smooth vertical guide, while end A moves in the smooth horizontal guide and compresses the spring as the bar falls. Plot the angular velocity of the bar and the velocities of A and B as a function of θ from 0 to 90°. The stiffness of the spring is 30 lb/in. Problem Formulation From the figure to the right, the lengths CB and CA are 2sinθ and 2cosθ respectively. Using the instantaneous center C we can write the two velocities in terms of the angular velocity ω. v A = CAω = 2ω cos θ v B = CBω = 2ω sin θ Now we need to divide the range for θ into two distinct intervals depending upon whether or not the spring has been engaged. Since the velocities for A and B are known in terms of ω and θ, we need to find only the angular velocity in these two intervals. From the diagram we see that A will first contact the spring at an angle θ = sin −1 (18 / 24) = 0.8481 rads (48.6°). (a) Before the spring is engaged (θ ≤ 48.6°). [T = 1 1 mv 2 + I ω 2 ] 2 2 ∆T = ( 1 40 (2ω cosθ )2 + 1  1 40 4 2 ω 2 2  12 32.2  2 32.2 ) ∆T = 0.8282 4 − 3 cos 2 θ ω 2 [ ∆V g = W∆h] ∆Vg = 40(2 cosθ − 2) = 80(cosθ − 1) We now substitute into the energy equation U 1′− 2 = ∆T + ∆V g = 0 , ( ) 0 = 0.8282 4 − 3 cos 2 θ ω 2 + 80(cos θ − 1) , from which we find PLANE KINETICS OF RIGID BODIES 121 ω = 9.829 1 − cos θ 4 − 3 cos 2 θ (b) After the spring is engaged (48.6 ≤ θ ≤ 90°). The kinetic and potential energies are the same as in part (a). At any angle θ, point A has moved 2sinθ feet to the left. Thus, the spring is compressed by 2sinθ − 18/12 feet. 1 [V e = kx 2 ] 2 2 18  1  lb  in  ∆Ve =  30 12  2 sin θ −  − 0 12  2  in  ft  3  ∆Ve = 180 2 sin θ −  2  2 Again, we substitute into the energy equation U 1′− 2 = ∆T + ∆V g + ∆Ve = 0 , ( ) 3  0 = 0.8282 4 − 3 cos 2 θ ω 2 + 80(cos θ − 1) + 180 2 sin θ −  2  ω = 2.457 2 216 sin θ − 16 cos θ − 144 sin 2 θ − 65 4 − 3 cos 2 θ MATLAB Scripts In the following, terms ending with “_a” are for part (a) where θ is between 0 and 48.6° while terms ending with “_b” are for part (b) where θ is greater than 48.6°. %%%%%%%%%%%%%%%% Script #1 %%%%%%%%%%%%%%%%% % This script solves the energy equation for omega % for cases a and b syms theta omega DT = 0.8282*(4-3*cos(theta)^2)*omega^2; DVg = 80*(cos(theta)-1); DVe = 180*(2*sin(theta)-3/2)^2; U12_a = DT + DVg U12_b = DT + DVg + DVe omega_a = solve(U12_a,omega) omega_b = solve(U12_b,omega) %%%%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%% 122 CH. 6 PLANE KINETICS OF RIGID BODIES Output of script #1 U12_a = (4141/1250-12423/5000*cos(theta)^2)*omega^2+80*cos(theta)-80 U12_b = (4141/1250-12423/5000*cos(theta)^2)*omega^2+80*cos(theta)80+180*(2*sin(theta)-3/2)^2 omega_a = 200/(-16564+12423*cos(theta)^2)*41410^(1/2)*((-4+3*cos(theta)^2)* (cos(theta)-1))^(1/2) -200/(-16564+12423*cos(theta)^2)*41410^(1/2)*((-4+3*cos(theta)^2)* (cos(theta)-1))^(1/2) It shouldn’t be surprising that MATLAB has found two solutions to the equation. As you can see, they have the same magnitude with one being positive and the other negative. The solution of interest is the positive one which turns out to be the second solution. This solution is copied and pasted into script #2. omega_b = 50/(4141+12423*sin(theta)^2)*(-(41410+124230*sin(theta)^2)* (16*cos(theta)+65+144*sin(theta)^2-216*sin(theta)))^(1/2) -50/(4141+12423*sin(theta)^2)*(-(41410+124230*sin(theta)^2)* (16* cos(theta)+65+144*sin(theta)^2-216*sin(theta)))^(1/2) In this case, the first solution is positive. %%%%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%% % This script plots omega versus theta. Note that % two range variables are used to divide the range % for theta into the two regions. % a = theta for theta < 48.6 deg. (.8481 rads) % b = theta for theta > 48.6 deg. % The results for omega are copied and pasted from % the output of script #1 a = 0:0.001:0.8481; b = 0.8481:0.001:pi/2; omega_a = -200./(-16564+12423*cos(a).^2)*41410^(1/2).*((4+3*cos(a).^2).*(cos(a)-1)).^(1/2); omega_b = 50./(4141+12423*sin(b).^2).*((41410+124230*sin(b).^2).*(16*cos(b)+65+144*sin(b).^2216*sin(b))).^(1/2); plot(a*180/pi, omega_a, b*180/pi, omega_b) PLANE KINETICS OF RIGID BODIES 123 xlabel('theta (degrees)') ylabel('omega (rad/sec)') %%%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% Script #3 %%%%%%%%%%%%%%%%% % This script plots vA and vB versus theta. a = 0:0.001:0.8481; b = 0.8481:0.001:pi/2; omega_a = -200./(-16564+12423*cos(a).^2)*41410^(1/2).*((4+3*cos(a).^2).*(cos(a)-1)).^(1/2); omega_b = 50./(4141+12423*sin(b).^2).*((41410+124230*sin(b).^2).*(16*cos(b)+65+144*sin(b).^2216*sin(b))).^(1/2); vA_a = 2*omega_a.*cos(a); vA_b = 2*omega_b.*cos(b); vB_a = 2*omega_a.*sin(a); vB_b = 2*omega_b.*sin(b); plot(a*180/pi, vA_a, b*180/pi, vA_b, a*180/pi, vB_a, b*180/pi, vB_b) xlabel('theta (degrees)') ylabel('velocity (ft/sec)') %%%%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%% 124 CH. 6 PLANE KINETICS OF RIGID BODIES For Further Study A striking feature of the velocity curves above is the sudden change in shape when the spring is engaged at about θ = 49°. The details depend very much upon the relative magnitudes of the weight and spring constant. To illustrate, the figure below shows the angular velocity for several different values of the spring constant k. Note that for stiff springs, the angular velocity goes to zero before reaching θ = 90°. The physical explanation for this is that, for a stiff spring, the bar will rebound before it reaches the horizontal position. PLANE KINETICS OF RIGID BODIES 125 6.6 Problem 6/206 (Impulse/Momentum) Determine the minimum velocity v that the wheel may have and just roll over the obstruction. The centroidal radius of gyration of the wheel is k, and it is assumed that the wheel does not slip. Plot v versus h for three cases: k = ½, ¾, and 1 m. For each case take r = 1 m. Problem Formulation During Impact: Conservation of Angular Momentum As usual, we neglect the angular impulse of the weight during the short interval of impact. With this assumption we have conservation of angular momentum about point A. Immediately before impact, the center of the wheel is not moving in a circular path about A and we need to use the formula for general plane motion. Note that I = mk 2 . H A = I ω + mv d = mk 2ω + mv(r − h ) = mk 2 v + mv(r − h ) r We will use primes to denote the state immediately after impact. Since the wheel now rotates about A we can use the simpler formula H A = I Aω . Note that, by the ( ) parallel axis theorem, I A = I + mr 2 = m k 2 + r 2 . ( ) v′ ′ H A = I Aω ′ = k 2 + r 2 r Setting H A = H A′ and solving yields, rh   v ′ = v 1 − 2  2  k +r  After Impact: Work-Energy ∆T + ∆V g = 0 = ( ) 1 I A 0 2 − ω ′ 2 + mgh 2 126 CH. 6 PLANE KINETICS OF RIGID BODIES ( ) 1  v′  m k 2 + r 2   = mgh 2 r 2 Substituting the result for v′ into the above equation followed by simplification yields, v= ( r 2 gh k 2 + r 2 ) k + r − rh 2 2 MATLAB Scripts %%%%%%%%%%%%%%% Script #1 %%%%%%%%%% % symbolic solution for v syms m g k r h v vp Ib = m*k^2; IA = m*(k^2+r^2); HA = Ib*v/r+m*v*(r-h); HAp = IA*vp/r; vp = solve(HA-HAp,vp) % conservation of angular momentum eqn = 1/2*IA*vp^2/r^2-m*g*h; % work/energy solve(eqn,v) %%%%%%%%%%% end of script %%%%%%%%%% Output of script #1 vp = v*(k^2+r^2-r*h)/(k^2+r^2) ans = (2*g*h*r^2+2*g*h*k^2)^(1/2)*r/(k^2+r^2-r*h) -(2*g*h*r^2+2*g*h*k^2)^(1/2)*r/(k^2+r^2-r*h) The first solution is the one we want. %%%%%%%%%%%%%%% Script #2 %%%%%%%%%% % This script plots v versus h for the specified values of k. % The result for v below was obtained by substituting g = 9.81 % and r = 1 in to the symbolic result found in script #1 v = inline('4.429*sqrt(h+h*k^2)/(1+k^2-h)'); v = vectorize(v); % puts “.”’s in expression h = 0:0.01:1; plot(h,v(h,1/2),h,v(h,3/4),h,v(h,1)) xlabel('h (m)') ylabel('v (m/s)') %%%%%%%%%%% end of script %%%%%%%%%% PLANE KINETICS OF RIGID BODIES 127 INTRODUCTION TO THREEDIMENSIONAL DYNAMICS OF RIGID BODIES This chapter presents a brief introduction to rigid body dynamics in three dimensions. In problem 7.1, the general 3-D motion of three connected bars is investigated. In particular, the angular velocities of two of the bars are plotted versus the length of the third bar. solve is used to solve four equations symbolically for four unknowns. In problem 7.2 we consider a bent plate rotating about a fixed axis. The problem illustrates a simplified version of what engineers might do in a real design situation. Two dimensions of the bent plate are left as variables and the objective of the problem is to find all suitable values of those dimensions which satisfy several constraints simultaneously. Here we illustrate a graphical approach to this type of design problem. 7 130 CH. 7 INTRODUCTION TO 3-D DYNAMICS OF RIGID BODIES 7.1 Sample Problem 7/3 (General Motion) Crank CB rotates about the horizontal axis with an angular velocity ω1 = 6 rad/s, which is constant for a short interval of motion that includes the position shown. Link AB has a balland-socket fitting on each end and connects crank DA with CB. Let the length of crank CB be d mm (instead of 100 mm as in the sample problem in your text) and plot ω2 and ωn as a function of d for 0 ≤ d ≤ 200 mm. ω2 is the angular velocity of crank DA while ωn is the angular velocity of link AB. Problem Formulation Our analysis will follow closely that in the sample problem in your text. vA = vB + ωn×rA/B where vA = 50ω2j vB = 6di rA/B = 50i + 100j + dk Substitution into the velocity equation gives i 50ω 2 j = 6di + ω nx 50 j k ω ny ω nz 100 d Expanding the determinant and equating the i, j, and k components yields the following three equations ( ) d 6 + ω ny − 100ω nz = 0 50(ω 2 − ω nz ) + dω nx = 0 2ω nx − ω ny = 0 At this point we have three equations with four unknowns. As explained in the sample problem in your text, the fourth equation comes by requiring ωn to be normal to vA/B ωn• rA/B = 50ωnx + 100ωny + dωnz = 0 These four equations will be solved simultaneously for ω2, ωnx, ωny, and ωnz. Once this is done, 3-D DYNAMICS OF RIGID BODIES 131 2 2 ω n = ω nx2 + ω ny + ω nz MATLAB Scripts %%%%%%%%%%%%%%%% Script #1 %%%%%%%%%%%%%%%%%% % This script solves four equations for % omega_2 (w) and the three components of % omega_n (x, y, and z). The magnitude of % omega_n is then found from its components. syms w x y z d eqn1 = d*(6+y)-100*z; eqn2 = 50*(w-z)+d*x; eqn3 = 2*x-y; eqn4 = 50*x+100*y+d*z; [w,x,y,z]=solve(eqn1,eqn2,eqn3,eqn4) omega_n = sqrt(x^2+y^2+z^2); omega_n = simplify(omega_n) %%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%% Output of script #1 w = 3/50*d (ω2) x = -3*d^2/(d^2+12500) y = -6*d^2/(d^2+12500) z = 750*d/(d^2+12500) omega_n = 3*5^(1/2)*(d^2/(d^2+12500))^(1/2) (ωn) %%%%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%%% % This script plots omega_2 and the magnitude of % omega_n as functions of d d = 0:0.5:200; omega_2 = 3/50*d; omega_n = 3*5^(1/2)*(d.^2./(d.^2+12500)).^(1/2); plot(d, omega_2, d, omega_n) xlabel('d (mm)') title('angular velocity (rad/s)') %%%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%% 132 CH. 7 INTRODUCTION TO 3-D DYNAMICS OF RIGID BODIES 7.2 Sample Problem 7/6 (Kinetic Energy) The bent plate has a mass of 70 kg per square meter of surface area and revolves around the zaxis at the rate ω = 30 rad/s. Let the dimensions of part B be a and b where a is the dimension parallel to the x-axis and b is the dimension parallel to the z-axis. Part A remains unchanged. (a) Find all suitable values for a and b which satisfy the following conditions: a ≤ 0.2 m, b ≤ 0.6 m, and 15 ≤ T ≤ 30 J where T is the kinetic energy of the plate. (b) Find a and b for the case where T = 40 J and H0 = 5 N•m•s where H0 is the magnitude of the angular momentum about O. Problem Formulation Substitution of ωx = 0, ωy = 0, and ωz = ω into equations 7/11 and 7/18 of your text yields. HO = ω(−Ixzi −Iyzj +Izzk) 2 2 2 2 H O = H Ox + H Oy + H Oz = ω I xz2 + I yz + I zz2 3-D DYNAMICS OF RIGID BODIES 133 T= 1 I zz ω 2 2 The moments and products of inertia for part A remain unchanged. For part B, mB = 70ab where a and b are in meters. The moments and products of inertia for part B are I zz = I zz + md 2 = 2  mB 2 a  2 a + m B  (.125) +     12  2     a  b  I xz = I xz + md x d z = 0 + m B     2  2  b I yz = I yz + md y d z = 0 + m B (0.125)  2 The total moment and products of inertia are found by adding the above to those found for part A (see the sample problem in your text). After substituting for mB and simplifying we have, I zz = 0.00456 + 1.094ab + 23.33a 3b I xz = 17.5a 2 b 2 I yz = 0.00273 + 4.375ab 2 Substitution of these results will give HO and T as functions of a and b. Part (a) The most efficient way to show the acceptable ranges for a and b is to find the required relationship between these two dimensions in order to satisfy the upper and lower bounds on T. This is accomplished by substituting these bounds for T in the equation above and then solving that equation for b as a function of a. To illustrate, consider the lower limit on T (15 J). Substituting T = 15 into the equation above gives 1 2 I zz (30 ) = 2.052 + 492.2ab + 10,500a3b 2 0.078921 b= (for T = 15 J) a 3 + 64a 2 T = 15 = Solving for b, ( ) A similar result can be obtained for the upper limit, 134 CH. 7 INTRODUCTION TO 3-D DYNAMICS OF RIGID BODIES b= 0.17035 a 3 + 64a 2 ( ) (for T = 30 J) Plotting these two functions defines the acceptable regions for a and b. Part (b) This is similar to (a) except that we solve two equations (T = 40 and H0 = 5) simultaneously for two unknowns, a and b. The result is a = 0.1211 m and b = 0.4852 m. MATLAB Scripts %%%%%%%%%%%%%%% Script #1 %%%%%%%%%%%%%%%%%%%%%% % The first part of this script solves the equations % T = 15 and T = 30 symbolically for b in terms of a. % These two results are named b15 and b30 respectively. % When using the symbolic solve it is best to rewrite % the equations in the form expression = 0. The "=0" % is automatically understood by MATLAB and doesn't % have to be typed in. In the following, our two % expressions are named eqn1 and eqn2. syms a b mB = 70*a*b; omega = 30; Izz = mB/12*a^2+mB*(0.125^2+(a/2)^2)+0.00456; Ixz = mB*a/2*b/2; Iyz = mB*0.125*b/2+0.00273; H0 = omega*sqrt(Ixz^2+Iyz^2+Izz^2) T = 1/2*Izz*omega^2 eqn1 = T - 15; eqn2 = T - 30; b15 = solve(eqn1,b) b30 = solve(eqn2,b) % This part of the script solves the two equations % for part b. eqn3 = T - 40; eqn4 = H0 - 5; [a,b]=solve(eqn3,eqn4) %%%%%%%%%%% end of script %%%%%%%%%%%%%%%%%%%%%% Output of script #1 H0 = 1/10000*(27562500000000*a^4*b^4+1722656250000*a^2*b^4+2149875000*a*b^2+25 42185+49000000000000*a^6*b^2+4593750000000*a^4*b^2+19152000000*a^3*b+107 666015625*a^2*b^2+897750000*a*b)^(1/2) 3-D DYNAMICS OF RIGID BODIES 135 T = 2625*a^3*b+31500*a*b*(1/64+1/4*a^2)+513/250 b15 = 8632/109375/a/(64*a^2+3) b30 = 18632/109375/a/(64*a^2+3) a= -.118672955347671873813795187 -.924961825359232183759243317e-1-.20405912363720769438094726*i -.924961825359232183759243317e-1+.20405912363720769438094726*i -.257514497287494286960690298e-2-.25491731461066579787563883*i -.257514497287494286960690298e-2+.25491731461066579787563883*i .349640543257908784471004656e-2-.265040090095767094717337313*i .349640543257908784471004656e-2+.26504009009576709471733731*i .903798369106820783468910094e-1-.204312206452426767071840559*i .9037983691068207834689100945e-1+.2043122064524267670718405*i .121063125678745863921655544 b= -.499591682557401531621796988 .286519678237340239067367533+.281089760611403745879347469*i .286519678237340239067367533-.281089760611403745879347469*i .64375775885513160849563080e-1-.778487038145858201006376746*i .64375775885513160849563080e-1+.77848703814585820100637674*i -.5367082781855869430616764e-1-.579398660286892007965618348*i -.5367082781855869430616764e-1+.57939866028689200796561834*i -.29810587868680635995547322+.28455901589818267092944534*i -.29810587868680635995547322-.284559015898182670929445341*i .485167563521434010754982656 We see from the above that MATLAB has found ten solutions to our equations. Of these only the first and last are real. The first solution has negative values for a and b and thus can be excluded. This leaves only the last solution, a = 0.1211 m; b = 0.4852 m %%%%%%%%%%%%%%% Script #2 %%%%%%%%%%%%%%%%%%%%%% % This script plots the two expressions for b % obtained with script #1 a = 0.01:0.001:0.2; b15 = 8632/109375./a./(64*a.^2+3); b30 = 18632/109375./a./(64*a.^2+3); plot(a, b15, a, b30) xlabel('a (m)') 136 CH. 7 INTRODUCTION TO 3-D DYNAMICS OF RIGID BODIES ylabel('b (m)') axis([0 0.2 0 0.6]) %%%%%%%%%%% end of script %%%%%%%%%%%%%%%%%%%%%% The two curves above represent the values of a and b for which T is exactly 15 or 30 J. Thus, the acceptable values of a and b satisfying the condition 15 ≤ T ≤ 30 J are all those combinations lying on or between the two curves. VIBRATION AND TIME RESPONSE This chapter considers an important class of Dynamics problems that involve linear or angular oscillations of a body or structure about some equilibrium position or configuration. Very few of the homework problems in your text require you to actually plot the oscillations of a body versus time. For this reason, all of the problems in this chapter will involve such plots. This is very useful in visualizing the time response of a vibrating system, especially for the case of damped or forced vibrations. Problem 8.1 looks at the effects of damping coefficient upon time response while Problem 8.3 considers the effects of initial conditions. 8 138 CH. 8 VIBRATION AND TIME RESPONSE 8.1 Sample Problem 8/2 (Free Vibration of Particles) The 8-kg body is moved 0.2 m to the right of the equilibrium position and released from rest at time t = 0. Plot the displacement as a function of time for three cases, c = 8, 32, and 56 N•s/m. The spring stiffness k is 32 N/m. Problem Formulation As in the sample problem, the natural circular frequency is ω n = k / m = 32 / 8 = 2 rad/sec. Now we find the damping ratio for each case (from ζ = c/2mωn) with the result ζ = 0.25, 1, and 1.75 for c = 8, 32, and 56 N•s/m respectively. (a) ζ = 0.25. Since ζ < 1, the system is underdamped. The damped natural frequency is ω d = ω n 1 − ζ 2 = 1.937 rad/sec. The displacement and velocity are x = Ce −ζω n t sin(ω d t + ψ ) = Ce −t / 2 sin(1.937t + ψ ) x& = −0.5Ce− t / 2 sin(1.937t + ψ ) + 1.937Ce −t / 2 cos(1.937t + ψ ) From the initial conditions x 0 = 0.2 and x& 0 = 0 we find C = 0.207 m and ψ = 1.318 rad. (b) ζ = 1. For ζ = 1, the system is critically damped. The displacement and velocity are x = ( A1 + A2 t )e −ω n t = ( A1 + A2 t )e −2t x& = A2 e −2t − 2( A1 + A2 t )e −2t From the initial conditions x 0 = 0.2 and x& 0 = 0 we find A1 = 0.2 m and A2 = 0.4 m/s. (c) ζ = 1.75. Since ζ > 1, the system is overdamped. The displacement and velocity are  −ζ + ζ 2 −1 ω t  n x = B1 e   −ζ − ζ 2 −1 ω t  n + B2 e  = B1 e −0.628t + B 2 e − 6.372 t VIBRATION AND TIME RESPONSE 139 x& = −0.628B1 e −0.628 t − 6.372 B 2 e −6.372t From the initial conditions x 0 = 0.2 and x& 0 = 0 we find B1 = 0.222 m and B2 = −0.0219 m/s. MATLAB Script %%%%%%%%%%%%%%% Script %%%%%%%%%%%%%%%%%%%%%% t=0:0.01:5; C = 0.207; psi = 1.318; xa = C*sin(1.937*t+psi).*exp(-t/2); A1 = 0.2; A2 = 0.4; xb = (A1+A2*t).*exp(-2*t); B1 = 0.222; B2 = -0.0219; xc = B1*exp(-0.628*t)+B2*exp(-6.372*t); x0 = 0.*t; plot(t,xa,t,xb,t,xc,t,x0) xlabel('time (sec)') ylabel('x (m)') %%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%%% 140 CH. 8 VIBRATION AND TIME RESPONSE 8.2 Problem 8/139 (Damped Free Vibration) The mass of a critically damped system having a natural frequency ω n is released from rest at an initial displacement x0. (a) Determine the time t required for the mass to reach the position x = 0.1x0 if ω n = 4 rad/s. (b) Plot the nondimensional displacement x/x0 for ω n = 2, 4, and 8 rad/s. Problem Formulation Start with the equation for the critically damped case on page 606 of your text. x = ( A1 + A2 t )e −ωnt First we have to evaluate the constants in terms of x0 using the initial conditions. This will require the derivative of x. x& = A2 e −ωnt − ( A1 + A2 t )ω n e −ωnt Now, substituting the initial conditions, x(t = 0 ) = x0 = A1 x& (t = 0) = 0 = A2 − A1ω n These two equations give A1 = x 0 and A2 = x0ω n x = ( x0 + x 0ω n t )e −ωnt or η = (1 + ω n t )e −ωnt where η is the non-dimensional displacement x/x0. Part (a) will require the solution to the equation VIBRATION AND TIME RESPONSE 141 0.1 = (1 + 4t )e −4t For part (b) we will simply plot η for three different natural frequencies ω n . MATLAB Worksheet and Script Part (a) EDU>> f = '(1 + 4*t)*exp(-4*t) = 0.1' f= (1 + 4*t)*exp(-4*t) = 0.1 EDU>> solve(f,'t') ans = -.24044468956330014201486061769676 .97243004246685726447599505623170 Thus, the time t required for the mass to reach the position x = 0.1x0 when ω n = 4 rad/s is t = 0.972 sec. Part (b) %========= script for part b ======== eta = inline('(1 + wn*t).*exp(-wn*t)') t = 0:0.05:3; plot(t,eta(t,2),t,eta(t,4),t,eta(t,8)) xlabel('time (sec)') title('non-dimensional displacement') %========= end of script ============ eta = Inline function: eta(t,wn) = (1 + wn*t).*exp(-wn*t) 142 CH. 8 VIBRATION AND TIME RESPONSE VIBRATION AND TIME RESPONSE 143 8.3 Sample Problem 8/6 (Forced Vibration of Particles) The 100-lb piston is supported by a spring of modulus k = 200 lb/in. A dashpot of damping coefficient c = 85 lb-sec/ft acts in parallel with the spring. A fluctuating pressure p = 0.625 sin(30t) (psi) acts on the piston, whose top surface area is 80 in2. Plot the response of the system for initial conditions x0 = 0.05 ft and x& 0 = 5, 0, and −5 ft/sec. Problem Formulation The particular (steady state) solution was found in the sample problem in your text, x p = X sin(ωt − φ ) where X = 0.01938 m, φ = 1.724 rad and ω = 30 rad/sec. Also from the sample problem, ω n = k / m = 27.8 rad/sec and ζ = c/2mωn = 0.492. The complete solution is found by adding the complementary (transient) and particular solutions. Since the system is underdamped (ζ < 1), the complementary solution is, x c = Ce −ζω n t sin(ω d t + ψ ) where ω d = ω n 1 − ζ 2 = 24.2 rad/sec. The displacement of the system is x = x c + x p = Ce −ζω n t sin(ω d t + ψ ) + X sin (ωt − φ ) x = Ce −13.68t sin( 24.2t + ψ ) + 0.01938 sin(30t − 1.724) The velocity is found by differentiating x, x& = Cω d e −ζω n t cos(ω d t + ψ ) − Cζω n e −ζω n t sin(ω d t + ψ ) + Xω cos(ωt − φ ) x& = Ce −13.68t (24.2 cos( 24.2t + ψ ) − 13.68 sin(24.2t + ψ ) ) + 0.5814 cos(30t − 1.724) The constants C and ψ are found from the initial conditions, 144 CH. 8 VIBRATION AND TIME RESPONSE x 0 = 0.05 = C sin ψ + X sin φ = C sinψ − 0.0192 x& 0 = Cω d cos ψ − Cζω n sin ψ + Xω cos φ = 24.2C cos ψ − 13.7C sin ψ − 0.0887 The first equation can be solved for C = 0.0692/sinψ . Substitution into the second equation gives ψ .  1.675     x& 0 + 1.035  ψ = tan −1  C= 0.0692 sin ψ This yields the following values for C and ψ . x0 = 0.05 ft, x& 0 = 5 ft/s, C = 0.259 ft, ψ = 0.271 rad x0 = 0.05 ft, x& 0 = 0 ft/s, C = 0.081 ft, ψ = 1.017 rad x0 = 0.05 ft, x& 0 = −5 ft/s, C = −0.178 ft, ψ = −0.400 rad MATLAB Script %%%%%%%%%%%%%%%% Script %%%%%%%%%%%%%%%%%%%%% % This script plots x versus t for the three % sets of initial conditions x = inline('C*sin(24.2*t+psi)*exp(-13.68*t)+.01938*sin(30*t-1.724)') x = vectorize(x); t = 0:0.001:0.5; plot(t,x(.259,.271,t),t,x(.081,1.017,t),t,x(-.178,-.4,t)) xlabel('t (sec)') ylabel('x (ft)') %%%%%%%%%%%%% end of script %%%%%%%%%%%%%%%%%% VIBRATION AND TIME RESPONSE 145 Notice how quickly the three cases converge to the steady state solution.