The document Solving Equations of Motion for Systems of Particles with MATLAB (Part -2) Civil Engineering (CE) Notes | EduRev is a part of the Civil Engineering (CE) Course Introduction to Dynamics and Vibrations- Notes, Videos, MCQs.

All you need of Civil Engineering (CE) at this link: Civil Engineering (CE)

**Example 3: Silly FE exam problem**

This example shows how polar coordinates can be used to analyze motion.

The rod shown in the picture rotates at constant angular speed in the horizontal plane. The interface between block and rod has friction coefficient Î¼ . The rod pushes a block of mass m, which starts at r=0 with radial speed V. Find an expression for r(t).

**1. Introduce variables to describe the motion **â€“ the polar coordinates r,Î¸ work for this problem

**2. Write down the position vector and differentiate to find acceleration** â€“ we donâ€™t need to do this â€“ we can just write down the standard result for polar coordinates

**3. Draw a free body diagram** â€“ shown in the figure â€“ note that it is important to draw the friction force in the correct direction. The block will slide radially outwards, and friction opposes the slip.

**4. Write down Newtonâ€™s law**

**5. Eliminate reactions**

**F=ma** gives two equations for N and T. A third one comes from the friction law T = Î¼N

The third solution can be rearranged into an equation of motion for r

**6. Identify initial conditions Here**, r=0 dr/dt=V at time t=0.

**7. Solve the equation:** If youâ€™ve taken AM33 you will know how to solve this equationâ€¦ But if not, or you are lazy, you can use MAPLE to solve it for you

This can be simplified slightly by hand

**Example 4: Motion of a pendulum (R-rated version)**

A pendulum is a ubiquitous engineering system. You are, of course, familiar with how a pendulum can be used to measure time. But itâ€™s used for a variety of other scientific applications. For example, Professor Criscoâ€™s lab uses pendulum to measure properties of human joints, see Crisco Joseph J; Fujita Lindsey; Spenciner David B, â€˜The dynamic flexion/extension properties of the lumbar spine in vitro using a novel pendulum system.â€™ Journal of biomechanics 2007;40(12):2767-73

In this example, we will work through the basic problem of deriving and solving the equations of motion for a pendulum, neglecting air resistance.

**1. Introduce variables to describe motion: **The angle Î¸(t) shown in the figure is a convenient variable.

**2. Write down the position vector as a function of the variables **We introduce a Cartesian coordinate system with origin at O, as shown in the picture.

Elementary geometry shows that r = l sinÎ¸ i-l cosÎ¸ j

Note that we have assumed that the cable remains straight â€“ this will be true as long as the internal force in the cable is tensile. If calculations predict that the internal force is compressive, this assumption is wrong. But there is no way to check the assumption at this point so we simply proceed, and check the answer at the end

**3. Differentiate the position vector to find the acceleration: **The computer makes this painless.

4. The free body diagram is shown in the figure. The force exerted by the cable on the particle is introduced as an unknown reaction force. The force vector is F = -R sinÎ¸ i + (R cosÎ¸ - mg) j

5. Newtonâ€™s laws of motion can be expressed as** F = ma **

Equating the i and j components gives two equations for the two unknowns

**6. Eliminate the reaction forces.** â€“ In this problem, it is helpful to eliminate the unknown reaction force R. You can do this on the computer if you like, but in this case it is simpler to do this by hand. You can simply multiply the first equation by cosÎ¸ and the second equation by sinÎ¸ and then add them. This yields

**7. Identify initial conditions**. Some calculations are necessary to determine the initial conditions in this problem. We are given that Î¸ = 0 at time t=0, and the horizontal velocity is V_{0} at time t=0, but to solve the equation of motion, we need the value of dÎ¸ / dt = 0 . We can find the relationship we need by differentiating the position vector to find the velocity

Setting v = v_{0}i and Î¸ = 0 at t=0 shows that

so dÎ¸ / dt = V_{0} /l

**8. Solve the equations of motion** This equation of motion is too difficult for MAPLE but actually the solution does exist and is very well known â€“ this is a classic problem in mathematical physics. With initial conditions Î¸ = 0, dÎ¸ / dt = V_{0} / l t = 0 the solution is

The first solution describes swinging motion of the pendulum, while the second solution describes the motion that occurs if you push the pendulum so hard that it whirls around on the pivot. The equations may look scary, but you can simply use MAPLE to calculate and plot them.

1. In the first equation, sn(x,k) is a special function called the `sin amplitude.â€™ You can think of it as a sort of trig function for adults â€“ in fact for k=0, sn(x,0) = sin(x) and we recover the PG version. You can compute it in Mupad using jacobiSN (x,k)

2. Similarly, am(x,k) is a function called the `Amplitude.â€™ You can calculate it in Mupad using jacobiAM (x,k). In Mupad, the am function has range -Ï€<am(x,k)<Ï€ , so the solution predicts that as the pendulum whirls around the pivot, the angle Î¸ increases from 0 to 2Ï€ , then jumps to -2Ï€ , increases to 2Ï€ again, and so on.

