Missile Lauch friction gamma factor - Output deviating from desired output, please help!

Goal:

I'm trying to implement friction into this script of a launch of a missile.

The friction force should be related to the velocity like: Ff = -gamma|v|v.

Subsequently we also have the gravity force acting on the missile: Fg = mg.

The values of ad should be as follows, but I am unable to get this desired output:

ad = [300332, 302763, 303903, 303771, 302384, 299763, 295927, 290896, 284693];

for degrees of launch 37, 39, 41, 43, 45, 47, 49, 51, 53 respectively.

Script:

clear;
close all;
clc;

tic

m = 500;

aT = [];
ad = [];

gamma = 1.0e-5

for i = [37:2:53]
    r = 6371 * 10^3;
    G = 6.674 * 10^-11;
    M = 5.972 * 10^24;
    g = (G * M)/(r^2);
    theta0 = i;
    ax = 0;
    ay = r;
    v0 = 2000;
    vx0 = v0*cosd(theta0);
    vy0 = v0*sind(theta0);
    x = 0;
    y = r;
    vx = vx0;
    vy = vy0;
    T = 0;
    dt = 0.01;
    at = 0;
    landed = 0;
    z = 1;
    while landed == 0
        z = z + 1;
        T = T + dt;
        xo = x;
        yo = y;
        x = x + vx * dt;
        y = y + vy * dt;
        d = sqrt(x^2 + y^2);
        alpha = atand(x/y);
        g = (G*M)/(d^2);
        gy = cosd(alpha) * g;
        gx = sind(alpha) * g;
        FgY = m * gy;
        FgX = m * gx;
        FricY = gamma*abs(vy)*vy;
        FricX = gamma*abs(vx)*vx;
        FnetY = FgY - FricY;
        FnetX = FgX - FricX;
        vy = vy - (gy * dt) - (FnetY / m)*dt; % This should substract gravity and friction from the velocity, not sure if this is correct.
        vx = vx - (gx * dt) - (FnetX / m)*dt; % Same for this line of code
        v = vx/sin(alpha);
        ax = [ax, x];
        ay = [ay, y];
        if d < r
            landed = 1;
        end
    end
    
    aT = [aT, T];
    distance = (alpha/360) * 2 * pi * r;
    ad = [ad, distance];
    
    fprintf('Checked for degree: %.0f\n', i)
    
end

toc

Output:

ad = [200358, 203662, 205972, 207257, 207539, 206800, 205056, 202324, 198615];

for degrees of launch 37, 39, 41, 43, 45, 47, 49, 51, 53 respectively.

Any help / suggestions are immensly appreciated!

[the cyclist]

I'm not able to run the code right now, but one suggestion I have is to run your code without the frictional force (for which you can calculate the correct answer by hand), and see if your code is giving the correct answer for that. That way you'll know if it is the friction calculation that is wrong, or something else.

Answer 1

[Jim Riggs]

Step 1: Draw a free body diagram that shows the body with all of the forces that act on it and the coordinate frame(s) that are being used.

Step 2: For a point-mass in 2 dimensions, you should be able to use the free body diagram (by inspection) to ensure that you are calculating the force components correctly; sum of the forces in X = m*Ax, sum of forces in Y = m*Ay.

Are forces defined in body or inertial axes? (The aerodynamic force is usually in body axes, but the gravity force is defined in Inertial axes)

Show us this drawing and we can work from there.

[Jim Riggs]

One observation:

You have defined the angle alpha as atand(x/y) - This angle is based on the position of the body (x, y).

Perhaps what you want is the flight path angle, which is based on the velocity of the body (Vx, Vy);

Also, the suggestion by "the cyclst" is an excellent one: get it to work with zero friction first, then add the friction term.

[El Pepeno]

Trust me on this I really really appreciate the effort and therefor I will accept your answer.

I still don't really know how to implement this into my script tbh. My mind is just full with things right now and I can't really think clearly.

Thank you for the help though, I appreciate it.

[Jim Riggs]

For problems like this, drawing the picture always helps a lot.

It will help you define your terms, the relationship between terms, sign conventions, etc.

It will be very usefull for you to sort out your thoughts.

[El Pepeno]

The code below is the code without friction which is working. Now I am trying to add friction to this code to get the desired values of the distance as posted in my question above, but I'm struggling alot:

clear;
close all;
clc;

