MATLAB Script for Vector Calculations and Equilibrium, Cheat Sheet of Mathematics

Download MATLAB Script for Vector Calculations and Equilibrium and more Cheat Sheet Mathematics in PDF only on Docsity! Using MATLAB for Statics and Dynamics by Ron Larsen and Steve Hunt 1. Resolving Forces, Calculating Resultants 2. Dot Products 3. Equilibrium of a Particle, Free-Body Diagrams 4. Cross Products and Moments of Forces 5. Moment of a Couple 6. Reduction of a Simple Distributed Loading 7. Equilibrium of a Rigid Body 8. Dry Friction 9. Finding the Centroid of Volume 10. Resultant of a Generalized Distributed Loading 11. Calculating Moments of Inertia 12. Curve-Fitting to Relate s-t, v-t, and a-t Graphs 13. Curvilinear Motion: Rectangular Components 14. Curvilinear Motion: Motion of a Projectile 15. Curvilinear Motion: Normal and Tangential Components 16. Dependent Motion of Two Particles 17. Kinetics of a Particle: Force and Acceleration 18. Equations of Motion: Normal and Tangential Components 19. Principle of Work and Energy 20. Rotation About a Fixed Axis 2.4 - 2.6 2.3 - 2.5 Hibbeler Bedford and Fowler 2.9 2.5 3.3 3.2 - 3.3 9.2 7.4 10.1 - 10.2 8.1 - 8.2 4.2 - 4.3 2.6, 4.3 4.6 4.4 4.10 7.3 St at ic s D yn am ic s 5.1 - 5.3 5.1 - 5.2 8.2 9.1 9.5 7.4 12.6 2.3 12.3 2.2 12.5 2.3 12.7 2.3 12.9 ch. 2 13.4 3.1 - 3.4 13.5 3.4 14.3 8.1 - 8.2 16.3 9.1 Resolving Forces, Calculating Resultants Ref: Hibbeler § 2.4-2.6, Bedford & Fowler: Statics § 2.3-2.5 Resolving forces refers to the process of finding two or more forces which, when combined, will produce a force with the same magnitude and direction as the original. The most common use of the process is finding the components of the original force in the Cartesian coordinate directions: x, y, and z. A resultant force is the force (magnitude and direction) obtained when two or more forces are combined (i.e., added as vectors). Breaking down a force into its Cartesian coordinate components (e.g., Fx, Fy) and using Cartesian components to determine the force and direction of a resultant force are common tasks when solving statics problems. These will be demonstrated here using a two-dimensional problem involving co- planar forces. Example: Co-Planar Forces Two boys are playing by pulling on ropes connected to a hook in a rafter. The bigger one pulls on the rope with a force of 270 N (about 60 lbf) at an angle of 55° from horizontal. The smaller boy pulls with a force of 180 N (about 40 lbf) at an angle of 110° from horizontal. a. Which boy is exerting the greatest vertical force (downward) on the hook? b. What is the net force (magnitude and direction) on the hook – that is, calculate the resultant force. 270 N 18 0 N -110° -55° Solution First, consider the 270 N force acting at 55° from horizontal. The x- and y-components of force are indicated schematically, as 1 Note: The angles in this figure have been indicated as coordinate direction angles. That is, each angle has been measured from the positive x axis. 270 N 18 0 N 77° F R FRx F R y 62 N 155 N » F_Rx = F_x1 + ( -F_x2 ) F_Rx = 93.3020 The minus sign was included before Fx2 because it is directed in the –x direction. The result is an x- component of the resultant force of 93 N in the +x direction. Once the x- and y-components of the resultant force have been determined, the magnitude can be calculated using 2 Ry 2 RxR FFF += The MATLAB calculation uses the built-in square-root function sqrt(). » F_R = sqrt( F_Rx ^ 2 + F_Ry ^ 2 ) F_R = 401.3124 The angle of the resultant force can be calculated using any of three functions in MATLAB: Function Argument(s) Notes atan(abs(Fx / Fy)) one argument: abs(Fx / Fy) Returns the included angle atan2(Fy, Fx) two arguments: Fx and Fy Returns the coordinate direction angle Angle value is always between 0 and π radians (0 and 180°) A negative sign on the angle indicates a result in one of the lower quadrants of the Cartesian coordinate system cart2pol (Fx, Fy) two arguments: Fx and Fy Returns the positive angle from the positive x-axis to the vector Angle value always between 0 and 2π radians (0 and 360°) An angle value greater than 180° (π radians) indicates a result in one of the lower quadrants of the Cartesian coordinate system The atan2() function is used here, and FRy is negative because it is acting in the –y direction. » F_Rx = 93.302; » F_Ry = -390.316; » theta = 180/pi * atan2( F_Ry, F_Rx ) theta = -76.5562 Answer, part b) The net force (magnitude and direction) on the hook is now known: FR = 401 N (about 90 lbf) acting at an angle 76.6° below the x-axis. F R FRx F R y 77° θ Annotated MATLAB Script Solution % Determine the x- and y-components of the two forces % (270 N at -55°, and 180 N at -110°) % % Note: These trig. Calculations use the included angles % (55° & 20°), with minus signs added to both y-component % equations to indicate the forces act in the -y direction, % and the F_x2 equation to show that this force acts in % the -x direction. % Calculate the x- and y- components of the first force (270 N) F_x1 = 270 * cos( 55 * pi/180 ); F_y1 = -270 * sin( 55 * pi/180 ); fprintf('\nF_x1 = %8.3f N\t F_y1 = %+9.3f N\n',F_x1,F_y1); % Calculate the x- and y- components of the first force (180 N) F_x2 = -180 * sin( 20 * pi/180 ); F_y2 = -180 * cos( 20 * pi/180 ); fprintf('F_x2 = %7.3f N\t F_y2 = %9.3f N\n',F_x2,F_y2); % Sum the y-components of the two forces to determine the % y-component of the resultant force F_Ry = F_y1 + F_y2; % Sum the x-components of the two forces to determine the % x-component of the resultant force F_Rx = F_x1 + F_x2; fprintf('F_Rx = %7.3f N\t F_Ry = %9.3f N\n\n',F_Rx,F_Ry); % Calculate the magnitude of the resultant force F_R = sqrt( F_Rx ^ 2 + F_Ry ^ 2 ); fprintf('F_R = %8.3f N\n',F_R); % Calculate the angle of the resultant force % (in degrees from the x-axis) theta = atan2( F_Ry, F_Rx ) * 180/pi; fprintf('theta = %7.3f N\n\n',theta); Note that the axes do not intercept at (0, 0, 0). Actually, the middle point on the plot is at (0, 0, 0). The ‘-o’ option in the plot3 function is a line specification that determines line type, and marker symbol of the plotted lines. Here, the ‘-‘ represents a solid line and the ‘o’ is used for circular marker symbols. The usefulness of the 3-d plot for visualization is the ability to rotate the graph to see how the plot looks from different angles. To rotate the plot, use the menu commands: Tools/Rotate 3D, then simply click inside the plot area, hold the left mouse button down, and drag the mouse pointer around the plot. The curve will respond to the location of the pointer. In the graph below, the plot has been rotated until the plane of the vectors is approximately parallel with the page. In this view, the calculated angle of 150° looks to be about right for these vectors. Annotated MATLAB Script Solution % Define the vectors A = [ 1 2 3]; B = [-1 -2 -1]; % Take the dot product dot(A,B) % Check MATLAB's dot product operator by calculating the dot product % explicitly... x = 1; y = 2; z = 3; % coordinate index definitions A(x) * B(x) + A(y) * B(y) + A(z) * B(z) % To calculate the angle between the vectors, first need the magnitude % of each vector. A_mag = sqrt( A(x)^2 + A(y)^2 + A(z)^2 ); B_mag = sqrt( B(x)^2 + B(y)^2 + B(z)^2 ); % Then calculate the angle theta = 180/pi * acos( dot(A,B)/( A_mag * B_mag ) ); % To visualize the vectors, create matrices describing the starting and % ending coordinates of each vector. X = [ 0 0 1 -1]; Y = [ 0 0 2 -2]; Z = [ 0 0 3 -1]; % Then create a 3-d Scatter Plot plot3(X, Y, Z, '-o') % Creates plot grid; % Draws gridlines xlabel('\bf{x-axis}'); % Labels x-axis ylabel('\bf{y-axis}'); % Labels y-axis zlabel('\bf{z-axis}'); % Labels z-axis view(29,20); % 3-D graph viewpoint specification Example: Finding the Components of a Force Vector A force acts on a fixed pin with a magnitude of 270 N (about 60 lbf) at an angle of 55° from horizontal. a. Find the horizontal and vertical components of the force using dot products. F = 27 0 N 55° x y Fx F y Solution To find the magnitude of the x-component of F, calculate the dot product of F and a unit vector in the x direction, ux. x x u )cos(FF xuF • =θ= Note: Most texts do not show ux in the denominator of the last term of this series of equalities, since the magnitude of a unit vector is known to be one. » F = 270; » theta = 55 * pi/180; % 55 deg converted to radians » F_x = F * cos( theta ) F_x = 154.8656 Similarly, the vertical component of F is calculated as: » F_y = F * cos( (90-55) * pi/180 ) F_y = 221.1711 Alternatively, the F and ux vectors can be written in vector form. » F(1) = 270 * cos(55 * pi/180); » F(2) = 270 * cos( (90-55) * pi/180 ); » F(3) = 0; » u_x = [ 1 0 0 ]; And the dot product can be used to calculate the horizontal component of F. » F_x = dot( F, u_x ) F_x = 154.8656 » m = 18; %Mass of light (kg) » g = 9.8; %Gravitational constant (m/s^2) » F_grav = -m * g F_grav = -176.4000 So the downward force exerted by each traffic light is 176.4 N. The minus sign has been included to show that the force is downward, i.e., in the –y direction. If the lights are not moving up or down, the cable must apply an equal but oppositely directed force on the light, since the vertical force components must sum to zero if the body is in equilibrium. Part b. – Free-Body Diagram for Light at B B Fgrav = -176.4 N FAB FBC Since cable section BC is horizontal, there is no vertical component of force in cable section BC. So, all of the weight of the traffic light at B must be carried by cable AB, more specifically, by the vertical component of force in cable AB. » F_ABy = -F_grav F_ABy = 176.4000 The horizontal component of force in cable AB can be determined using the specified angle (cable AB is 5° from horizontal, or 175° from +x). F_ABx = F_ABy / tan( 175 * pi/180) F_ABx = -2016.3000 The horizontal component is 2016 N (about 450 lbf) acting in the –x direction. The force acting in the direction of the cable has a magnitude of 2024 N. » F_AB = sqrt( F_ABx ^ 2 + F_ABy ^ 2 ) F_AB = 2024.0000 Finally, if the light at point B is not moving to the left or right, the sum of the horizontal forces acting on point B must be zero, so the horizontal component of force in section BC is +2016 N. Since there is no vertical component of force in section BC, this is also the total force on section BC at point B. » F_BCx = -F_ABx; % Since the sum of the x-components of force is zero » F_BCy = 0; » F_BC = sqrt( F_BCx ^ 2 + F_BCy ^ 2 ) F_BC = 2016.3000 Part b. – Free-Body Diagram for Light at C CFBC = 2016.3 N FCD Fgrav = -176.4 N The free body diagram for the light at point C is constructed using the following concepts: • If the light at point C is not moving left or right, then the horizontal component of force FCD must be +2016 N. • If the light at point C is not moving up or down, then the vertical component of force FCD must be +176 N. This is essentially the same as the free-body diagram at point B, just flipped left to right. So the magnitude of force FCD is 2024 N, and acts at 5° from horizontal. This can be verified as follows: » F_CDx = -F_BCx F_CDx = 2016.3000 » F_CDy = -F_grav F_CDy = 176.4000 » F_CD = sqrt( F_CDx ^ 2 + F_CDy ^ 2 ) F_CD = 2024.0000 » theta_CD = atan(F_CDy/F_CDx) * 180/pi theta_CD = 5.0000 Note: The value of FBC used in this section is –2016.3 N, as shown in the free-body diagram for point C. See the Annotated MATLAB Script Solution to see how this sign change is handled for the complete problem solution. Part b. – Free-Body Diagram for Light at D FCD = 2024 N D FDE Fgrav = -176.4 N The weight of the traffic light at point D exerts a downward force, Fgrav, of 176.4 N. In addition, force FCD, acting on point D, has a vertical component of 176.4 N, also acting downward. If the light at point D is not moving up or down, then these downward forces must be counterbalanced by the vertical component of force FDE. » F_DEy = -( F_grav + F_CDy ) F_DEy = 352.8000 The horizontal component of force FDE must be equal to –FCD_x if the light at point D is to be stationary. F_DEx = 2016.3000 Once the horizontal and vertical components are known, the angle (in degrees) of cable DE can be determined. » theta_DE = atan(F_DEy/F_DEx) * 180/pi theta_DE = 9.9250 A = (1, 2, 3) x y zB = (-1, -2, -1) Then A×B can be calculated like this (bold letters have not been used for matrix names in MATLAB): » A = [1 2 3]; » B = [-1 -2 -1]; » A_x_B = cross(A,B) %Uses MATLAB’s cross product function A_x_B = 4 -2 0 To verify this result, we can do the math term by term… » x = 1; y = 2; z = 3; %Coordinate index definitions » [ A(y)*B(z)-A(z)*B(y) -(A(x)*B(z)-A(z)*B(x)) A(x)*B(y)-A(y)*B(x) ] ans = 4 -2 0 Or, we can use the trig. form of the cross product. First, we calculate the norm, or magnitude, of the A and B vectors using the norm function in MATLAB… » A_mag = norm(A) A_mag = 3.7417 » B_mag = norm(B) B_mag = 2.4495 Then find the angle between vectors A and B using MATLAB’s acos() function. » theta = 180/pi * acos(dot(A,B)/(A_mag * B_mag)) theta = 150.7941 The magnitude of the C matrix can then be calculated… » C_mag = A_mag * B_mag * sin(theta * pi/180) C_mag = 4.4721 The direction of C is perpendicular to the plane formed by A and B, and is found using the cross product. To obtain the direction cosines of C, divide the cross product of A and B by its magnitude. » cross(A,B)/norm(cross(A,B)) ans = 0.8944 -0.4472 0 » alpha = acos(0.8944) * 180/pi %from x+ alpha = 26.5685 » beta = acos(-0.4472) * 180/pi %from y+ beta = 116.5642 » gamma = atan(0) * 180/pi %from z+ gamma = 0 This is, of course, equivalent to » C = cross(A,B); » C/C_mag ans = 0.8944 -0.4472 0 The vectors can be graphed to see how the cross product works. The plot on the left shows the original plot, with the axes oriented in the same way as the drawing on the second page. In the plot on the right the axes have been rotated to show that vector C is perpendicular to the plane formed by vectors A and B. x = [ 0 0 0; y = [ 0 0 0; z = [ 0 0 0; 1 -1 4]; 2 -2 2]; 3 -1 0]; Annotated MATLAB Script Solution %Define the vectors A = [1 2 3]; B = [-1 -2 -1]; %Take the cross product A_x_B = cross(A,B); fprintf('A x B = [ %1.4f %1.4f %1.4f]\n\n', A_x_B) %Check MATLAB's cross product operator by calculating the cross product explicitly... x = 1; y = 2; z = 3; % Define coordinate index definitions A_x_B_exp = [A(y)*B(z)-A(z)*B(y) -(A(x)*B(z)-A(z)*B(x)) A(x)*B(y)-A(y)*B(x)]; fprintf('A x B calculated explicitly= [ %1.4f %1.4f %1.4f]\n\n', A_x_B_exp) %Use the trigonometric form of the cross product operator to find the magnitude of the C vector. % First, find the magnitude of the A & B vectors using the norm function. A_mag = norm(A); B_mag = norm(B); fprintf('Magnitude of vector A = %1.4f \n', A_mag) fprintf('Magnitude of vector B = %1.4f \n', B_mag) % Then, find the angle between vectors A and B. theta = 180/pi * acos(dot(A,B)/(A_mag * B_mag)); fprintf('Angle between vectors A and B = %1.4f deg\n', theta) % Finally, solve for the magnitude of the C vector. C_mag = A_mag * B_mag * sin(theta * pi/180); fprintf('Magnitude of vector C = %1.4f \n\n', C_mag) %Solve for the direction of the C vector cross(A,B)/norm(cross(A,B)); % or C/C_mag where C = cross(A,B) alpha = acos(0.8944) * 180/pi; beta = acos(-0.4472) * 180/pi; gamma = atan(0) * 180/pi; fprintf('alpha = %1.4f deg from +x\n', alpha) fprintf('beta = %1.4f deg from +y\n', beta) fprintf('gamma = %1.4f deg from +x\n', gamma) M_mag = 64.8062 However, the moment of force F about the x axis requires an additional dot product with a unit vector in the x-direction, and is found as » u_x = [1 0 0]; » M_L = dot(u_x, cross(r,F)) M_L = -62.8000 The minus sign indicates that the moment is directed in the –x direction. Solution Using Scalars The moment of force F about the x axis can also be determined by multiplying the y-component of F and the perpendicular distance between the point at which F acts and the x axis. ML = Fy d Since the component of F in the y-direction is known (157 N), and the perpendicular distance is 0.4 m, the moment can be calculated from these quantities. » d = 0.40; » M_L = F(2) * d M_L = 62.8000 The direction must be determined using the right-hand rule, where the thumb indicates the direction when the fingers are curled around the x axis in the direction of the rotation caused by Fy. Annotated MATLAB Script Solution %Define the vectors r = [0 0 0.4]; F = [-40 157 118]; %Take the cross product of r with F to get the moment about point 0 (the origin); M_o = cross(r,F); M_mag = norm(M_o); fprintf('M_o = r x F = [ %1.4f %1.4f %1.4f]\n', M_o) fprintf('M_mag = |M_o| = %1.4f\n', M_mag) %Note: This is not the solution to the stated question. % The question asks for the moment about the x axes. % That will be calculated next %Declare a unit vector in the x-direction in order to calculate the moment about the x axis. u_x = [1 0 0]; %Calculate the moment about the x axis M_L = dot(u_x, cross(r,F)); fprintf('Moment about the x axis = %1.4f\n', M_L) %Check the results using scalar math d = r(3); M_L = F(2) * d; fprintf('Moment about the x axis (with scalar math) = %1.4f\n', M_L) Moment of a Couple Ref: Hibbeler § 4.6, Bedford & Fowler: Statics § 4.4 A couple is a pair of forces, equal in magnitude, oppositely directed, and displaced by perpendicular distance, d. FB (= -FA) d FA Since the forces are equal and oppositely directed, the resultant force is zero. But the displacement of the force couple (d) does create a couple moment. The moment, M, about some arbitrary point O can be calculated. d FA FB O rA rB ( )ABAA BBAA FrFr FrFrM −×+×= ×+×= If point O is placed on the line of action of one of the forces, say FB, then that force causes no rotation (or tendency toward rotation) and the calculation of the moment is simplified. FA FBO r AFrM ×= This is a significant result: The couple moment, M, depends only on the position vector r between forces FA and FB. The couple moment does not have to be determined relative to the location of a point or an axis. 5 Annotated MATLAB Script Solution %Solution 1 % %Define the vectors F_A = [ 106.06 -106.06 0]; %Newtons r_A = [ 28.284 28.284 0]; %centimeters F_B = [-106.06 106.06 0]; %Newtons r_B = [-28.284 -28.284 0]; %centimeters %Compute the moment using MATLAB's cross product function M = cross(r_A,F_A) + cross(r_B,F_B); %Newton centimeters %Convert to Newton meters M = M/100; %Newton meters fprintf('The resulting couple moment is = [ %1.2f %1.2f %1.2f ] (Nm)\n', M) %Solution 2 % %Define the vectors F_A = [ 0 -150 0]; %Newtons r = [ 80 0 0]; %centimeters %Compute the moment M = cross(r,F_A); %Newton centimeters %Convert to Newton meters M = M/100; %Newton meters fprintf('The resulting couple moment is = [ %1.2f %1.2f %1.2f ] (Nm)\n', M) %Solution 3 % %Define the applied force and the perpendicular distance between the forces F = 150; %Newtons d = 80; %centimeters %Compute the moment M= d * F; %Newton centimeters %Convert to Newton meters M = M/100; %Newton meters fprintf('The resulting couple moment is = %1.2f (Nm)\n', M) fprintf(' The direction must be determined using the right-hand rule.\n\n') Example B: Moment from the Small Hand Wheel The moment resulting from applying the same forces to the smaller hand wheel can be determined using any of the three solution procedures outlined above. Solution 2 is shown here. FA = 150 N FB = 150 N r x y %Example B: Moment from the Small Hand Wheel % The moment resulting from applying the same forces to the % smaller hand wheel can be determined using any of the three % solution procedures above. Solution 2 is shown here. % %Define the vectors F_A = [ 0 -150 0]; %Newtons r = [ 30 0 0]; %centimeters %Compute the moments M = cross(r,F_A); %Newton centimeters %Convert to Newton meters M = M/100; %Newton meters fprintf('The resulting couple moment for the smaller wheel is = [ %1.2f … %1.2f %1.2f ] (Nm)\n', M) Reduction of a Simple Distributed Loading Ref: Hibbeler § 4.10, Bedford & Fowler: Statics § 7.3 Loads are often not applied to specific points, but are distributed across a region. Design calculations can be simplified by finding a single equivalent force acting at a point; this is called reduction of a distributed loading. For loads distributed across a single direction (beam loading, for example) we need to find the both the magnitude of the equivalent force, and the position at which the equivalent force acts. Example 1: Load Distribution Expressed as a Function of Position, x The load on a beam is distributed as shown below: p(x) = b0 + b1 x + b2 x2 + b3 x3 + b4 x4 where b0 through b4 are coefficients obtained by fitting a polynomial to data values (see Example 2). The x values range from 0.05 to 0.85 meters. The pressure is expressed in Pa. The values of the coefficients are tabulated below. b x10-4 b0 0.107 b1 -1.68 b2 11.9 b3 -19.9 b4 9.59 The data represent the pressure on the floor beneath a pile of 12-foot 2x4’s, so y = 12 feet or 3.66 m. Determine: a) the two-dimensional loading function, w(x) b) the magnitude and position of the equivalent force. Solution 1 The two-dimensional loading function is obtained from the pressure function using the length of the boards, y = 3.66 m. ( )( )66.3xbxbxbxbb y)x(p)x(w 4 4 3 3 2 210 ++++= = The units on w are Pa.m. The magnitude of the resultant force is obtained by integration. 6 The integral for force is approximated as a summation: ∑ = ∆≈ N 1i iR xwF where ∆x is the width of a stack (0.10 m) and N is the number of stacks (9). » delta_x = 0.10; %meters » F_R = sum(w .* delta_x) %Newtons F_R = 6277 Similarly, the centroid is approximated as follows: » x_loc = sum(x .* w .* delta_x) / F_R %meters x_loc = 0.4900 Annotated MATLAB Script Solution function w = LoadingFunction(x, y, b) %Saved as LoadingFunction.m in the MATLAB search path. %Pressure Function (Pa) p = b(1) + b(2).*x + b(3).*x.^2 + b(4).*x.^3 + b(5).*x.^4; %Loading Function (Pa m) w = p .* y; function mean = FirstMoment(x, y, b) %Saved as FirstMoment.m in the MATLAB search path. %Pressure Function (Pa) p = b(1) + b(2).*x + b(3).*x.^2 + b(4).*x.^3 + b(5).*x.^4; %Loading Function (Pa m) w = p .* y; mean = w .* x; %Length of Boards (meters) y = 3.66; %Coefficients obtained by fitting a polynomial to data values b = [ 1.07E3 -1.68E4 1.19E5 -1.99E5 9.59E4 ]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Solution 1 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %The magnitude of the resultant force is calculated. F_R = quad('LoadingFunction',0.05,0.85,[],[], y, b); fprintf('\nThe magnitude of the resultant force is = %1.0f N\n', F_R) %The location of the resultant force is calculated. x_loc = quad(' FirstMoment',0.05,0.85,[],[], y, b) ./ F_R; fprintf('The location of the resultant force is = %1.2f m\n', x_loc) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Solution 2 Using Tabulated Data % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x = [ 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85]; %meters p = [ 343 1029 1372 2401 3087 4116 2744 1372 686]; %Pa %Calculate w values from pressure values using the length of the boards. w = p .* y; %Calculate the magnitude of the resultant force (the use % of a summation instead of an integration % makes this an approximate result). delta_x = 0.10; % meters F_R = sum(w .* delta_x); fprintf('\nThe approximated magnitude of the resultant force is = %1.0f N\n', F_R) %Calculate the location of the resultant force x_loc = sum(x .* w .* delta_x) / F_R; fprintf('The approximated location of the resultant force is = %1.2f m\n', x_loc) A Final Note Both of these solutions are approximations. The second method approximates the integrals using summations, while the first method approximates the solution because the polynomial regression is a “best fit” but not a perfect fit to the data. In the graph shown here, the circles represent the tabulated values in Example 2, while the pressures predicted by the polynomial are indicated by the line. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 500 1000 1500 2000 2500 3000 3500 4000 4500 x (m) P re ss ur e (P a) Tabulated Predicted » Fsp_y = F_spring * sin(alpha) %Newtons Fsp_y = 17.1569 At equilibrium the sum of the moments at A must be zero. We can use this to solve for the applied force, F. » F = -(abs(Fsp_x) * 16 + Fsp_y * 120) / 150 %Newtons F = -17.1673 Note: The absolute value function was used on Fsp_x since the direction of rotation was accounted for as the equilibrium equation was written. Here, counter-clockwise rotation was assumed positive. The equilibrium relationships for the x and y components of force can be used to determine the reactions at A. First, x-component equilibrium requires that the sum of the x components of force be zero. This relationship is used to solve for Ax. » A_x = -Fsp_x %Newtons A_x = 32.2674 Finally, the equilibrium relationship for the sum of the y components of force is used to calculate Ay. » A_y = -Fsp_y – F %Newtons A_y = 0.0105 Annotated MATLAB Script Solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Spring Support Problem % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Create conversion factor from degrees to radians deg2rad = pi /180; %Find the total extended lngth of the spring assembly. L_ext = 120 / cos( 28 * deg2rad ); %mm fprintf('\nTotal extended length = %1.0f mm\n', L_ext) %Calculate the spring's extension. L_unstretched = 75; %mm L_stretch = L_ext - L_unstretched; %mm fprintf('Spring''s extension = %1.0f mm\n\n', L_stretch) %Calculate spring force. k_spring = 600; %N/m F_spring = k_spring * L_stretch / 1000; %Newtons fprintf('Spring force = %1.1f N\n\n', F_spring) %Find x and y components of spring force. alpha = (180 - 28) * deg2rad; %angle from +x axis Fsp_x = F_spring * cos(alpha); %Newtons Fsp_y = F_spring * sin(alpha); %Newtons fprintf('x-component of spring force = %1.1f N', Fsp_x) fprintf('\t<= acting in -x direction\n') fprintf('y-component of spring force = %1.1f N\n\n', Fsp_y) %Use the sum of the moments about A to find applied force F %We know => 16 |Fsp_x| + 120 Fsp_y+ 150 F = 0 %Solving for F F = -(abs(Fsp_x) * 16 + Fsp_y * 120) / 150; %Newtons fprintf('Applied force = %1.1f N', F) fprintf('\t<= acting to cause clockwise rotation\n\n') %Use x and y component force balences to solve for the reactions at A %We know => A_x + Fsp_x = 0 %Solving for A_x A_x = -Fsp_x; %Newtons fprintf('A_x = %1.1f N', A_x) fprintf('\t<= acting in +x direction\n') %We know => A_y + Fsp_y +F = 0 %Solving for A_y A_y = -Fsp_y - F; %Newtons fprintf('A_y = %1.1f N\n', A_y) Dry Friction Ref: Hibbeler § 8.2, Bedford & Fowler: Statics § 9.1 To move a heavy crate, such as the one illustrated in the example below, you have to push on it. If you push hard enough to overcome friction, the crate will slide along the floor. But there is the danger that the crate will tip, especially if you push too high on the side of the crate. To better understand the calculations involved with this dry friction problem, and to see how you can check for tipping, we will solve the “sliding a crate” problem for several cases, varying the push height between 0.2 and 1 m, and the magnitude of the pushing force between 100 and 120 N. Example: Pushing a Crate A force, P, is applied to the side of a crate at height yP from the floor. The crate mass is 35 kg, and is 0.5 m wide and 1.0 m in height. The coefficient of static friction is 0.32. For each of the test cases tabulated below, determine: a) the magnitude of the friction force, F, b) the magnitude of the resultant normal force, NC, and c) the distance, x, (measured from the origin indicated in the illustration below) at which the resultant normal force acts. Once these results have been calculated, check for tipping and/or slipping. yP P 0.25 m 0.25 m 1 m F NC x y O W x Test Cases P yP 1 100 N 0.2 m 2 100 N 0.4 m 3 100 N 0.8 m 4 100 N 1.0 m 5 120 N 0.2 m 6 120 N 0.4 m 7 120 N 0.8 m 8 120 N 1.0 m 8 Note: The if statements are certainly not required here – the test could be performed by inspection – but it is convenient to let MATLAB do the checking for each of the eight cases. From these tests we can see that a push of 100 N applied at a height of 0.2 m will not cause the crate to tip, but it will not cause the crate to move (slip), either. Solution, Case 2 The only change between Case 1 and Case 2 is the height at which the push is applied. In Case 2 yP = 0.4 m. This value is change at the beginning of the MATLAB script, and the script is simply re- executed to calculate the values for Case 2. The complete, annotated script solution for Case 2 is shown below. Annotated MATLAB Script Solution – Case 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Pushing a Crate Problem – Case 2 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Data from the problem statement... % for specific case: Case = '1'; P = 100; %Newtons y_p = 0.4; %meters % for all cases: M = 35; %kg x_crit = 0.25; %meters mu_s = 0.32; g = 9.807; %meters/s^2 %Calculate the weight of the crate. W = -M .* g; %Newtons fprintf('\nW = %1.1f N\t\t', W) %EQUILIBRUM: The sum of the x components of force is zero. %P + F = 0 so... F = -P; %Newtons fprintf('F = %1.1f N\t\t', F) %EQUILIBRUM: The sum of the y components of force is zero. %W + N_c = 0 so... N_c = -W; %Newtons fprintf('N_c = %1.1f N\n', N_c) %EQUILIBRUM: The sum of the moments about 0 is zero. %-P * y_p + N_c * x = 0 so... x = P .* y_p ./ N_c; %meters fprintf('x = %1.3f m\t\t', x) %Checking the Results... F_max = mu_s .* N_c; %Newtons fprintf('F_max = %1.1f N\n\n', F_max) if (abs(F) > F_max) fprintf(['Case ', Case, ' is slipping\n']) else fprintf(['Case ', Case, ' is NOT slipping\n']) end if (x > x_crit) fprintf(['Case ', Case, ' is tipping\n']) else fprintf(['Case ', Case, ' is NOT tipping\n']) end Summary of Results for All Eight Cases By repeatedly changing the assigned values for P and yP in the MATLAB script, the results for all cases can be determined. Those results are summarized here. The value NC = 343.2 N is not case dependent in this problem. Test Cases P yP F x Slip? Tip? 1 100 N 0.2 m -100 N 0.058 m No No 2 100 N 0.4 m -100 N 0.117 m No No 3 100 N 0.8 m -100 N 0.233 m No No 4 100 N 1.0 m -100 N 0.291 m* No Yes 5 120 N 0.2 m -120 N 0.070 m Yes No 6 120 N 0.4 m -120 N 0.140 m Yes No 7 120 N 0.8 m -120 N 0.280 m* Yes Yes 8 120 N 1.0 m -120 N 0.350 m* Yes Yes * The calculated result for x has no meaning except to indicate that the crate is tipping if x > xcrit. A faster approach to solving the problem is to define P and yP as vectors and having MATLAB perform all the calculations simultaneously. By default MATLAB performs matrix calculations when using the *, /, and ^ operators. To force MATLAB to perform element-by-element operations place a period in front of the operator. To display the results we can then take advantage of the MATLAB for() statement to loop through the solution vectors and display the results by Case. Annotated MATLAB Script Solution – All Cases %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Pushing a Crate Problem - All Cases % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Data from the problem statement... % for specific case: Case = [ 1 2 3 4 5 6 7 8 ]; P = [ 100 100 100 100 120 120 120 120]'; %Newtons y_p = [ 0.2 0.4 0.8 1.0 0.2 0.4 0.8 1.0]'; %meters % for all cases: M = 35; %kg x_crit = 0.25; %meters mu_s = 0.32; g = 9.807; %meters/s^2 %Calculate the weight of the crate. W = -M .* g; %Newtons %EQUILIBRUM: The sum of the x components of force is zero. %P + F = 0 so... F= -P; %Newtons %EQUILIBRUM: The sum of the y components of force is zero. %W + N_c = 0 so... N_c = -W; %Newtons %EQUILIBRUM: The sum of the moments about 0 is zero. %-P * y_p + N_c * x = 0 so... x = P .* y_p ./ N_c; %meters %Display and Check the Results... F_max = mu_s .* N_c; %Newtons fprintf('F_max = %1.1f N\t', F_max) fprintf('\nW = %1.1f N\t', W) fprintf('N_c = %1.1f N\n\n', N_c) %Display results for each case. for i = 1 : length(Case) fprintf('\nCASE %1.0f:\n',Case(i)) fprintf('P = %1.1f N\t\t', P(i)) fprintf('y_p = %1.2f m\n', y_p(i)) fprintf('F = %1.1f N\t\t', F(i)) fprintf('x = %1.3f m\n', x(i)) if (abs(F(i)) > F_max) fprintf(['Case ', num2str(Case(i)), ' is slipping\n']) else fprintf(['Case ', num2str(Case), ' is NOT slipping\n']) end if (x > x_crit) fprintf(['Case ', num2str(Case(i)), ' is tipping\n']) else fprintf(['Case ', num2str(Case(i)), ' is NOT tipping\n']) end end yx z Vi = π ri 2 ∆z ri ∆z zi Here, the value of r at the current (point “i”) value of z has been used to calculate the volume of the disk. Since R = 7 cm in this example, we might try a ∆z of 1 cm, and calculate the volumes of seven disks (N = 7). For each disk, the value of zi can be calculated using zizi ∆= The value of the corresponding radius is then 2 i 2 i zRr −= The formula for volume can then be written in terms of z, as ( ) zzRV 2 i 2 i ∆−π= The sum of all of the volumes of all seven disks is an approximation of the volume of the hemisphere. ( )[ ]∑∑ == ∆∆−π=≈ N 1i 22 N 1i i zziRVV In MATLAB, this calculation looks like this: » R = 7; %Raidus » N = 7; %Number of intervals » deltaZ = R/N; %Interval spacing » %Sum (pi * (R^2 - (i*deltaZ)^2) * deltaZ) for i = 1,2,...,N » V_predicted = 0; %Initialize predicted volume to zero » for i = 1:N V_predicted = V_predicted + pi * (R^2 - (i*deltaZ)^2) * deltaZ; end » V_predicted V_predicted = 637.7433 We can see how good (or bad) the approximation is by calculating the actual volume of the hemisphere. V_actual = 2/3 * pi * R^3 V_actual = 718.3775 The approximation using only seven disks is not too good. If we use more disks, say N = 70, the approximation of the volume is quite a bit better. » R = 7; %Raidus » N = 70; %Number of intervals » deltaZ = R/N; %Interval spacing » %Sum (pi * (R^2 - (i*deltaZ)^2) * deltaZ) for i = 1,2,...,N » V_predicted = 0; %Initialize predicted volume to zero » for i = 1:N V_predicted = V_predicted + pi * (R^2 - (i*deltaZ)^2) * deltaZ; end » V_predicted V_predicted = 710.6440 To calculate the centroid using numerical methods, simply replace the ratio of integrals by the numerical approximation: ( )[ ] ( )[ ]∑ ∑ ∫ ∫ = = ∆∆−π ∆∆−π ≈= N 1i 22 N 1i 22 V V zziR zziRz dV dVz z The extra “z” has been included in the numerator as (I * deltaZ). In MATLAB, the calculation of the centroid looks like this: » numerator = 0; %Initialize numerator to zero » denominator = 0; %Initialize denominator to zero » for i = 1:N %Sum both numerator and denominator numerator = numerator + (i*deltaZ) * pi * (R^2 - (i*deltaZ)^2) * deltaZ; denominator = denominator + pi * (R^2 - (i*deltaZ)^2) * deltaZ; end » z_bar = numerator / denominator z_bar = 2.6530 We can use the analytical result to check our answer. » z_bar_act = 3/8 *R z_bar_act = 2.6250 Annotated MATLAB Script Solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Define the System %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R = 7; %Radius N = 70; %Number of disks deltaZ = R/N; %Calculate disk thickness fprintf('SYSTEM\n') fprintf('Raidus = %3.1f cm\t\tN = %3.0f\t\tdeltaZ = %3.2f\n\n',R,N,deltaZ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Calculate the Volume %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Sum (pi * (R^2 - (i*deltaZ)^2) * deltaZ) for i = 1,2,...,N V_predicted = 0; %Initialize predicted volume to zero for i = 1:N V_predicted = V_predicted + pi * (R^2 - (i*deltaZ)^2) * deltaZ; end fprintf('CALCULATE THE VOLUME\n') fprintf('Predicted Volume = %3.3f cm^3\n',V_predicted) V_actual = 2/3 * pi * R^3; fprintf('Actual Volume = %3.3f cm^3\n\n',V_actual) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Calculate the Centroid %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% numerator = 0; %Initialize numerator to zero denominator = 0; %Initialize denominator to zero for i = 1:N %Sum both numerator and denominator numerator = numerator + (i*deltaZ) * pi * (R^2 - (i*deltaZ)^2) * deltaZ; denominator = denominator + pi * (R^2 - (i*deltaZ)^2) * deltaZ; end z_bar = numerator / denominator; fprintf('CALCULATE THE CENTROID\n') fprintf('Calculated z_bar = %3.3f cm\n',z_bar) z_bar_act = 3/8 *R; fprintf('Actual z_bar = %3.3f cm\n',z_bar_act) Annotated MATLAB Script Solution function p = PressureFunction(x, y) %Saved as PressureFunction.m in the MATLAB search path. %Pressure Function (Pa) p = 1000 + 230 .*x - 210 .*x.^2 + 120 .* y - 70 .* y.^2; function FM = FirstMomentX(x, y) %Saved as FirstMomentX.m in the MATLAB search path. %Pressure Function (Pa) FM = x .* (1000 + 230 .*x - 210 .*x.^2 + 120 .* y - 70 .* y.^2); function FM = FirstMomentY(x, y) %Saved as FirstMomentY.m in the MATLAB search path. %Pressure Function (Pa) FM = y .* (1000 + 230 .*x - 210 .*x.^2 + 120 .* y - 70 .* y.^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Resultant of a Generalized Distributed Loading % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %The magnitude of the resultant force is calculated. F_R = dblquad('PressureFunction', 0, 1, 0, 1);%Newtons fprintf('\nThe magnitude of the resultant force is = %1.0f N\n', F_R) %The location of the resultant force is calculated. x_loc = dblquad('FirstMomentX', 0, 1, 0, 1) ./ F_R;%meters y_loc = dblquad('FirstMomentY', 0, 1, 0, 1) ./ F_R;%meters fprintf('The location of the resultant force is = (%1.3f, %1.3f)m\n', x_loc, y_loc) Calculating Moments of Inertia Ref: Hibbeler § 10.1-10.2, Bedford & Fowler: Statics § 8.1-8.2 Calculating moments of inertia requires evaluating integrals. This can be accomplished either symbolically, or using numerical approximations. Mathcad’s ability to integrate functions to generate numerical results is illustrated here. Example: Moment of Inertia of an Elliptical Surface Determine the moment of inertia of the ellipse illustrated below with respect to a) the centroidal x’ axis, and b) the x axis. The equation of the ellipse relative to centroidal axes is 1 14 'y 8 'x 2 2 2 2 =+ C O d = dy = 16 y' -x x = x' dy' y x' x In this problem, x and y have units of cm. Solution The moment of inertia about the centroidal x axis is defined by the equation ∫= A 2 'x dA'yI where dA is the area of the differential element indicated in the figure above. 'dyx2dA = So, the integral for moment of inertia becomes 11 ∫= A 2 'x 'dyx2'yI Furthermore, x (or x’) can be related to y’ using the equation of the ellipse. Note: Because of the location of the axes, x = x’ in this example.         −== 2 2 2 14 'y18'xx The equation for the moment of inertia becomes: ∫−         −= 8 8 2 2 22 'x 'dy 14 'y182'yI To perform this integration we need to place the integrand in an m-file function and call MATLAB’s quad() function on the m-file. function Ix_integrand = Moment_Of_Inertia_Integrand(y_prime) %Saved as Moment_Of_Inertia_Integrand.m in the MATLAB %search path. x = sqrt(8.^2 .* (1 - y_prime.^2./14^2)); Ix_integrand = y_prime.^2 .* 2 .* x; » Ix_prime = quad('Moment_Of_Inertia_Integrand',-8,8) %cm^4 Ix_prime = 4890 The moment of inertia relative to the original x axis can be found using the parallel-axis theorem. 2 y'xx dAII += Where A is the area of the ellipse, and dy is the displacement of the centroidal y axis from the original y axis. The required area can be calculated by integration in the same fashion as before. function A = Area_Integrand(y_prime) %Saved as Area_Integrand.m in the MATLAB search path. x = sqrt(8.^2 .* (1 - y_prime.^2 ./ 14^2)); A = 2 .* x; » A = quad('Area_Integrand',-8,8) %cm^2 A = 241.2861 Then the moment of inertia about x can be determined. » dy = 16; %cm 0 2 4 6 8 10 12 14 16 18 20 0 10 20 30 40 50 60 70 80 t(sec) v (m ph ) The graph shows the increasing velocity for the first 14 seconds then the velocity drops off again. It looks like a smooth curve, and clean (not noisy) data. A commonly used approach to fitting a curve to a set of data values is polynomial regression. We might try, for example, the following second-order polynomial.1 32 2 1 btbtbvp ++= MATLAB’s Statistics Toolbox provides a function for fitting linear2 and nonlinear models such as this polynomial, the nlinfit() function. This function calculates the coefficient values that provide the best fit of the model equation to the data. The regression model is described by either an inline function or an m-file containing the desired function of the independent variable (t, in this case) and coefficients . For our second-order polynomial, the vector could be written as follows: » F = inline('b(1).*t.^2 + b(2).*t ', 'b', 't') F = Inline function: F(b,t) = b(1).*t.^2 + b(2).*t Note: We want the velocity to be zero when t = 0. The way to force this to happen is to leave off the constant term in the F(b,t) function. Then, nlinfit() uses the t and v values (as vectors), the F(b,t) function shown above, and a vector of initial guesses for each coefficient in F to calculate the b values that produce the “best” fit of the polynomial to the data. » b = nlinfit(t,v,F,[0 0]); » b b = 1 A second-order polynomial is not a good choice for this data. A parabola is a second-order equation, and the data shows a lot more curvature than a simple parabola. It will require a higher-order polynomial to adequately fit this data – but the goal of the first attempt is to demonstrate that the “best fit” polynomial may still be a very poor fit to the data. 2 For curve-fitting, the variables are the model coefficients, the b values. All of the times are known from the data set. This polynomial is linear in the b values, so it is a linear model for curve fitting. -0.2706 8.2273 These coefficients tell us that the regression model is 02273.82706.0 2 ++−= ttvp The subscript p indicates that the equation is used to calculate predicted values of velocity. The calculated coefficients provide the best possible fit of the model equation to the data, but that does not guarantee that the equation is a good fit. You should always calculate predicted values with the regression equation and compare them to the original data to visually determine if the regression equation is actually fitting the data. » v_p = b(1).*t.^2 + b(2).*t; » plot(t,v,'o',t,v_p); » xlabel(' t(sec) '); » ylabel(' v(mph) '); » legend('v','v_p',2) Clearly, the regression curve (solid line) is not doing a good job of fitting the data. Solution – Second Attempt We need a better regression equation. We’ll try a third-order polynomial. 43 2 2 3 1 btbtbtbvp +++= The inline function F(b,t) gets another element… » F = inline('b(1).*t.^3 + b(2).*t.^2 + b(3).*t','b','t'); …and nlinfit() now finds three b values… » b = nlinfit(t,v,F,[0 0 0]); » b b = -0.0493 1.0762 -0.0403 0 2 4 6 8 10 12 14 16 18 20 0 10 20 30 40 50 60 70 80 t(sec) v (m ph ) v vp Again, we calculate predicted velocities to check the fit. » v_p = v_p = b(1).*t.^3 + b(2).*t.^2 + b(3).*t; » plot(t,v,'o',t,v_p); » xlabel(' t(sec) '); » ylabel(' v(mph) '); » legend('v','v_p',2) The new model is doing a good job of fitting the data. The regression equation that fits the data is 32 p t049.0t076.1t04.00v −+−= This equation can be used to generate the s-t and a-t graphs. MATLAB can integrate the vp equation to get position, s, and differentiate the equation for acceleration, a, using the polyint() and polydiff() functions respectively. MATLAB represents polynomial’s as vectors of coefficients. The first vector element is the coefficient on the highest ordered term and the last vector element is the constant. Therefore, we need to add a zero to the end of our b vector so it is treated as a 3rd order polynomial by MATLAB. » b(4) = 0 b = -0.0493 1.0762 -0.0403 0 » b(4) = 0; » b_a = polyder(b) b_a = -0.1480 2.1524 -0.0403 » b_s = polyint(b) b_s = -0.0123 0.3587 -0.0202 0 0 %Calculate predicted velocity values (at each t value) to see if the regression equation actually fits the data. v_p = b(1).*t.^3 + b(2).*t.^2 + b(3).*t; subplot(1,2,2); plot(t,v,'o',t,v_p); xlabel('{\bf t(sec)}'); ylabel('{\bf v(mph)}'); legend('v','v_p',2); %Convert b to a MATLAB polynomial vector b(4) = 0; %Take the derivative of b to find the equation for a. b_a = polyder(b); %Calculate a at each value of t a = b_a(1).*t.^2 + b_a(2).*t + b_a(3); %Simplify the units by building in the conversion from hours to seconds. a = a ./ 3600; %Integrate b to find the equation for s. b_s = polyint(b) %Calculate s at each value of t s = b_s(1).*t.^4 + b_s(2).*t.^3 + b_s(3).*t.^2 + b_s(4)*t; %Simplify the units by building in the conversion from hours to seconds. s = s ./ 3600; %Plot the results in a new figure. figure; subplot(1,2,1); plot(t,a,'o'); xlabel('{\bf t}'); ylabel('{\bf a}'); subplot(1,2,2); plot(t,s,'o'); xlabel('{\bf t)}'); ylabel('{\bf s}'); Curvilinear Motion: Rectangular Components Ref: Hibbeler § 12.5, Bedford & Fowler: Dynamics § 2.3 A fixed x, y, z frame of reference is commonly used to describe the path of a particle in two situations: 1. The path “fits” the Cartesian coordinate system well (tends to follow straight lines, or simple paths.) 2. The path does not “fit” any commonly used coordinate system well. The example used here falls into the latter category. Example: Finding the Velocity and Acceleration of a Particle The complex path of a particle for times between 0 and 10 seconds has been fit using multivariable regression to the following functions: 3 32 t14.04.2)t(z t031.0t7.23.1)t(y )tcos(2)t(x += −+= = The position of the particle at any point in time for which the functions are valid (0 ≤ t ≤ 10 sec.) can be determined as kjir zyx ++= Given these functions for x, y, and z determine the position, and the magnitudes of the velocity and acceleration of the particle at t = 7 sec. Solution The first thing we might want to do is graph this particle path, just to see what it looks like. MATLAB can help here. First, declare the three functions of time that describe the particle’s path: » x = inline('2 .* cos(t)', 't') x = Inline function: x(t) = 2 .* cos(t) » y = inline('1.3 + 0.27.*t.^2 - 0.031.*t.^3', 't') y = Inline function: y(t) = 1.3 + 0.27.*t.^2 - 0.031.*t.^3 » z = inline('2.4 + 0.14.*t.^3', 't') z = Inline function: z(t) = 2.4 + 0.14.*t.^3 13 Then ask MATLAB to plot the particle path by placing the three function names into the scatter3() function. » t = 0 : 0.1: 10; %sec - time over which to plot » scatter3(x(t),y(t),z(t),'o') » xlabel('x-axis') » Ylabel('y-axis') » zlabel('z-axis') Solving for Position Declare r(t) as a three-component matrix composed of the three functions. The position is described by function r(t), so finding the position at t = 7 sec. Simply requires evaluating r(t) at this time. » r = inline('[ [2 .* cos(t)] [1.3 + 0.27. *t.^2 - 0.031 .*t.^3] [2.4 + 0.14. * t.^3] ]', 't') r = Inline function: r(t) = [[2.*cos(t)] [1.3+0.27.*t.^2-0.031.*t.^3] [2.4+0.14.*t.^3]] » r(7) ans = 1.5078 3.8970 50.4200 So, at t = 7 seconds, the particle is located at x = 1.5078, y = 3.8970, and z = 50.4200. Solving for the Magnitude of the Velocity For velocity we need to take the derivative of the x, y, and z functions. Here, the derivatives are obtained using MATLAB’s symbolic math toolbox and cutting-and-pasting into three new inline functions: vx(t), vy(t), and vz(t). Note: Capital letters are used to represent symbolic variables while lower case letters are assigned a value or to a function. » syms T %Define a symbolic variable T » X = 2 .* cos(T); %Define X » diff(X) %Differentiate on X ans = -2*sin(T) v_y = inline('27/50 - 93/1000 .*t.^2', 't'); %Assign v_y(t) Z = 2.4 + 0.14.*T.^3; %Define Z diff(Z) %Differentiate on Z v_z = inline('21/50 .* t.^2', 't'); %Assign v_z(t) %Calculate the magnitude of the velocity at 7 seconds. v = sqrt(v_x(7).^2 + v_y(7).^2 + v_z(7).^2); fprintf('Magnitude of the velocity at 7 seconds = %1.3f \n', v) %Take the time derivatives of the velocity functions. diff(-2*sin(T), T) %Differentiate on v_x a_x = inline('-2.*cos(t)','t'); %Assign a_x(t) diff(27/50 .* T - 93/1000 .* T.^2, T); %Differentiate on v_y a_y = inline('0.54 - 0.186.*t','t'); %Assign a_y(t) diff(21/50 .* T.^2, T) %Differentiate on v_z a_z = inline('0.84.*t','t'); %Assign a_z(t) %Calculate the magnitude of the acceleration at 7 seconds. a = sqrt(a_x(7).^2 + a_y(7).^2 + a_z(7).^2); fprintf('Magnitude of the acceleration at 7 seconds = %1.3f \n\n', a) Curvilinear Motion, Motion of a Projectile Ref: Hibbeler § 12.6, Bedford & Fowler: Dynamics § 2.3 Rectilinear motion refers to motion in a straight line. When a particle follows a non-straight path, it’s motion is termed curvilinear. Projectile motion is typically curvilinear, although a projectile fired straight up (in the absence of a crosswind), or moving along a straight track would be rectilinear motion. A projectile’s motion can be broken down into three phases: an acceleration phase where the actuator (gun, catapult, golf club, etc.) gets the projectile moving. The second phase of motion is after the projectile leaves the actuator, when the only acceleration acting on it is the acceleration due to gravity. Note: A common assumption that simplifies the problem considerably, but is not altogether accurate, is that the frictional drag between the projectile and the fluid through which it moves is negligible. This assumption is more reasonable for small, smooth, slow-moving particles through low-viscosity fluid than for large, irregularly shaped particles moving at high speeds through highly viscous fluids. The third phase of a projectile’s motion is after impact. The particle may roll, continue moving through a very different medium (water or earth), or break up. Here we will consider only the second phase of the projectile’s motion: the projectile has already been accelerated by an actuator, and is flying through the air. Example: Slingshot Contest Projectile throwing contests are pretty common at engineering schools. There are several ways to run the contest: highest, farthest, most accurate, etc; and several ways to propel the projectile: slingshot, catapult, etc. This example focuses on a tennis ball slingshot competition. Tennis Ball Data (approximate) Diameter 2.5 inches (6.4 cm) Weight 2 oz (57 g) A team of mechanical engineering students has built a slingshot that allows precise control of the pre- release tension and angle. They also connected the release mechanism to a digital camera to take a series of photos in 1millisecond intervals just as the tennis ball is released. The photos are used to calculate the velocity (speed and angle) at which the tennis ball leaves the sling. 14 A test run with the camera operational gave the following set of photos (superimposed). 25° 0 20 40 60 mm Part 1. Determine: a. The initial velocity of the tennis ball as it leaves the sling. b. The predicted time of flight for the ball. c. The predicted horizontal travel distance for the ball. Part 2. The team that gets most balls into a basket set 30 meters from the launch site wins. If they can keep the initial speed constant (at the test run value of 22.1 m/s), what angle should they use to shoot the tennis balls into the basket? Part 1. Solution The initial velocity is calculated from the 60 mm horizontal travel distance observed in the photos using the camera’s 1 millisecond interval between snapshots. » x_test = 60 / 1000; %Horizontal travel distance for all four frames (m) » theta = 25 * pi/180; %Angle of trajectory (radians) » d_test = x_test / cos(theta) %Distance ball traveled in direction of motion (m) d_test = 0.0662 » dt_pic = 0.001; %Interval between snapshots (s) » dt_test = 3 * dt_pic; %3 time intervals between the four photos » v_int = d_test / dt_test %Initial velocity of the tennis ball v_int = 22.0676 The horizontal and vertical components of the initial velocity will be useful for later calculations. » v_yint = v_int * sin(theta) %y component of initial velocity Annotated MATLAB Script Solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Projectile Motion: Slingshot Contest %Calculate Initial Velocity from Test Run %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Horizontal travel distance for all four frames (m) x_test = 60 / 1000; %Angle of trajectory (radians) theta = 25 * pi/180; %distance ball traveled in direction of motion (m) d_test = x_test / cos(theta); fprintf('Distance ball traveled in direction of motion = %1.4fm\n',d_test) %Interval between snapshots (s) dt_pic = 0.001; %3 time intervals between the four photos dt_test = 3 * dt_pic; %Initial velocity of the tennis ball v_int = d_test / dt_test; fprintf('Initial velocity = %3.1f m/s\n',v_int) %y component of initial velocity v_yint = v_int * sin(theta); fprintf('\ty component = %3.1f m/s\n',v_yint) %x component of initial velocity v_xint = v_int * cos(theta); fprintf('\tx component = %3.1f m/s\n\n',v_xint) %x component of velocity is constant (if air resistance is ignored) v_x = v_xint; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Calculate Time of Flight %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %The constant acceleration in this problem is due to gravity a_c = -9.8; t_flight= -2 * v_yint / a_c; fprintf('Flight time (Ignoring y_o) = %3.3f s\n',t_flight) y_o = 0; %Initial height y_final = 0; %Final height y_max = y_o + v_yint * t_flight/2 + a_c/2 * ( t_flight/2 )^2; fprintf('y_max = %3.1f m\n\n',y_max) continues… %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Alternative Solution Method - Include y_o and use MATLAB's root finder %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %The predicted time of flight can be calculated using % y_max = y_o + v_yint * t_flight + a_c/2 * t_flight^2 % %This equation is a second order polynomial in t_flight and can be %rewritten in the form % 0 = a2 * x^2 + a1 * x + a0 % with % a2 = a_c/2 % a1 = v_yint % a0 = y_o - y_final % %Polynomials are represented in MATLAB by an array of the %coefficients in descending order. Therefore, this polynomial %would be represented by the following % coefficients = [a2 a1 a0] % or % coefficients = [( a_c/2 ) ( v_yint ) ( y_o - y_final )]; y_o = 20 / 100; %Initial height (m) y_final = 0; %Final height (m) coeffs = [( a_c/2 ) ( v_yint ) ( y_o - y_final )]; %Polynomial t_flight= roots(coeffs); %Roots t_flight= t_flight(1); %Choose the positive root fprintf('Alternative Solution Method - Include y_o and use MATLAB''s root finder\n') fprintf('Flight time = %3.3f s\n',t_flight) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Calculate Horizontal Motion %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x_o = 0; x = x_o + v_x * t_flight; fprintf('Horizontal Motion = %3.3f m\n',x) Part 2. With the MATLAB Script, we can simply try some different angles until we calculate a predicted horizontal travel distance of 30 meters. First, try 20 degrees. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Calculate Inital Velocity Components %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Initial velocity of the tennis ball v_int = 22.1; %Angle of trajectory (radians) theta = 20 * pi/180; fprintf('Initial velocity = %3.1f m/s\t\t',v_int) fprintf('Angle of trajectory = %3.2f deg\n',theta * 180/pi) %y component of initial velocity v_yint = v_int * sin(theta); fprintf('\ty component = %3.1f m/s\n',v_yint) %x component of initial velocity v_xint = v_int * cos(theta); fprintf('\tx component = %3.1f m/s\n\n',v_xint) %x component of velocity is constant (if air resistance is ignored) v_x = v_xint; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Calculate Time of Flight %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %The constant acceleration in this problem is due to gravity a_c = -9.8; t_flight= -2 * v_yint / a_c; fprintf('Flight time (Ignoring y_o) = %3.3f s\n',t_flight) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Calculate Horizontal Motion %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x_o = 0; x = x_o + v_x * t_flight; fprintf('Horizontal Motion = %3.3f m\n',x) MATLAB Output Initial velocity = 22.1 m/s Angle of trajectory = 20.00 deg y component = 7.6 m/s x component = 20.8 m/s Flight time (Ignoring y_o) = 1.543 s Horizontal Motion = 32.035 m With a little trial and error, the predicted angle should be 18.5 degrees. The skateboarder’s acceleration can be written in terms of normal and tangential velocities, as ( ) nt nt uu uua ρ += ρ += 2 s m s m 2 2.3 8.7 vv 2 & The radius of curvature, ρ, of the path at the point of interest (X=2.828 m, Y = 2 m) can be calculated as: 2 2 2 dx yd dx dy1 2 3              + =ρ In MATLAB, this equation is evaluated as follows: >> y = r - sqrt( r^2 - x^2 ); %Define Symbolic Equation >> dydx = diff(y,x) %First Derivative dydx = 1/(r^2-x^2)^(1/2)*x >> DYDX = inline('1/(R^2-x^2)^(1/2)*x','x','R'); %Create inline function DYDX >> dy2dx2 = diff(diff(y,x),x) %Second Derivative dy2dx2 = 1/(r^2-x^2)^(3/2)*x^2+1/(r^2-x^2)^(1/2) %Create inline function DY2DX2 >> DY2DX2 = inline('1/(r^2-x^2)^(3/2)*x^2+1/(r^2-x^2)^(1/2)','x','r'); >> RHO = (1 + DYDX(X,R)^2)^(3/2) / abs(DY2DX2(X,R)) RHO = 3.0000 The calculated radius of curvature is 3 m – which should come as no surprise since the path is a circle of radius 3 m. The equation for the acceleration of the skateboarder in terms of tangential and normal components is now: ( ) nt nt uu uua m3 2.3 8.7 vv 2 s m s m 2 2 += ρ += & The magnitude of the acceleration is found with the following calculation: θ A ut un a φ >> a = sqrt(7.8^2 + ( 3.2^2 / 3 )^2) %meters per second squared a = 8.5142 While the angle (in degrees), φ, is found using the atan() function and multiplying by the conversion factor 180o/π. >> PHI = atan(v_dot / ( v^2/rho )) *180/pi PHI = 66.3655 The angle of the acceleration in degrees from the positive x axis is: >> ANGLE = THETA + 90 + PHI ANGLE = 226.8941 Annotated MATLAB Script Solution %Define the function on which to use fzero. function F = FindMinimum(x) R=3; %3 meters radius Y=2; %2 meters off the ground %The equation describing the scateboarders path is that of % and quarter circle => x^2 + y^2 = R^2 % or similarly => y = R - sqrt(x^2 + y^2) % % Then setting the equation equal to some arbitrary variable % and solving for the root F = R - sqrt(R^2 - x^2) - Y; %Call fzero on the FindMinimum function with an initial guess of 2 X = fzero('FindMinimum',2); %Guess: 2 fprintf('X = %1.4f m .\n', X) %Define all other variables using uppercase letters to indicate the variable has been assigned a value. V_dot = 7.8; V = 3.2; R = 3; Y = 2; %Define the x, y , and r as symbolic variables. %Note: Lowercase variables are all symbolic variables while % uppercase variables are assigned a numeric value. syms x y r; %Find symbolic result to first and second derivatives y = r - sqrt( r^2 - x^2 ); %Define Symbolic Equation dydx = diff(y,x) %First Derivative DYDX = inline('1/(R^2-x^2)^(1/2)*x','x','R'); dy2dx2 = diff(diff(y,x),x) %Second Derivative DY2DX2 = inline('1/(r^2-x^2)^(3/2)*x^2+1/(r^2-x^2)^(1/2)','x','r'); %Calculate the angle at point A. THETA = atan(X) *180/pi; fprintf('The angle at point a is %1.4f degrees.\n', THETA) %Calculate the radius of curvature at X = 2.828 m and y = 2 m. RHO = (1 + DYDX(X,R)^2)^(3/2) / abs(DY2DX2(X,R)); fprintf('The radius of curvature at X = 2.828 m and y = 2 m is %1.4f m.\n', RHO) %Calculate the magnitude of acceleration. A = sqrt(V_dot^2 + (V^2/RHO)^2); fprintf('The magnitude of acceleration is %1.4f m/s^2.\n', A) %Calculate the algle of acceleration relative to the t,n coordinate system. PHI = atan(V_dot / ( V^2/RHO )) *180/pi; %Calculate the angle of acceleration relative to the positive x-axis. ANGLE = THETA + 90 + PHI; fprintf('The angle of acceleration relative to the positive x-axis is %1.4f degrees.\n', ANGLE) The length of the cord segment connected to weight A is labeled SA. There are four cord segments around the four pulleys that do not change as the weights move. These segments are shown in red in the figure, and their lengths are combined and called Lconst. The total cord length, then, is totalconstAB LLss4 =++ Differentiating this equation with respect to time yields the time rates of change of position, or velocities of each weight. AB AB vv4 00 dt ds dt ds 4 −= =++ So, if weight B moves down at 0.3 m/s, weight A will move up at 1.2 m/s. So far, these problems require only simple calculus, and there is little need for MATLAB’s calculational abilities. However, when there are several cords involved, MATLAB’s matrix math capabilities can come into play. Example: Three Cord System In the example shown below, three cords are connected around a complicated pulley system. The length of each of the cords is constant, so the sum of the segment lengths can be written for each cord. B A 0.3 m/s To simplify things a bit, we first determine which segments are of constant length, shown in gray in the following figure. Then the segments that are changing length are labeled (S1 through S6, as shown in the figure.) B A 0.3 m/s DATUM S1 S2 S3 S4 LLINK S5 S6 The following relationships involving the changing segment lengths can be written: 1) Segments S1 and S2, plus a fixed length, FG, sum to the total length of the green cord. LG. 2) Segments S3 and S4, plus a fixed length, FB, sum to the total length of the blue cord. LB. 3) Segments S5 and S6, plus a fixed length, FY, sum to the total length of the yellow cord. LY. Two additional relationships can be written: 4) Segments S1, S3, S5 and the length of the floating link, LLINK, sum to an unknown, but constant (fixed) length, F1. 5) Segments S4, S5 and the length of the floating link, LLINK, sum to an unknown, but constant (fixed) length, F2. In summary, the following five equations are written: 2LINK54 1LINK531 YY65 BB43 GG21 FLSS FLSSS LFSS LFSS LFSS =++ =+++ =++ =++ =++ Taking time derivatives, the segment lengths become velocities, and the constant terms go to zero. 00vv 00vvv 00vv 00vv 00vv 54 531 65 43 21 =++ =+++ =++ =++ =++ Since v6 is known to be 0.3 m/s (downward), these 5 equations in 5 unknowns can be solved. Including the known value, the equations become… 0vv 0vvv 3.0v 0vv 0vv 54 531 5 43 21 =+ =++ −= =+ =+ To solve these simultaneous linear equations using matrix methods in MATLAB, we rewrite the equations as a coefficient matrix, C, and a right-hand-side vector, r. Note that the apostrophe on the r vector transposes r from a row vector into a column vector. >> C = [ 1 1 0 0 0; 0 0 1 1 0; 0 0 0 0 1; 1 0 1 0 1; 0 0 0 1 1; ]; >> r = [ 0 0 -0.3 0 0 ]'; %Units => m/s We then invert the coefficient matrix, and multiply the inverted coefficient matrix and the r vector to obtain the unknown velocities. >> C_inv = inv(C) C_inv = 0 -1 -2 1 1 1 1 2 -1 -1 0 1 1 0 -1 0 0 -1 0 1 0 0 1 0 0 >> v = C_inv * r %Units => m/s v = 0.6000 -0.6000 -0.3000 0.3000 -0.3000 These velocities represent the rate at which segment lengths S1 through S5 are changing with time, and are not (generally) velocities relative to the datum. In order to determine the velocity of point A relative to the datum, the distance of point A relative to the datum (called LA in the figure below) must be written in terms of segment lengths. Kinetics of a Particle: Force and Acceleration Ref: Hibbeler § 13.4, Bedford & Fowler: Dynamics § 3.1-3.4 In an earlier example (#8) we looked at “pushing a crate” and investigated the effect of changing the height of the push and the magnitude of the push on the rigid body. Specifically, we tested to see if the crate would tip or slip. In this example, we will determine the initial acceleration of the crate and the velocity after a short time at that acceleration. Example: Pushing a Crate A crate weighs 35 kg, and is 0.5 m wide and 1.0 m in height. A force, P, is applied to the side of the crate at height yP = 0.4 m from the floor. From the earlier example (#8) we know that this force is sufficient to overcome friction, and when applied at this height will not tip the crate. The coefficient of kinetic friction is 0.24. Determine the initial acceleration and, assuming the acceleration remains constant, the expected velocity of the crate after 3 seconds. yP = 0.4 m P F NC x y O W Solution First, a free-body diagram is drawn. 17 yP P F NC x y O W We start the solution by assigning values given in the problem statement to variables: » M = 35; %kg » P = 120; %Newtons » y_p = 0.4; %meters » g = 9.807; %meters/s^2 » mu_k = 0.24; Then we calculate the weight of the crate. » W = -M .* g %Newtons W = -343.2450 We can use the equation of motion for the y-components of force to determine the resultant normal force, NC (but the result is not much of a surprise.) » %EQ of MOTION: y components of force. » % W + N_c = 0 so... » N_c = -W %Newtons - acts in the +x direction N_c = 343.2450 Knowing NC, the friction force can be determined. » F= mu_k .* N_c; %Newtons - this equation does not account for direction » F = -F %Newtons - sign changed since F acts in -x direction F = -82.3788 Next, we use the equations of motions with the information on the free-body diagram to solve for the acceleration, ax. » %EQ of MOTION: x components of force. » % P + F = M * a_x so... » a_x = (P + F) ./ M %m / s^2 a_x = 1.0749 Now the velocity after 3 seconds can be determined. » %Calculate the velocity after three seconds. » v_o = 0; %m / s » delta_t = 3; %s » a = a_x; %m / s^2 » v = v_o + a .* delta_t %m / s v = 3.2247 Annotated MATLAB Script Solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Pushing a Crate - Initial Acceleration % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Data from the problem statement... M = 35; %kg P = 120; %Newtons y_p = 0.4; %meters g = 9.807; %meters/s^2 mu_k = 0.24; %Calculate the weight of the crate. W = -M .* g; %Newtons fprintf('\nW = %1.1f N\t\t', W) %EQ of MOTION: y components of force. %W + N_c = 0 so... N_c = -W; %Newtons fprintf('N_c = %1.1f N\n', N_c) %Calculate teh friction force, F. F= mu_k .* N_c; %Newtons - this equation does not account for direction F = -F; %Newtons - Sign changed since F acts in -x direction fprintf('F = %1.1f N\t\t', F) %EQ of MOTION: x components of force. %P + F = M * a_x so... a_x = (P + F) ./ M; %m/s^2 fprintf('a_x = %1.3f m/s^2\n\n', a_x) The tangential acceleration at this location was given in the problem statement, as at = 7.8 m/s2. The normal acceleration can be calculated from the stated velocity (3.2 m/s) and the radius of curvature (3 m). » a_n = v.^2 ./ rho %m / s^2 a_n = 3.4133 The normal component of the equation of motion can be used to find NA, but first we need the normal component of the skateboarder’s weight. » W = M .* g %Newtons W = 539.3850 » W_n = -W .* sin(alpha) %Newtons - minus sign accounts for direction W_n = -179.1627 » W_t = W .* cos(alpha) %Newtons W_t = 508.7602 Then find the normal force, NA. » % N_A + W_n = M * a_n so... » N_A = M .* a_n - W_n %Newtons N_A = 366.8961 Next, use the tangential component of the equation of motion to determine the friction force. » % W_t + F = M * a_t so... » F = M .* a_t - W_t %Newtons F = -79.7602 The friction force, F, and the normal force, NA, can be used together to calculate the apparent coefficient of kinetic friction for this problem. » mu_k = abs(F) ./ N_A %need only the magnitude of F here mu_k = 0.2174 Annotated MATLAB Script Solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Friction Force on a Skateboarder % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Set problem parameters as stated in the problem. v = 3.2; %m/s a_t = 7.8; %m/s^2 M = 55; %kg g = 9.807; %m/s^2 alpha = 19.4 * pi /180; %radians - angle from vertical rho = 3; %meters - circular path, constant rho %Solve for the normal acceleration. a_n = v.^2 ./ rho; %m/s^2 fprintf('\nNormal Acceleration = %1.1f M/s^2 \n', a_n) %Find the normal and tangential components of the skateboarder’s weight W = M .* g; %Newtons W_n = -W .* sin(alpha); %Newtons - minus sign accounts for direction W_t = W .* cos(alpha);%Newtons fprintf('Skateboarder''s Weight = %1.0f N\n', W) fprintf('\tNormal component = %1.0f N\n', W_n) fprintf('\tTangential component = %1.0f N\n', W_t) %EQ of MOTION: normal component %N_A + W_n = M * a_n so... N_A = M .* a_n - W_n; %Newtons fprintf('Normal Force = %1.2f N\n', N_A) %EQ of MOTION: tangential component %W_t + F = M * a_t so... F = M .* a_t - W_t; %Newtons fprintf('Friction Force = %1.1f N\n', F) %Calculate the coefficient of kinetic friction mu_k = abs(F) ./ N_A; %need only the magnitude of F here fprintf('Coefficient of Kinetic Friction = %1.2f\n\n', mu_k) Principle of Work and Energy Ref: Hibbeler § 14.3, Bedford & Fowler: Dynamics § 8.1-8.2 The principle of work and energy states that the work done by all of the external forces and couples as a rigid body moves between positions 1 and 2 is equal to the change in the body’s potential energy. Hibbeler writes the resulting equation as 2211 TUT =+∑ − Hibbeler (14-7) The same equation is written as 1212 TTU −= Bedford and Fowler: Dynamics (8.4) by Bedford and Fowler. The difference between these two equations is simply nomenclature. Both ∑ −21U and 12U represent the sum of work done by all external forces and couples on the body. We can use the principle of work and energy to solve problems involving force, displacement, and velocity. Example: Skid-to-Stop Braking Distance The driver of a 1600 kg passenger car hits the brakes and skids to a stop to try to avoid hitting a deer. When the deer suddenly appeared in the headlights [at an estimated distance of 300 ft (91 m)] the car was moving at 75 mph (120 km/hr). After the incident long skid marks were visible on the pavement, and the deer was seen running over a nearby hill. Did the driver get the car stopped in time, or did the deer manage to get out of the way? Assumptions: Assume a 1.5 second reaction time – that is, it took 1.5 seconds from the time the driver saw the deer until he or she pressed the brake pedal. Assume a dry, flat road with a coefficient of friction, µk = 0.80. Solution A free body diagram of this system shows the various forces acting on the car. W NA FA 19