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

Re-entry Projectile Analysis Program

  1. Apr 23, 2015 #1
    Hi all:

    I'm writing a program based on simple re-entry physics and can't quite get the algorithm to match the theoretical results. A little background:

    The equations I'm comparing this to are derivations from the equations of motion and the Sutton-Graves equation for peak heat flux at the stagnation point. The main body of the program is:

    // ****** KRUPS RE-ENTRY ANALYSIS OF MOTION ***** //
    // ********************************************** //
    //
    #include <fstream>
    #include <iostream>
    #include <sstream>
    #include <stdlib.h>
    #include <math.h>
    using namespace std;

    #define PI 3.14159265

    int main(){

    //User Input
    double reVel = 7500; //m/s, re-entry velocity
    double dragCo = 1.0625; //drag coefficient
    double vehicleMass = 10; //kg
    double cross = 0.06131; //m^2, cross-sectional area
    double positionY = 120000;
    double positionX = 0; //m
    double theta = 1.5; //degrees, re-entry angle
    double noseRadius = 0.06985; //m

    //System Variables
    double earthRadius = 6371000; //Earth mean radius
    double standardG = 9.806;
    double earthGravity; //m/s^2
    double decelX; //m/s^2, deceleration due to drag
    double decelY;
    double velocityX; //m/s
    double velocityY;
    double time = 0.01; //seconds, time increment
    int timecount = 1; //Number of iterations
    double rho; //kg/m^3, atmospheric density
    double totalVelocity;
    double totalAccel;
    double heatFlux; //W/m^2
    double heatLoad = 0; //J/m^2
    double earthG = 0;
    double dragX;
    double dragY;

    //Open File
    ofstream myfile("analysis.csv");
    if (myfile.is_open()) myfile<<"Time Passed,X Position,Y Position,X Velocity,Y Velocity,Total Velocity,Deceleration X,Deceleration Y,Total Deceleration,Earth G's,Angle,Atmospheric Density,Heat Flux,Heat Load"<<endl;

    velocityX = cos(theta*PI/180)*reVel;
    velocityY = sin(theta*PI/180)*reVel;
    while (positionY>0){
    //Calculate Velocity Components from Initial Velocity
    rho = 1.225*exp(-0.0001378*positionY);
    earthGravity = standardG*pow((earthRadius/(earthRadius+positionY)),2);
    dragX = cos(theta*PI/180)*rho*0.5*pow(velocityX,2)*dragCo*cross/vehicleMass;
    dragY = sin(theta*PI/180)*rho*0.5*pow(velocityY,2)*dragCo*cross/vehicleMass;
    decelX = dragX;
    decelY = dragY-earthGravity;

    // cout<<"Position X = "<<positionX<<endl;
    // cout<<"Position Y = "<<positionY<<endl;
    // cout<<"Velocity X = "<<velocityX<<endl;
    // cout<<"Velocity Y = "<<velocityY<<endl;
    // cout<<"Drag X = "<<dragX<<endl;
    // cout<<"Drag Y = "<<dragY<<endl;
    // cout<<"Decel X = "<<decelX<<endl;
    // cout<<"Decel Y = "<<decelY<<endl;
    // cout<<"Angle = "<<theta<<endl;
    // cout<<"Density = "<<rho<<endl;
    // cin.get();

    positionX = positionX + velocityX*time - 0.5*decelX*pow(time,2);
    positionY = positionY - velocityY*time - 0.5*decelY*pow(time,2);
    velocityX = velocityX - decelX*time;
    velocityY = velocityY - decelY*time;
    totalVelocity = pow(pow(velocityX,2)+pow(velocityY,2),0.5);
    totalAccel = pow(pow(decelX,2)+pow(decelY,2),0.5);
    earthG = totalAccel/earthGravity;
    heatFlux = (0.000183*pow(rho/noseRadius,0.5)*pow(totalVelocity,3));
    heatLoad += heatFlux*time/timecount;
    // theta = atan(velocityY/velocityX)*180/PI;
    ++timecount;

    //Output Results to CSV Format

    if (myfile.is_open()){
    myfile<<time*timecount<<","<<positionX<<","<<positionY<<","<<velocityX<<","<<velocityY<<","<<totalVelocity<<","<<decelX<<","<<decelY<<","<<totalAccel<<","<<earthG<<","<<theta<<","<<rho<<","<<heatFlux/10000<<","<<heatLoad<<endl;
    }
    }
    myfile.close();
    return 0;
    }

    The program seems to be over-predicting the heat flux and the deceleration.
    The theoretical values are:

    -3.8g's MAX from a_max = 7500^2*0.0001378*sin(1.5)/2e = 37.32m/s^2 = 3.8g

    q_max = 1.833E-4*(6345m/s)^3*SQRT(0.023787) = 393.9714W/cm^2

    I'm having a problem finding the bug, but my output for the vehicle I'm analyzing looks to be about 10 times higher for deceleration and about 200 W/cm^2 higher than the theoretical.

    I do KNOW that the assumptions will cause a much different profile since the angle changing will increase the deceleration considerably, but I can't be certain if it's a bug or not even though I've gone through the code a hundred times. Also I should mention that I know that the problem has to be within the equations of motion because I'm using the Sutton-Graves within the program to predict stagnation heating points. So if the velocity is off, then the q will be off as well.

    I even tried to remove the change of theta and the gravity term from the equations of motion but still do not get matching numbers?

    If anyone can help it would be greatly appreciated.

    Thanks
    JMC
     
  2. jcsd
  3. Apr 23, 2015 #2
    What kind of class is this for just curious?
     
  4. Apr 23, 2015 #3
    Thanks for the quick reply.

    This is for a capstone project for a senior design.
     
  5. Apr 24, 2015 #4
    Ok for the debug you need to analyse each intermediate variable and compute value by hand. I see this "heatLoad += heatFlux*time/timecount;" Be careful because i don't remember but i think that dividing a double by an integer give a integer value. You need to do this : heatLoad += heatFlux*time/((double) timecount))".

    For you C++ Skill :
    Pi is in math.h as M_PI.
    pow(a,2) is very slow, prefer a*a
    pow(a,0.5) is very slow, prefer sqrt()
    Don't check file each time by adding a else to the first check
    if(!myfile.open()) etc
    else
    - throw exeception
    - abort(),
    - return code;

    good luck for the bug.
     
  6. Apr 24, 2015 #5
    Thank you for the all the suggestions! I have iterated hand values, unsuccessfully. I'm a new programmer to c++ so thank you for the input. The Heat LOAD is wrong, as it was just changed to:

    heatLoad += heatFlux*time;

    Since I am after a graph of heatFlux as it progresses through time. The Total Heat Flux can be summed from the outputs in MATLAB.

    Thanks again
    JMC
     
  7. Apr 24, 2015 #6

    Mark44

    Staff: Mentor

    No. Dividing a double by an int (or short, long, char, long long, etc.) causes the int value to be promoted to double, and then floating point division is performed. It's only when both operands are an integral type that integer division is performed.
    I agree completely. I never use pow() to compute the square of a number
     
  8. Apr 24, 2015 #7

    Mark44

    Staff: Mentor

    No. Dividing a double by an int (or short, long, char, long long, etc.) causes the int value to be promoted to double, and then floating point division is performed. It's only when both operands are an integral type that integer division is performed.
    I agree completely. I never use pow() to compute the square of a number
     
  9. Apr 26, 2015 #7

    mfb

    User Avatar
    2016 Award

    Staff: Mentor

    You seem to take the angle into account twice here, first in theta and then again in the x^2 and y^2 expressions.
    Also, the angle will change over time.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Re-entry Projectile Analysis Program
  1. Explain this program. (Replies: 3)

  2. Programming languages (Replies: 16)

  3. HMI programming (Replies: 1)

Loading...