tic
tic
aT = 0;
ad = 0;
format long
for i = [37:2:53]
    r = 6371 * 10^3;
    G = 6.674 * 10^-11;
    M = 5.972 * 10^24;
    g = (G * M)/(r^2);
    theta0 = i;
    ax = 0;
    ay = r;
    v0 = 3000;
    vx0 = v0*cosd(theta0);
    vy0 = v0*sind(theta0);
    x = 0;
    y = r;
    vx = vx0;
    vy = vy0;
    T = 0;
    dt = 0.01;
    at = 0;
    landed = 0;
    z = 1;
    while landed == 0
        z = z + 1;
        T = T + dt;
        %at = [at, T];
        xo = x;
        yo = y;
        x = x + vx * dt;
        y = y + vy * dt;
        d = sqrt(x^2 + y^2);
        alpha = atand(x/y);
        g = (G*M)/(d^2);
        gy = cosd(alpha) * g;
        gx = sind(alpha) * g;
        vy = vy - (gy * dt);
        vx = vx - (gx * dt);
        v = vx/sin(alpha);
        ax = [ax, x];
        ay = [ay, y];
        if d < r
            landed = 1;
        end
    end
    
    aT = [aT, T];
    distance = (alpha/360) * 2 * pi * r;
    ad = [ad, distance];
    
    figure(1)
    
    th = 2.8*pi/6:pi/500:3.1*pi/6;
    xunit = r * cos(th) + 0;
    yunit = r * sin(th) + 0;
    
    figure(2)
    hold on
    plot(ax, ay)
    toc
    
    fprintf('Degree processed: %.0f\n', i)
end

th = 2.8*pi/6:pi/500:3.1*pi/6;      %use for small section of earth comparison
%th = 0:pi/500:2*pi;                %use for entire earth comparison
xunit = r * cos(th) + 0;
yunit = r * sin(th) + 0;

plot(xunit, yunit)
hold off

legend('37 degrees', '39 degrees', '41 degrees', '43 degrees', '45 degrees', '47 degrees', '49 degrees', '51 degrees', '53 degrees', 'earth')
title('The trajectory of the different missiles')
xlabel('x(m)')
ylabel('y(m)')
toc

[Jim Riggs]

The friction force always opposes the velocity vector.

To add friction, calculate the total friction force using your physical model for the friction.

Then, apply this to the body such that it is in the opposite direction to the velocity.

If F is the total force on the body, calculate the unit velocity vector Uvx = Vx / |V|, Uvy = Vy / |V|,

now the friction force is applied using Fx = -F*Uvx, Fy = -F * Uvy.

[El Pepeno]

What is the total force in this case though and what direction does it have?

[Jim Riggs]

I am not sure because I do not know the details of the problem.

If the friction you want to consider is the aerodynamic drag in a projectile moving through the atmosphere, then this type of force can be characterized in aerodynamic terms as

DragForce = 1/2 * AirDensity * Airspeed^2 * DragCoefficient * referenceArea.

The reference area is a constant (depending on the geometry of the body) And for a simple model, you may consider the drag coefficient as a constant too. The details of the aerodynamic force are then a function of the air density. So you need some kind of model for air density (typically a function of the altitude above sea level of the projectile).

For a more sophisticated model, the drag coefficient could be a function of the velocity as well. You would need more information on the aerodynamics of your projectile to include this detail.

As for the direction, it is opposite the velocity vector;

Fx = -DragForce * Uvx

Fy = -DragForce * Uvy

For the aerodynamic force, the velocity used is the velocity of the projectile relative to te air. If there is no wind, this is the same as the velocity of the projectile relative to the earth.

[El Pepeno]

friction is determined by the following relation: Ff = -gamma * | v | * v

but I have no idea how to implement this into the script

only gravity and friction should be acting on the missile, no thrust is present, only an initial starting velocity of 2000 m/s

[El Pepeno]

but do I plug these Fx and Fy into the velocities vx and vy? because I also need there to be a gamma factor in there somewhere.. (and to convert Fx to velocity I should do vx = vx - (gx * dt) - (Fx / m)dt right?, but then still theres no gamma)

[Jim Riggs]

OK. If you compare your definition of the Friction force to my equation for the aerodynamic drag, you get

-gamma * |v| * v = 1/2 AirDensity * Airspeed^2 * DragCoefficient * referenceArea

The |v| * V term is a way of computing Airspeed^2 and preserving the sign, so it can be positive or negative. If we let |v| * v = Airspeed^2 in my equation, we are left with

-gamma = 1/2 AirDensity * DragCoefficient * referenceArea.

So gamma represents an aerodynamic factor which is good for a constant air density, and a constant drag coefficient (the reference area is always constant).

To apply this, you are almost there, with one small change.

The drag force in the X direction is DragForce * Uvx

Note that Uvx = Vx / | V | and V^2 = |V | * | V |, so

If DragForce = 1/2 * AirDensity * | v | * | v | * DragCoefficient * referenceArea
Fx = 1/2 AirDensity  | v| * | v| * DragCoefficient * referenceArea * Vx / | v|
substitute -gamma = 1/2 AirDensity * DragCoefficient * referenceArea
Fx = -gamma * | v | * | v | * Vx / | v |
Fx = -gamma * | v | * Vx

Compare this with your equation

FrictX = gamma * abs(vx) * vx;

So, you need to make 2 changes:

1) You need a - sign on gamma

2) the abs(Vx) should be abs(V) the magnitude of the total velocity vector.

[El Pepeno]

Thank you so much for the help, I will try to figure it out now.

Much appreciated, have a nice week!