Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Have a look at this code snippet - aerodynamics

  1. Aug 11, 2005 #1
    Hello everybody

    This snippet works if you insert 0 as spin around the z -axis, i.e. the angular velocity around the z-axis perpendicular to the x,y-plane (or dot product to linear velocity of the ball).

    Unfortunately it all comes crashing down when I insert an angular velocity, which in return gives me a magnus force excerted on the ball due to pressure differences on the surface of the ball.

    This is the code that does the math:
    Code (Text):
            private void btnStart_Click(object sender, System.EventArgs e)
                wO("Initiating the ball object & setting parameters...");
                Ball b = new Ball(Convert.ToDouble(txtInitV.Text),
                b._X = 0;
                b._Y = 0;
                b._Z = 0;
                b.Mass = Convert.ToDouble(txtBallM.Text);
                b.Radius = Convert.ToDouble(txtBallR.Text)/1000;
                b.C_d = Convert.ToDouble(txtC_d.Text);                                  //get this experimentally!!
                b.XAngVelocity = Convert.ToDouble(txtTurnsX.Text)*2*Math.PI; //w = 2PI * n (rpm: w = [2*PI*n]/60)
                b.YAngVelocity = Convert.ToDouble(txtTurnsY.Text)*2*Math.PI;
                b.ZAngVelocity = Convert.ToDouble(txtTurnsZ.Text)*2*Math.PI; //perpendicular to normal flight path

                wO("Projection angle: " + txtInitAng.Text + " degrees");
                wO("Timespan: " + txtSimTime.Text + " milliseconds");
                wO("Initial velocity: " + Math.Round(b.GetVelocity(), 2).ToString() + " m/s");
                wO("Initial horizontal velocity: " + Math.Round(b.XVelocity, 2).ToString() + " m/s");
                wO("Initial vertical velocity: " + Math.Round(b.YVelocity, 2).ToString() + " m/s");
                wO("Ball radius: " + txtBallR.Text + " mm");
                wO("Ball mass: " + txtBallM.Text + " kg");
                wO("Z-axis spin: " + txtTurnsZ.Text + " rps");

                if (mnuAirResistance.Checked) wO("Drag is on");
                else wO("Drag is off!");

                wO("Starting simulation...");

                double end = Convert.ToDouble(txtSimTime.Text)/1000;
                double jump = 0.008; //and step

                double xOld, yOld, zOld;
                double xDrag, yDrag, zDrag; //Forces
                double xMagnus, yMagnus, zMagnus;
                double theta;

                for (double step=0; step < end; step += jump)
                    [B]xOld = b._X; //get old value
                    b._X += b.XVelocity*jump;[/B] //vt = d
                    //now give it a new velocity and calculate the horizontal component of drag.
                    [B]theta = Math.Atan(b.YVelocity/b.XVelocity);
                    xDrag = b.GetDrag()*Math.Cos(theta);
                    xMagnus = b.GetMagnusApprox()*Math.Cos(Math.PI - (theta+Math.PI/2));[/B]
                    //the drag de-acc and the magnus de-acc for 'step' amount of time
                    [B]b.XVelocity = (b._X - xOld)/jump + ((xDrag + xMagnus)/b.Mass)*jump; [/B]
                    //no magnus
    //              xDrag = b.GetDrag()*Math.Cos(Math.Atan(b.YVelocity/b.XVelocity));
    //              b.XVelocity = (b._X - xOld)/jump + (xDrag/b.Mass)*jump; //the drag de-acc for 'step' amount of time

                    [B]yOld = b._Y;
                    b._Y += b.YVelocity*jump;
                    yDrag = b.GetDrag()*Math.Sin(theta);
                    yMagnus = b.GetMagnusApprox()*Math.Sin(Math.PI - (theta+Math.PI/2));[/B]
                    //if the YVelocity is negative, then angle is negative, and drag * angle is positive.
                    //Hence a positive amount is added to the negative speed (downwards) and hence decreasing fall speed
                    //g is negative
                    [B]b.YVelocity = (b._Y - yOld)/jump + (g + ((yDrag + yMagnus)/b.Mass))*jump;[/B]
                    //no magnus
    //              yDrag = b.GetDrag()*Math.Sin(Math.Atan(b.YVelocity/b.XVelocity));
    //              b.YVelocity = (b._Y - yOld)/jump + (g + (yDrag/b.Mass))*jump;
                    writeC(Math.Round(b._X, 4).ToString() + "\t"
                        + Math.Round(b._Y, 4).ToString() + "\t"
                        + Math.Round(b._Z, 4).ToString() + "\t"
                        + Math.Round(xDrag, 4).ToString() + "\t"
                        + Math.Round(yDrag, 4).ToString() + "\t"
                        + Math.Round(b.XVelocity, 4).ToString() + "\t"
                        + Math.Round(b.YVelocity, 4).ToString());
                wO("Simulation ended");
    What I'm doing is that I'm using Newton's second (I think) law, that F=ma, and this gives me that a=F/m and since change in distance during a very short period of time is a*t I get:
    Code (Text):
    b.YVelocity = (b._Y - yOld)/jump + (g + ((yDrag + yMagnus)/b.Mass))*jump;
    that the new velocity is the change in distance over time (i.e. the regular linear velocity for "jump" seconds) plus the acceleration of gravity added together with the drag acceleration (deacc.) and magnus acceleration (can be both). These accelerations times time, gives me delta v. Hence I add delta v to v_0, initial velocity (initial for this loop revlovement at least).

    It ought to be easy, being high-school mechanics, but it won't work.
    The whole program for you to compile is attached as a zip file.

    ...and this is the code that gives the magnus force and velocity etc:
    Code (Text):
    public double Area()
                return Math.Round(this.Radius*this.Radius*Math.PI, 4);
            public double GetDrag() //Returns F_d
                return -0.5 * this.C_d * this.airDensity * Area()* this.GetVelocity();

            public double GetVelocity()
                return Math.Sqrt(   this.XVelocity*this.XVelocity +
                                    this.YVelocity*this.YVelocity +
            public double GetRenolds()
                return (airDensity*this.GetVelocity()*this.Radius*2)/airViscosity;
            public double GetRenolds(double v) //Returns constant depending on v
                return (airDensity*v*this.Radius*2)/airViscosity;
            public double GetMagnusApprox() //returns F_m
                return (Math.PI*this.airDensity*this.GetVelocity()*this.Radius*this.Radius*this.Radius*this.Radius*this.ZAngVelocity)/(2*this.Radius);
    I really would love to get this program working, and will also be glad to hear improvements on it. As stated; program attached. Thanks!


    Attached Files:

  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?