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-26 kg – average mass of the molecule of the air,

=3.65×10-10 m – average collision diameter of the molecule of the air

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 92 km, below 72 km a usual Sutherland formula is used. Between 72 and 92 km the density is assumed to be constant and equal to 1.425×10-5 kg/s/m, which is close to the average value for 82 km, given by different standard atmosphere. The reason for such sophisticated procedure is that numerous sources are pointing that Sutherland formula is valid only below some certain altitude, though different sources claim different altitudes (the lowest is above 72 km and the highest is below 92 km). In the present version the value of molecular viscosity is used only for Reynolds number calculation, which is used only for reference at such high altitudes. But in the future probably it will be necessary to use viscosity to calculate some other parameters, for example transition criteria for switch between different heat transfer formulas. In this case it may be useful to use some variation of viscosity with altitude in the intermediate range.

 

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 Stanton number,  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:

  1. 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:

  1. For temperature below 1400K a 4-th degree polynomial is used

  1. There are only 4 values of reference pressure: 0.001, 0.01, 0.1 and 1.0 bar.

 

 

Stanton number in free-molecular flow is calculated as

,

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 Stanton number is given by the same formula, but recovery coefficient is no more constant, but calculated as

,

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 Stanton number a simplified Fay-Riddle formula is used

,

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 Stanton number as product with density, and pressure and density – as ratio.

 

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 Stanton number is introduced, which was originally proposed by Chang basing on the analysis of the experimental data. The correction factor is calculated as

 

,

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 Stanton number is calculated using so called “friction analogy”:

.

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 Stanton number is also different:

 

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 130 km it is extremely small comparing to convective heat flux and can be neglected. Radiation from the shock wave is also rather small because of low density.

 

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 85 km and 76 km, while REST gives 1.30 and 1.02, which is quiet reasonable accuracy.

 

The detailed datasheet with results of Fotino-M2 re-entry simulation, starting from 270 km (separation from Foton) is presented in FotinoM2_270k_latest.xls file. Initial parameters for the simulations are stores in FotinoM2_270km.rest file.

 

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 Newton theory. With application of empirical correction factors modified Newton theory, which is usually valid only for , 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 Stanton number calculation it would be better to use local Knudsen number, based on local meen free path near the wall, or some other local criteria. Unfortunately it is hard to derive such criteria for re-entry case. It is much easier to trace parameter changes when Knudsen number increases, than when its decrease. For example it is much easier to check whether shock relations are STILL valid, than to check if they are ALREADY valid, because transitional formulas doesn’t give evidently wrong results when Knudsen number decreases (they just predict higher heat flux, resulting in higher wall temperature).

 

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.