You might have solved the pendulum problem already in an elementary physics course, and might remember a different solution. This is because you probably only derived an approximate solution, by assuming that the angle Î¸ remains small. This occurs when the initial velocity satisfies which case the solution can be approximated by

**3.3.3 Numerical solutions to equations of motion using MATLAB**

In the preceding section, we were able to solve all our equations of motion exactly, and hence to find formulas that describe the motion of the system. This should give you a warm and fuzzy feeling â€“ it appears that with very little work, you can predict everything about the motion of the system. You may even have visions of running a consulting business from your yacht in the Caribbean, with nothing more than your chef, your masseur (or masseuse) and a laptop with a copy of MAPLE.

Unfortunately real life is not so simple. Equations of motion for most engineering systems cannot be solved exactly. Even very simple problems, such as calculating the effects of air resistance on the trajectory of a particle, cannot be solved exactly.

For nearly all practical problems, the equations of motion need to be solved numerically, by using a computer to calculate values for the position, velocity and acceleration of the system as functions of time. Vast numbers of computer programs have been written for this purpose â€“ some focus on very specialized applications, such as calculating orbits for spacecraft (STK); calculating motion of atoms in a material (LAMMPS); solving fluid flow problems (e.g. fluent, CFDRC); or analyzing deformation in solids (e.g. ABAQUS, ANSYS, NASTRAN, DYNA); others are more general purpose equation solving programs.

In this course we will use MATLAB, which is widely used in all engineering applications. You should complete the MATLAB tutorial before proceeding any further.

In the remainder of this section, we provide a number of examples that illustrate how MATLAB can be used to solve dynamics problems. Each example illustrates one or more important technique for setting up or solving equations of motion.

**Example 1: Trajectory of a particle near the earthâ€™s surface (with air resistance)**

As a simple example we set up MATLAB to solve the particle trajectory problem discussed in the preceding section. We will extend the calculation to account for the effects of air resistance, however. We will assume that our projectile is spherical, with diameter D, and we will assume that there is no wind. You may find it helpful to review the discussion of aerodynamic drag forces in Section 2.1.7 before proceeding with this example.

**1. Introduce variables to describe the motion:** We can simply use the Cartesian coordinates of the particle (x(t), y(t), z(t))

**2. Write down the position vector in terms of these variables: **r = x(t)i + y(t) j + z(t)k

**3. Differentiate the position vector with respect to time to find the acceleration. Simple calculus gives**

**4. Draw a free body diagram.** The particle is now subjected to two forces, as shown in the picture. Gravity â€“ as always we have Fg = -mgk.

Air resistance.

The magnitude of the air drag force is given by , where

- p is the air density,
- C
_{D}is the drag coefficient, - D is the projectileâ€™s diameter, and
- V is the magnitude of the projectileâ€™s velocity relative to the air. Since we assumed the air is stationary, V is simply the magnitude of the particleâ€™s velocity, i.e.

The Direction of the air drag force is always opposite to the direction of motion of the projectile through the air. In this case the air is stationary, so the drag force is simply opposite to the direction of the particleâ€™s velocity. Note that v /V is a unit vector parallel to the particleâ€™s velocity. The drag force vector is therefore

The total force vector is therefore

5. Write down Newtonâ€™s laws of motion.

It is helpful to simplify the equation by defining a specific drag coefficient , so that

The vector equation actually represents three separate differential equations of motion

**6. Eliminate reactions** â€“ this is not needed in this example.

**7. Identify initial conditions.** The initial conditions were given in this problem - we have that

**8. Solve the equations of motion.** We canâ€™t use the magic â€˜dsolveâ€™ command in MAPLE to solve this equation â€“ it has no known exact solution. So instead, we use MATLAB to generate a numerical solution.

This takes two steps. First, we must turn the equations of motion into a form that MATLAB can use. This means we must convert the equations into first-order vector valued differential equation of the general form dt/dt = f(t,y). Then, we must write a MATLAB script to integrate the equations of motion.

Converting the equations of motion: We canâ€™t solve directly for (x,y,z), because these variables get differentiated more than once with respect to time. To fix this, we introduce the time derivatives of (x,y,z) as new unknown variables. In other words, we will solve for (x, y, z, v_{x} , v_{y} , v_{z}) , where

These definitions are three new equations of motion relating our unknown variables. In addition, we can re-write our original equations of motion as

So, expressed as a vector valued differential equation, our equations of motion are

MATLAB script. The procedure for solving these equations is discussed in the MATLAB tutorial. A basic MATLAB script is listed below.

This produces a plot that looks like this (the plotâ€™s been edited to add the grid,etc)

Offer running on EduRev: __Apply code STAYHOME200__ to get INR 200 off on our premium plan EduRev Infinity!