The order of calculations and the formulas for meen free path and molecular viscosity were changed. Now first the meen free path is calculated as
Note: In reality due to the changes in chemical composition of atmosphere with altitude these two parameters are not constant
- density, kg/m3, obtained from standard atmosphere subroutine
After this molecular viscosity is calculated as
where - thermal speed of air molecules or “meen molecular speed”,
- Boltzman constant,
T – temperature, obtained from standard atmosphere subroutine
The above formula for molecular
viscosity is used for altitudes higher then
Speed of sound, used for Mach number calculations, is now calculated with more accurate account for the influence of temperature. In a well-known formula
(specific heat ratio) and (gas constant for air) are no longer constants but functions of temperature. Gas constant is expressed as
, where and are specific heats at constant pressure and constant volume. In turn
To calculate and special functions are used, which are covered in section describing the new heat flux model.
A simple aerodynamic model, based on experimental results for the sphere, was implemented. For free-molecular flow the value of drag coefficient is constant and equal to 2.20.
Note: The real value of drag coefficient depends upon the wall temperature and wall roughness. For the sphere it can vary from 2.00 to 2.60
Flow is considered to be free-molecular when Knudsen number, calculated as
where is reference length (i.e. diameter of the capsule),
is above the certain limit. Though this limit was left as user input, most sources recommend for the sphere the value of 10.
For transitional flow the drag coefficient varies with Knudsen number. The following formula was use to approximate this variation
where 1.55 is the average between the drag coefficient of the sphere in free-molecular and in continuum flow, 0.65 is the difference between this value and drag coefficient in free-molecular flow, is the abovementioned limit and is another limit, below which the flow can be considered as continuum. This limit was also left as user input, for the sphere most sources recommend the value of 0.001.
For continuum flow, i.e. when
the drag coefficient depends upon Mach number. For it is constant and equal to 0.90. For the drag coefficient increases linearly with the decrease of Mach number from 0.9 to its maximum value 1.00, thus given by the formula
Note: Some sources give higher values for maximum drag coefficient of the sphere, up to 1.1 or 1.2
For drag coefficient remains constant and equal to its maximum value of 1.00. For drag coefficient drops down to its subsonic subcritical value of 0.45. This variation is given by the formula
For an annular supersonic region exists at the surface of the sphere, which is closed by the shock located near the equator. This shock causes the separation of the boundary layer, no matter what’s the Reynolds number is. But for lower Mach numbers the shock vanishes and the type of boundary layer governs the location of the separation line, i.e. by the value of Reynolds number. If at the moment when the Reynolds number is lower than some critical value, which was left as user input, then the flow is considered to be subcritical (i.e. laminar boundary layer with early separation) and the drag coefficient remains constant and equal to 0.45, since .
Note: There is a wide spread of the values for subcritical drag coefficient of the sphere. 0.45 is a kind of average one.
When Reynolds number reaches the critical value, a transition occurs which corresponds to the drop of drag coefficient to supercritical value of 0.20.
Note: Once again, there is no certainty for this value, though the spread is not as big as for subcritical one.
Usually this drop is considered to be sharp and is highly influenced by flow turbulence and surface roughness. However, the available experimental data for free-falling spheres show that for them the situation is different, and instead of a sharp drop one can expect a rather smooth decline. It usually starts at , which is the default and recommended value of in the present version, and ends with . The variation of drag coefficient for this range is given by the formula
If it happens so that at the moment when the Reynolds number is already higher than the user supplied value of , then the decline of the drag coefficient starts right from that moment, with .
For drag coefficient remains constant and equal to 0.20.
A simple attitude model was added, which considers oscillations in pitch only (i.e. in plane) of the trajectory. The change of the angle of attack is governed by an ordinary differential equation
where is aerodynamic moment, [N×m], for the sphere with the center of mass displaced from the geometrical center by it is equal to , where is drag,
is the moment of inertia, [kg×m2], referenced to the center of mass; considering that the mass of the sphere is uniformly distributed it’s equal to ,
is the slope of moment curve; taking into the account the expression for the moment above and the standard definition of moment coefficient, ,
is dynamic pressure, is projected surface of the sphere.
This equation is integrated together with equations of motion of the capsule. The coupling of attitude simulation with the other parts of REST is one-way, i.e. values of are taken from REST for the calculations, but the change of angle of attack doesn’t affect the motion in any way.
Note: As the angle of attack and its time derivative are fast varying oscillating variables, the automatic adjustment of time step of the integrator can slow down the whole simulation.
To integrate the above equation initial values of angle of attack and angular rate, together with the displacement of the center of mass for moment of inertia calculations, are required. User in the auxiliary form, which is called by “Simple attitude” button, can provide these values. The results are written into separate file, “RESTangle####.txt” – time, angle [deg] and angular rate [deg/s].
Note: In the present version the same output frequency is used for the results of the simple attitude simulation and for other results. Choosing this frequency too low (lower than double frequency of the oscillations) can give confusing and misleading output.
In the present version local convective heat flux in the stagnation point is calculated with a universal formula, which is valid for any kind of flow:
where is recovery factor, is atmospheric temperature, is wall temperature, is stagnation pressure.
The implemented formula allows calculating both heating and cooling, contrary to some simple engineering formulas.
The approach for calculation of recovery enthalpy is conservative (in case of proper choice of recovery factor, which account for entropy rise in continuum compressible flow and for energy transfer between different degrees of freedom in free-molecular flow). It assumes that the wall is fully catalytic, i.e. even if there is dissociation in the shock wave (in supersonic continuum flow), which drains part of the energy from the flow, all the energy spent on dissociation is transferred to the wall because a recombination occurs on it. Fully catalytic approach means that all the species produced in dissociation recombinate back to their original state, disregarding the wall temperature.
For the purpose of heat flux calculations specific heat is not considered constant but varies with temperature and pressure, due to excitations of additional degrees of freedom and chemical reactions in the air. Polynomial approximations were used to fit the available experimental data for real air. The procedure to calculate specific heat is the following:
2. The range from 1600K to 5800K is divided into 3 subranges: from 1600K to 3000K, from 3000K to 4400K and from 4400K to 5800K. Polynomials of the 9-th degree are used for each subrange. First two reference values of pressure are chosen, above and below the provided pressure: . The reference values are 0.001, 0.002, 0.004, 0.01, 0.02, 0.04, 0.1, 0.2, 0.4, 0.6 and 1.0 bar. For both reference values specific heat is calculated as
where are centering and normalizing constants.
3. Calculated values of specific heat are interpolated to the provided pressure
where is converted first from Pa to bars.
The same kind of procedure is implemented for calculation of specific heat ratio, with few exceptions:
where is energy accommodation factor, a user input (recommended value for engineering practice is 0.9, in reality it depends upon the material and roughness of the surface).
Recovery factor for free molecular flow is constant and equal to 1.1667 (or ). For the purpose of heat flux calculations the flow is considered to be free-molecular is Knudsen number based on nose radius, , is higher than , where factor of 3 was adjusted by trial and error to provide smooth translation form free-molecular to transitional flow. The disadvantages of such approach and the possible ways to improve it are discussed in the attachment.
For transitional flow
where is the value of recovery factor for continuum flow.
Continuum flow starts when . Factor of 500 was also adjusted by trial and error. It accounts to the fact that though from the general point of view the flow past the capsule exhibits free-molecular nature till quiet low Knudsen numbers, in the vicinity of the stagnation point it becomes dense enough for continuum effects to play a major role even at pretty high Knudsen numbers. Though factor of 500 is an overestimation, it is corrected during the calculations with special empirical adjustment factor, accounting for the rarefication.
For continuum supersonic (i.e. ) flow the recovery factor is constant and equal to 0.842 (or , where is Prandtl number, for general engineering estimations it’s usually considered to be equal to 0.71 for laminar flow).
where subscript ‘st’ marks stagnation parameters behind the shockwave, and
for the sphere.
Stagnation parameters behind the shockwave are calculated using the following relations
is calculated basing upon with Sutherland formula if , of with the formula of the same structure but with higher power (0.74 instead of 0.5); the choice of the coefficients in the second formula provides a smooth transition from Sutherland law.
Subscript ‘2’ in the above formula refers to the parameters behind the shock-wave. These are calculated using Rankine-Hugoniot relations
For the sake of simplicity it was considered that . This of course is not true for high-speed rarefied flows, also Rankine-Hugoniot relations are not valid for this case. In reality the temperature rise in shock wave results in a dissociation, which has two consequences – part of energy is drained fro the flow and its thermophysical properties are changed due to the change in chemical composition. These issues are discussed in more details in the appendix.
Utilization of Rankine-Hugoniot
relations results in the overestimation of temperature behind the shock wave
and underestimation of density and pressure. But this is compensated by the
fact that viscosity (which is proportional to temperature and so also
overestimated) appears in the formula for
The original Fay-Riddle formula was
derived for continuum flows, where is a distinction exists between the shock
and the boundary layer. In REST it is implemented for the flows where the
shock-wave and boundary layer are merged together due to rarefication (so
called “merged layer” or “shock layer”). Original formula also accounts for
dissociation effects by multiplier containing Lewis number, which is very
difficult to estimate within simple engineering calculations. To count the two
above-mentioned effects (rarefication and dissociation) a special correction to
Note: In case of modifying the model (e.g. using another transition criteria or switching thresholds) it’s necessary to check the value of - the Fay-Riddle formula with this correction can be used only if , below this value a formula for transitional flow should be used.
For subsonic laminar (subcritical)
For compressible flow (i.e. ) the effective Reynolds number is based at the stagnation parameters:
the former are calculated as
(used for calculation of with Sutherland formula).
For incompressible flow .
Turbulent flow is characterized by
higher value of Prandtl number () and thus higher
value of recovery factor (). The formula for
Radiative heat flux also plays an important role in the estimation of stagnation point temperature. It can be divided into two parts – radiation to the capsule and from the capsule. Unfortunately the radiation heat exchange is a very complex phenomena by itself. One if the simplest model is used for its estimation in REST – a black body radiation model.
The first part of radiative heat flux is calculated in two different ways: 1) for free-molecular, transitional and subsonic continuum flow it is given by standard black body formula
where is the absorbtivity of the surface, is Stephan-Boltzmann constant and for free-molecular and transitional flow (which is quiet a conservative assumption) and for subsonic continuum flow;
2) for supersonic continuum flow the radiation of the hot dissociating gas in shock wave is calculated with the approximate engineering formula
Note: Incoming radiative heat flux plays
significant role only in free-molecular flow. For altitudes lower than
The second part of radiative heat flux is given also by black body formula
where is emissivity of the surface, for black body . 0.9 is quiet a usual engineering value for it and 0.8 is pretty conservative. Though some sources warn that due to heating and burning this value can degrade down to 0.4, which has a sufficient impact on radiative cooling, and radiative cooling is a key factor for the survival of the re-entry capsule. Without radiative cooling the surface temperature can reach 6000K!
The cumulative heat flux
is then used to calculate the wall temperature by integration of ordinary differential equation
where is the surface density of the skin [kg/m2] and is specific heat of the material, both values are user input.
The initial condition for integration is .
Coupling with Andrea’s TPS tool (v. 3.0beta)
In order to increase accuracy for wall temperature predictions and extend the capabilities of REST a full coupling was made to Andrea’s TPS tool. For this purpose some modifications were made both to REST and to TPS tool. A 3-layer version of the tool was chosen for coupling. Calculation of wall heat flux with simple engineering formulas was replaced with input from REST. Boundary conditions for the last node (inside the capsule) were changed to zero temperature gradient. From REST side all heat flux calculations were doubled (for comparison purposes) – separate calculations are made using wall temperature predicted with simple model and wall temperature predicted by TPS tool.
To provide user input for TPS simulation a “TPS” button was added to the main REST form, which invokes the TPS parameters form. Because this code is on beta stage, only thickness of the layer and number of nodes for it can be changed. Changing material type doesn’t have any effect. For release version a material database from TPS v.10 can be used.
Note: Andrea’s tool performs internal integration with the time step provided by REST code. There is no back influence at integration time step, i.e. it is not adjusted by REST integrator to temperature variation. Because the integration procedure implemented in TPS tool is rather simple, providing too high values of initial time step can result in loss of accuracy and stability and abnormal program termination. IT IS STRONGLY RECOMMENDED TO KEEP INITIAL TIME STEP BELOW 0.1 S!
Though no extensive validation for the new versions of REST were done except numerous test runs, it should be noticed that the results obtained are at least look physical, i.e. there are no unexplainable peaks or pitfalls in parameter profiles, no oscillations for the parameters which are expected to vary monotonously and etc. Some comparison was made with the results of Marco’s simulations using Fastran CFDRC tool (a finite volume CFD code for continuum flow with the account for chemical reactions). Marco’s simulations for 0.4-m sphere with the conditions corresponding to the expected point of maximum heat flux (altitude 76...77km, =5520m/s, =1700K) predict convective heat flux at stagnation point about 0.6 MW/m2, while REST simulations predict 0.58 MW/m2 (the lower value in REST can be explained by the fact that for this moment predicted wall temperature is higher and makes 1900K). Simulations with DSMC code DSWT by Bird gives for the same conditions but a bit higher wall temperature (=1800K) the value of convective heat flux at stagnation point about 0.59 MW/m2.
Concerning predictions of drag
coefficient DSWT code gives the values of 1.25 and 1.03 for conditions
corresponding to altitudes of
The detailed datasheet with results of
Fotino-M2 re-entry simulation, starting from
Appendix: Drawbacks of the present models and possible ways of improvement
The present aerodynamic and simple
attitude models allow calculations for SPHERICAL CAPSULES ONLY! Aerodynamics
forces and moments for other shapes can be calculated by 3D panel model for
free-molecular, transitional and hypersonic continuum flow. For free-molecular
and transitional flow rarefied gas theory can be use, and for hypersonic
continuum – modified
To switch between different formulas for
To calculate parameters behind the shock wave in a more accurate way one can start from conservation equations with the account of dissociation:
together with state equation written in the form
Here specific heats and specific heat ratio are functions of temperature and pressure (usually some kind of approximations for real gas experimental data), is mass concentration of atomic oxygen behind the shock wave, is mass concentration of atomic nitrogen behind the shock wave, both also functions of temperature and pressure (they also can be found as approximations of experimental data), and are specific dissociation enthalpies for oxygen and nitrogen (constants).
This system of equations is non-linear, but some sources point that it’s possible to build an iterative procedure for its solution. For known parameters behind the shock wave stagnation enthalpy can be found as
Recovery enthalpy than can be calculated as
where are recombination factors for oxygen and nitrogen, depending upon wall temperature and wall catalysity.
To increase the accuracy of TPS simulation it is possible to turn to axysymmetrical problem, using Fourie equation instead of simple model implemented in Andrea’s tool:
where is thermal conductivity (a function of coordinates in case of multi-layer TPS and also temperature for some materials).
Boundary conditions for this equation are of so-called Neumann type, with temperature gradient at the wall (outside boundary) proportional to local heat flux and zero temperature gradient inside the capsule. For the sphere local heat flux distribution is given by the formula
, where is polar angle ( at the nose).
For other simple shapes different engineering estimations are also available.