Changes implemented in REST
v. 2.02sf2 -> v. 2.02sfa2 -> v. 3.0beta
Post-graduate student Anton Lyaskin
General atmospheric calculations
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
,
where 
=4.81×10-
=3.65×10-
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.
Aerodynamics calculation
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.
Simple attitude simulation
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.
New heat flux model
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 
are recovery and wall
enthalpies, given by formulas
,
,
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:
- For
temperature below 1600K a linear variation with temperature is assumed
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:
- For
temperature below 1400K a 4-th degree polynomial is used
- There
are only 4 values of reference pressure: 0.001, 0.01, 0.1 and 1.0 bar.
,
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).
To calculate
,
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
,
where 
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)
flow
.
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!
Validation and testing, results of Fotino-M2 simulations
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 
, can be extended nearly to the whole supersonic region.
Unfortunately there is no simple but universal way to predict aerodynamic
performance of complex-shape capsule in transonic and subsonic flows. For this
purpose CFD modeling is required.
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:
(mass conservation)
(momentum
conservation)
(energy conservation)
where 
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.