Hi all:(adsbygoogle = window.adsbygoogle || []).push({});

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

**Physics Forums | Science Articles, Homework Help, Discussion**

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Re-entry Projectile Analysis Program

**Physics Forums | Science Articles, Homework Help, Discussion**