Sign in to follow this  

Shuttlecock trajectory

This topic is 4351 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

Hello, I would like to simulate a badminton shuttlecock trajectory. I have found a phD thesis about this and i got the following Equation of Motion : m PLUS Dcos(teta) PLUS Lsin(teta) = 0 mÿ PLUS Dsin(teta) MINUS Lcos(teta) PLUS mg = 0 WHERE L Lift force and L = 1/2(ro*v²*S*Cd) D Drag force and D = 1/2(ro*v²*S*Cl) teta Angle between velocity vector and horizontal m Shuttlecock mass y Vertical distance ro Air density S Square of maximum shuttlecock skirt diameter Cd Aerodynamic drag coefficient Cl Aerodynamic lift coefficient v Shuttlecock velocity How can I get coordinate values(x,y) of the shuttlecock over the time ? I know it's about differential equation, but because of my poor mathematics knowledges I can't find the solution. Maybe somebody can guide me o how to proceed ? Thank you Fransoa

Share this post


Link to post
Share on other sites
Well, fortunately for you your solution is computational.

You'll probably use discrete increments of t. You'll perform the following steps:
1. Calculate forces
2. Apply forces to velocity
3. Apply velocity to position

#3 is the easiest. You'll have three equations:

x += vx * t;
y += vy * t;
z += vz * t;

#2 isn't too hard either:

vx += fx/m * t;
vy += fy/m * t;
vz += fz/m * t;

#1 is the one that confuses most people:
fx = fy = fz = 0;
fx = Cd * vx * vx;
fz = Cd * vz * vz;
fy = m * g + (Cl + Cd) * vy * vy;

where
x,y,z are positions
vx,vy,vz are velocities
fx,fy,fz are forces
m is the mass
g is gravitational acceleration. If your units are in meters, this is -9.8
Cd is the drag constant. It is the air density * area of the shuttlecock * some tuning factor / 2. It should be less than 0.
Cl is the lift constant. It is the same thing as Cd, except positive, and with a different tuning factor. It should be very close to abs(Cd).

x,y,z will start at 0.
vy will start > 0
vx,vz will be whatever you like.

Typically, Cd and Cl should be small.

Share this post


Link to post
Share on other sites
First, thank you for your very fast answer.

I am trying to understand how you can get to that result, sorry to bother you...

Could you give me some more explanations about :

"#1 is the one that confuses most people:
fx = fy = fz = 0;
fx = Cd * vx * vx;
fz = Cd * vz * vz;
fy = m * g + (Cl + Cd) * vy * vy;"

Why does Cl appears only for fy ?
In my equation there is teta, because incidence of the shuttle has to be taken into account, is it included in your "some tuning factor" ?

More I also have an equation for teta....

Thanks
Fransoa


Share this post


Link to post
Share on other sites
Well, Cl is the coefficient of lift, and it only really affects your height (therefore y) in this case. It is really the same principle as drag, which is why the equations are so similiar.

Theta is the angle of inclination from the ground at which the shuttlecock is launched. This would be incorporated into the starting values for vx,vy, and vz. You could decide the original vx,vy,vz by setting
vx=v*cos(theta)*cos(phi),
vy=v*sin(theta)*cos(phi),
vz=v*sin(phi)

or in 2 dimensions
vx=v*cos(theta)
vy=v*sin(theta)

Where v is the total velocity

Share this post


Link to post
Share on other sites
Sorry but i am a little bit confused.....reading your advice and reading the thesis. You can find it here : http://www.thnwb.com/cookeassociates.com/thesis.zip

Could you have a look at the Chapter 3 ???

Are you ok to start from the beginning ?
How is it possible to get the first equation ?

m + Dcos(p) + Lsin(p) = 0

L Lift force and L = 1/2(ro*v²*S*Cd)
D Drag force and D = 1/2(ro*v²*S*Cl)
m mass
p angle of inclination

Could you explain me how to get the following one ?
mÿ +Dsin(p) +Lcos(p) +mg = 0

Thanks
Fransoa

Share this post


Link to post
Share on other sites
These two equations, if you notice at the bottom, are for steady state. That means that no additional forces are acting on the shuttlecock, and the the forces that are acting on the shuttlecock sum to 0. Therefore, the shuttlecock continues at its unaltered velocity. This isn't necessarily what you want, but it gives you a good starting point.

