Re-entry Projectile Analysis Program

1. Apr 23, 2015

justin cooper

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

//System Variables
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 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);
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;
// theta = atan(velocityY/velocityX)*180/PI;
++timecount;

//Output Results to CSV Format

if (myfile.is_open()){
}
}
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. Apr 23, 2015

Futurestar33

What kind of class is this for just curious?

3. Apr 23, 2015

justin cooper

This is for a capstone project for a senior design.

4. Apr 24, 2015

kroni

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.

5. Apr 24, 2015

justin cooper

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:

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

6. Apr 24, 2015

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

7. Apr 24, 2015

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 26, 2015

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.