- #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
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