The equations start with a tallying of forces:
my:, where y: is the second derivated of y wrt t. It is equivalent to accleration a.

D sin(theta) is the drag force, which they take to only act in the forward direction. Theta is the angle of inclination, and because this is a sine factor, then this force is maximize when motion is all forward, and 0 when the motion is straight up.

L cos(theta) is the lift force, which acts in the up direction. It is really the same thing as drag, only it forces the shuttlecock up due to its shape. Because it is a cosine, then lift is maximize when motion is upward, and zero when motion is forward.

mg is the force of gravity.

The equation is set to 0 to find a stable condition.

The first equation can be found by integrating this equation several times to find appropriate starting conditions.

The equation of motion you want is this:
my - Dsin@ + Lcos@ = 0
y = (D/m)sin@ - (L/m)cos@

[Edited by - erissian on January 11, 2006 5:09:02 PM]

Share this post


Link to post
Share on other sites
Ok I better understand. But with this equation I can only get y. How can i make x and t appear in an equation from y = (D/m)sin@ - (L/m)cos@ ???

I know it's basic question, but I need to understand....

Also, if you look at the chapter 4.2 The Trajectory Prediction, it says they created a computer program to plot the trajectory, using Runge-Kutta algorithm based on those equations i told you. Why are they doing this if those equations are wrong ? So what are they doing exactly ?

And also there is an angular equation of motion, do you think it could be useful or should i use the velocity vector to set the orientation of the shuttlecock ?

Thanks for the course
Fransoa

Share this post


Link to post
Share on other sites
Well, I would treat this with vectors.
I would sum the forces:
fy = -m*g - D*vy*vy + L*vy*vy;
fx = -D*vx*vx; // No gravity or lift

Then I would adjust the velocities:
vy += fy * t / m;
vx += fx * t / m;

Then I would adjust the positions:
y += vy * t;
x += vx * t;

A good number for these would be around
D = 0.5 * air_density * area * coefficient_of_drag
= 0.5 * 1.2 * 0.0025 * 0.48;
= 7e-4;
L should be about 70% of that, so
L = 5e-4;

So you could rewrite fy as:
fy /*= (L - D)*vy*vy - mg */ = -2e-4*vy*vy - mg;
fx = -7e-4*vx*vx;

Share this post


Link to post
Share on other sites
Ok I have used this method and the shape of the trajectory I plot looks good, however for a high clear starting at 47 m.s-1, the birdie comes to y=0 for x=66 meters, isn't that strange ???

Does that come from the drag force that should be higher ??

Here are the values I have used :

Cd 0,48
Cl 0,35
D 0,001331712
L 0,00097104
t 0,005 s
@ 0,628318531 rads
v 47 (m.s-1)
m 0,049 (kg)
g 9,8
x0 0
y0 3

Thanks
Fransoa

Share this post


Link to post
Share on other sites
Well, you can always solve it without drag and lift to get a sanity check.

Using x = x0 + v0t + at2

Solving for t using y:
0 = 3 + 47 * sin(0,628318531) * t - 9.8 * t2
0 = -9.8t2 + 27.63t + 3

Using the quadratic equation:
t = [-27.63 - sqrt(27.632 + 9.8 * 4 * 3)]/[-9.8*2]
t = 2.92 s

This is the time when the shuttlecock lands. To find out how far it goes:
x = x0 + v0t + at2
x = 0 + 47 * cos(0,628318531) * (2.92) + 0 * 2.922
x = 111 m

So with those starting conditions, neglecting air effects, the shuttlecock would travel 111 meters. So 66 meters with air resistance seems reasonable after all.

[Edited by - erissian on January 12, 2006 3:15:29 PM]

Share this post


Link to post
Share on other sites
Well, the math doesn't lie. If it seems unrealistic, then you might consider launching the shuttlecock at less than 100 miles per hour (47 m/s).

You can work backwards, decide what a realistic range is and then work out initial velocity from there.

Share this post


Link to post
Share on other sites
Ok. But I didn't invent that the shuttlecock is launched at that speed, it has been measured, so this initial speed is realistic. I believe that maths don't lie. So where is the error ? :)

Fransoa

Share this post


