Re-entry Projectile Analysis Program

In summary: 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. 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*aI agree completely. I never use pow() to compute the square of a number pow(a,0.5) is very slow, prefer sqrt()Don't check file each time by adding a else to the first
  • #1
justin cooper
3
0
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
 
Technology news on Phys.org
  • #2
What kind of class is this for just curious?
 
  • #3
Futurestar33 said:
What kind of class is this for just curious?
Thanks for the quick reply.

This is for a capstone project for a senior design.
 
  • #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.
 
  • #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
 
  • #6
kroni said:
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.
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.
kroni said:
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
I agree completely. I never use pow() to compute the square of a number
kroni said:
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.
 
  • #7
kroni said:
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.
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.
kroni said:
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
I agree completely. I never use pow() to compute the square of a number
kroni said:
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.
 
  • #8
justin cooper said:
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;
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.
 

Related to Re-entry Projectile Analysis Program

1. What is a Re-entry Projectile Analysis Program?

A Re-entry Projectile Analysis Program is a computer program designed to simulate the flight and impact of a projectile that re-enters the Earth's atmosphere. It takes into account various factors such as the shape and composition of the projectile, the atmospheric conditions, and the trajectory to accurately predict its behavior.

2. How does a Re-entry Projectile Analysis Program work?

The program uses mathematical models and equations to calculate the forces and factors that affect the trajectory of a projectile during re-entry. It takes into account the mass, velocity, and shape of the object, as well as the atmospheric conditions such as air density, temperature, and wind speed. Using these inputs, it can accurately predict the flight path and impact location of the projectile.

3. What types of projectiles can be analyzed using this program?

A Re-entry Projectile Analysis Program can be used to analyze a wide range of projectiles, including spacecraft, missiles, and other man-made objects that re-enter the Earth's atmosphere. It can also be used to study natural phenomena such as meteorites and asteroids.

4. What are the benefits of using a Re-entry Projectile Analysis Program?

Using this program can help scientists and engineers understand the behavior of projectiles during re-entry, which is crucial for designing safe and efficient space missions. It can also aid in predicting the potential impact and damage of a projectile, allowing for better risk assessment and mitigation strategies.

5. Is a Re-entry Projectile Analysis Program accurate?

Yes, a Re-entry Projectile Analysis Program is highly accurate when used with correct input data. However, due to the complex nature of atmospheric conditions and varying factors, there may be some margin of error in the predictions. It is important to continuously validate and refine the program to improve its accuracy.

Hot Threads

Recent Insights

Back
Top