Link to post
Share on other sites
Well, what would be an acceptable maximum range for a shuttlecock? I read an article found on google that a shuttlecock can frequently reach up to around 150 mph (70 m/s), so that speed is ok. Then again, I also read that a shuttlecock can go as far as 150 yards (140 m), so maybe our numbers aren't too far off.

The range should maximize at 45 degrees, discounting air effects, so perhaps a higher angle is better? I don't play the game myself. If you want to limit the motions, then decide what your constraints are, and we'll work out how to stay within them.

Share this post


Link to post
Share on other sites
Ok.

Speed : from 0 to 70 m/s in the case of power smash or a power clear (out of the field).

Angle : there is no limitation, I would say from -85 to +85, because when you are very close to the net you can play a net shot(+85 max) or a net shot attack(-85). But If you consider you can miss the shot I would say -90/+90....

distance : the field is 13.4 meters long. If the powerful player in the world hit the shuttle as fast as he can from x=0, I would say he can reach 15.5m (far from 140... :) where did you read that ?!! I would like to see the article)

Let's have a look at the chapter 4 of the thesis : 4.3.3.5,there are some examples for different kind of shot. One example is :
v0=27.25m/s, @ = 20°,y0 = 1, Range = 8.5m, max height = 2.35m

I am all listening now ;)

Fransoa

Share this post


Link to post
Share on other sites
Turns out 150 yards was from a shuttlecock style bullet :)

Anyways, I think the problem is in the implementation somewhere. I wrote a quick program and I got a range of about 5 meters for the values you said you were using. Could you post your code?

Here's the heart of what I used:


while ( y >= 0 ) {
vx += -drag * vx * abs(vx) * timestep / mass;
vy += -g * timestep + (lift - drag) * vy* abs(vy) * timestep / mass;
x += vx * timestep;
y += vy * timestep;
time += timestep;
printf("%0.1fs : %0.2f, %0.2f m : %0.2f, %0.2f m/s\n",time,x,y,vx,vy);


[Edited by - erissian on January 13, 2006 1:12:37 PM]

Share this post


Link to post
Share on other sites
I used an excel spreadsheet to plot the trajectory....
But my code would have been :


while ( y >= 0 ) {
vx += -drag * vx * abs(vx) * timestep / mass;
vy += -g * timestep + (lift - drag) * vy* abs(vy) * timestep / mass;
x += vx * timestep;
y += vy * timestep;
}


So I misunderstood something about time !!!

Quote:

Well, I would treat this with vectors.
I would sum the forces:
fy = -m*g - D*vy*vy + L*vy*vy;
fx = -D*vx*vx; // No gravity or lift

Then I would adjust the velocities:
vy += fy * t / m;
vx += fx * t / m;

Then I would adjust the positions:
y += vy * t;
x += vx * t;


I thought t was always timestep. Why do you use time to compute vx and not timestep as for vy ?


Share this post


Link to post
Share on other sites
Umm, call it an error in implementation :)

I fixed the code in that post, and now my result is around 14 meters.

t is 'timestep', and I'm using 'time' as a total running time. I didn't label either one 't' because I thought it would be too ambiguous.

Share this post


Link to post
Share on other sites
mmmh, that is very strange.... I have checked my formulas in the excel spreadsheet and it looks the same as your code. However I don't get 14 meters(it would be a pretty good result).
So what I suppose is that we are not using the same values :

Cd = 0,48
Cl = 0,35
D = 0.5 * 1.2 * 0.068 * 0.068 * Cd = 0,001331712
L = 0.5 * 1.2 * 0.068 * 0.068 * Cl = 0,00097104
t = 0,005
@ = 0,628318531
v = 47
m = 0,049
g = 9,8


If I use D = 0,01331712 instead of 0,001331712 i get a result around 12m....

Fransoa

Share this post


Link to post
Share on other sites
Well, here's the whole thing:


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc,char* argv[])
{
float timestep = 0.005; // s
float velocity = 47; // m/s
float inclination = 36; // degrees
float x = 0; // meters
float y = 3; // meters
float lift = 0.00097104;
float drag = 0.001331712;
float mass = 0.0049;
float g = 9.8;
float time = 0;

float vy = velocity * sin(inclination * M_PI / 180);
float vx = velocity * cos(inclination * M_PI / 180);

printf("t: x,y vx,vy\n");

while ( y >= 0 ) {
vx += -drag * vx * abs(vx) * timestep / mass;
vy += -g * timestep + (lift - drag) * vy* abs(vy) * timestep / mass;
x += vx * timestep;
y += vy * timestep;
time += timestep;
printf("%0.1fs : %0.2f, %0.2f m : %0.2f, %0.2f m/s\n",time,x,y,vx,vy);
}

return 0;
}

Share this post


Link to post
Share on other sites
Ok sorry my mistake, i used 0.049 instead of 0.0049 ... shame on me !

However the shape doesn't look so good. The shuttle goes too high, more than 10 meters...But in distance that's ok !

I have read in the thesis that there are different values for Cl :
0.07 @ 10°
0.21 @ 20 °
0.35 @ 30°

Do you think I should adjust Cl all along the motion ?

And also, What is the difference between your approach and this one :
Quote:

4.2.1 The Main Program

The format of the main program can be seen in Figure 4.19. It reads in the necessary data, some of which is input directly by the user. The rest is read from data files which contain the experimental design parameters measured as described in Chapters 2 and 3. The main program calls the necessary subroutines for the trajectory prediction and writes the resulting trajectory data into a file, i.e. the time, horizontal and vertical distances, resultant and component velocities, incidences and direction of the velocity vector.

4.2.2 Subroutine for 2D Particular Motion; TRAJEC

The flow chart shown in Figure 4.20 describes the subroutine program which is used for predicting the shuttlecock trajectory with and without incidence. (For particular motion, with no incidence, TRAJEC is similar but does not call the 2D angular subroutine, TURN).

4.2.2.1 Interpolation of the Drag Coefficient, Cd

Two standard subroutines, SPLINE and SPLINT, from Press et al [27] were used to carry out a cubic spline interpolation method for setting the value of the drag coefficient Cd at a particular velocity. This was required as the coefficient varies with Reynolds number and therefore varies throughout the trajectory. (This is described in Chapter 2.)
• SPLINE: given the two arrays of data for the drag parameter at the corresponding velocity and the first derivatives of the interpolating function at the first and last data points, this subroutine computes an array of the same length containing the second derivatives of the interpolating function at the tabulated points.
• SPLINT: given the two arrays of data for the drag parameter at the corresponding velocity and the array of second derivatives computed by SPLINE, this subroutine will return a value of the drag parameter at a given velocity.

4.2.2.2 Integration of the Differential Equations of Motion

Two standard Fortran subroutines from Press et al [27] were used to carry out a fourth-order Runge-Kutta integration on the equations of motion:
• DERIVS4: given the horizontal and vertical velocities at a certain time this subroutine computes the accelerations at that time. In this case the equations that were used were the shuttlecock equations of motion, Equations 3.1 and 3.2.
• RK44: given the horizontal and vertical accelerations of the shuttlecock with the corresponding velocities at a certain time this subroutine uses the fourth-order Runge-Kutta method to advance the solution over a given time step and to compute the incremented velocities at the new time.
A suitable time step size for the Runge-Kutta method was found to ensure an accurate prediction with no divergence. It was found in the following way. A theoretical value for the range of a shuttlecock in a vacuum given certain initial conditions was calculated. This assumed no aerodynamic forces. By varying the time step, the optimum value was chosen such that the computed range agreed with the theoretical range as accurately as possible. This is shown in Tables 4.2 and 4.3. The optimum time step for the computer program proved to be 0.005 seconds.


Umm, i know it's a long post. Just tell me if you want to stop that thread..

thx
Fransoa

Share this post


Link to post
Share on other sites
Since this isn't thesis work, I don't do much of that at all :)

This situation is meant to be more of an introduction to this kind of situation.

Also, I've added in this bit, and I'm calculating max range and max height from inclinations of 30 - 45 degrees. Those are the values between which range should be the greatest.


float air_density = 1.225; // kg/m^3
float skirt_radius = 0.02; // m
float Cl = 0.2;
float Cd = 1.0;

float lift = Cl * M_PI * air_density * skirt_radius * skirt_radius;
float drag = Cd * M_PI * air_density * skirt_radius * skirt_radius;


I've found that these yield much more reasonable ranges. (11,7)

Share this post


Link to post
Share on other sites

This topic is 4351 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this