Physical and computational aspects of chemically reacting hypersonic flows

Computer Methods in Applied Mechanics and Engineering North-Holland

111 (1994)

1-35

Physical and computational aspects of chemically reacting hypersonic flows Ioannis

John Argyris St. Doltsinis, Heinz F’riz, Jiirgen Urban

Institute for Computer Applications, Stuttgart,

Germany

Received 8 December 1992

In continuation of the research and development work at the ICA in the field of hypersonic reentry flow simulation, this paper proposes techniques for a reduction of the computational effort which is associated with chemical reactions as well as viscous effects in the flow field. First, the numerical effort for the computation of steady-state solutions of inviscid and viscous flows may be reduced by application of a local time-stepping technique. A further reduction is obtained by a rigorous adaption of the integration scheme for the chemical kinetics to the characteristics of the flow problem. The number of numerical operations may be reduced on the one hand by adaptive techniques which restrict the viscous as well as the reactive region of computation and, on the other hand, by a simplified description of the physical transport coefficients. Appropriately selected numerical examples demonstrate the efficiency of the proposed techniques.

1. Introduction Within the Research and Development_ Program (R&D) and the Research for Qualification Program (RfQ) for the Aerodynamics of the European Space Project “Hermes”, the Institute for Computer Applications (ICA) is entrusted with the investigation of problems related to the numerical simulation of chemically reacting hypersonic viscous flow. For this purpose, specific computational tools have been developed at the ICA in the context of a general CFD methodology. A first stage of the research work was documented in [l], and comprises the solution of the compressible Euler and Navier-Stokes equations as well as the consideration of chemical reactions in mixtures of inviscid gases. An investigation of the numerical scheme for the time integration of the chemical kinetics and the extension of the numerical method to the viscous reactive case was presented in [2]. For a large-scale application of the computational techniques developed so far for the numerical simulation of hypersonic flow, we considered at the ICA [3] the symmetric re-entry maneuver of a complete three-dimensional configuration of the Hermes space-shuttle as shown in Fig. 1.1. The maneuver conditions are specified by the Mach number M, = 25, the angle of incidence of the longitudiual axis of the shuttle (Y = 30” and the atmospheric conditions at the altitude H = 75 km. An impression of the finite element discretisation of the surface of the object, the plane of symmetry and the outflow plane is given in Fig. 1.2. They comprise 9,907 finite surface elements defining a three-dimensional space covered by 110,488 hexahedral finite elements. The discretisation mesh has a total of 117,973 nodal points corresponding to 589,865 unknown field variables in the discretised Euler equations, excluding chemical reactions. The computation up to steady-state conditions in the flow field applies the technique of the local timestepping and requires 2,600 increments of time. The CPU-time on the Cray 2 operated at the University of Stuttgart was 8 hours. As a result of the computation, Fig. 1.3 shows the distribution of the Mach number in the flow field around the shuttle. Figure 1.4 demonstrates the variation of the temperature in 0045-7825/94/$07.00

@ 1994 Elsevier Science B.V. All rights reserved

J. Argyris

2

et al., Physical

and computational

Fig. 1.1. Hermes space shuttle.

Fig. 1.3. Mach number distribution

Fig. 1.5. Temperature

distribution

(3-D simulation).

(2-D simultation).

aspects

of reacting

Aows

Fig. 1.2. Layout of finite element mesh.

Fig. 1.4. Temperature

Fig. 1.6. Temperature (2-D simultation).

field (3-D simulation).

distribution

for reacting

flow

J. Argyris et al., Physical and computational

aspects of reacting flows

3

the field indicating local values in the range of 25,000 K. Such temperature values pose of course serious problems to the structure of the shuttle. They are, however, not realistic due to the assumption of air as an ideal gas. Computations taking into account real gas effects as by chemical reactions between different components of the air considered as a gas mixture were performed for a two-dimensional variant of the problem modelling the plane of symmetry of the shuttle. Here, the computation of the Euler flow again required 2,500 incremental steps up to the steady-state solution, whilst the consideration of chemical kinetics increased that number by a factor of 2 to 5,000 steps. The temperature distribution for the ideal gas is depicted in Fig. 1.5 and maximum values range about 26,000 K. Figure 1.6 shows the temperature in the field for the chemically reacting two-dimensional flow and documents that maximum temperatures now are in the range of merely 12,000 K. The above suggests that real gas effects as chemistry and also viscosity have to be accounted for in the numerical simulation and they reduce the requirements for the realisation of the shuttle. On the other hand, they increase the necessary computer capacity enormously and demand a refined simplification of the physical modelling and the numerical computation. The incorporation of chemical reactions and viscous effects into the flow simulation leads to a considerable increase in the computational effort. As a matter of fact, the consideration of reaction kinetics in the three-dimensional case considered previously would double the number of unknowns to 1,179,730 and also requires about twice the time steps of the non-reactive Euler computation. In addition viscous effects demand an appropriate refinement of the discretisation mesh. In the present account, techniques for the reduction of the computational effort are discussed. For this purpose, the application of local time-stepping is investigated for the chemically reacting inviscid flow. Numerical studies on the flow around a double ellipse demonstrate the improvement of the rate of convergence of the numerical solution. The treatment of chemical reactions in the flow simulation requires additional computational time for a number of reasons. First, the computation of the production rates for the species and the incremental integration requires considerable computational effort. In addition, the occurrence of chemical reactions may increase the stiffness of the flow problem and complicate the numerical procedure. According to the characteristics of chemically reacting hypersonic flows, two approaches are considered in order to reduce the computational effort. The first approach concerns an adaptive application of different integration methods for the chemical kinetics to the instant state of the reacting flow field. The second approach is based on the fact that areas where chemical reactions markedly affect the flow field are restricted to high-temperature regions following shock waves. To achieve a reduction of the number of numerical operations, the integration of the chemical kinetics may be restricted to a certain domain with considerable chemical production rates. Studies on the twodimensional problem of an inviscid chemically reacting flow over a wedge demonstrate the efficiency of the considered techniques. In principle, local time-stepping may be also applied to viscous flow problems. In the context of the Taylor-Galerkin method, however, this affects the convergence behaviour as well as the steadystate solution of the problem. In general, the incorporation of viscous phenomena requires additional computational effort in the numerical treatment of the flow problem. Additional operations have to be carried out in order to obtain the diffusive contributions in the finite element equations. Since areas where diffusive transport markedly affects the flow field are restricted to near-wall regions, a specific adaptive activation of the respective modules will be considered aiming at a reduction of the number of numerical operations. The consequences of the application of local time-stepping and the adaptive activation of physical operators on the flow field is demonstrated on the chemically reacting viscous flow over a wedge. An important difference between inviscid and viscous flow concerns the specification of particular boundary conditions for the latter. In this connection, a brief survey of appropriate boundary conditions at solid walls is presented. The employment of the physical models for the transport coefficients for the numerical treatment of a flow problem increases the computational requirements because of the complex structure of the expressions. In order to reduce the number of numerical operations, it appears advantageous to replace the complex formulations by simplified expressions adopting similarity considerations. The consequences of a reduced complexity in modelling and the effect of the different solid wall models on a viscous boundary layer are studied for the case of the flow over a flat plate.

4

J. Argyris

et al., Physical and computational

aspects of reacting Aows

2. Inviscid hypersonic flow with chemical reactions The basic equations governing the hypersonic flow within the compressible regime are the well-known Euler equations. They may be expressed in the conservative form

au

wl af2 af3 at ax1 ax2 ax3

(1)

-+-+-+----cl.

In this formulation, t denotes the time and x1,22, x3 the spatial Cartesian coordinates. The vector u comprises the five conservation quantities Q, evl, QQ, ~~13,Q& as the independent field variables; i.e.: u = {Q, QUA,p2, ~213,QE}. These incorporate the gas density e, the components of the velocity vector v and the total specific energy & of the gas. The flux vectors fi, f2, f3 contain the dependent field variables as detailed in [11. The occurrence of chemical reactions in a multicomponent gas mixture necessitates additional considerations. In order to evaluate the convective fluxes in the Euler equations, the caloric and thermal properties of the gas mixture have to be determined. As indicated in [l], the determination of the above properties requires a knowledge of the local gas composition which has to be deduced from the solution of the equations for the conservation of species. In collective form,

aw

ah4

+ a(V24 + ab34

at+-ax1 - ax2

3g-=w

.

where w = (~1 ~2.. . pr . . . pN} and cj = {ijl L&J.. .& . . . ijN} are N x 1 vector columns referring to the N species. The mass conservation is expressed in terms of the local time rate of the specific masses w of the individual substances, the convective mass transfer and the production rates per unit volume ti. Due to the fact that N, chemical elements have to be conserved, only N - N, reaction rates tir are linearly independent. Thus, 5 + N - N, differential conservation equations have to be solved in conjunction with N, algebraic equations for the conservation of elements. As the species equations are strongly coupled with the flow equations (1) via the velocity field, they have to be solved simultaneously with the latter. The combined system of differential equations may be presented in a single matrix form

:[:I+&[$w] +&Jufw] +&[vfw] =[l]. In order to integrate eqs. (3), it appears advantageous to distinguish between flow and chemical phenomena because they occur in general on significantly different time scales. In this connection, the chemical reactions are first, assumed frozen (Lj = o) in order to obtain an initial estimate for the flow field in the course of an incremental integration procedure. The change in the gas composition due to chemical reactions is performed separately by integration of the rate equations ;( = F(T, w). The ultimate incremental solution for the reacting gas mixture is then determined by

(4) in which the index f refers to a frozen solution. Some appropriate numerical methods for the integration of the rate equations as well as an investigation of their numerical stability have been presented in [2]. An improved version of these methods is introduced in section 3.3.

3. Improved solution methods for inviscid chemically reacting flows 3.1. Iterative

solution

procedure

The numerical integration of (3) is performed in an explicit two-step finite element procedure. Here, the governing equation (3) is first discretised with respect to time by a Taylor-series expansion, and is subsequently approximated in the space domain via the Galerkin technique [l].

J. Argyris et al., Physical and computational

Application of the finite element Galerkin technique equation governing the entire finite element domain MUAf

= ?-[R -

5

aspects of reacting Bows

to the spatial discretisation

leads to a global

s] .

(5)

Here, M denotes the consistent “msss”matrix of the problem. The vector UAf = A Uf = bUf - “U comprises incremental changes of the independent field variables at the nodal points of the finite element mesh. These variations arise from the numerical integration of the flow equations within the small interval in time r = bt - at, in which the chemical reactions are first assumed to be frozen. Correspondingly, the “force”vector R and the “stress”vector S comprise the resultants of the individual element contributions at the mesh nodal points. A solution of (5) for the unknown vector bU may be obtained via the iterative scheme uAf,j+1 = uAf,j + Mil

{r [R - S] - MA Uf,j} .

In the recurrence formula (6), we apply the diagonal lumped “mass”matrix MJ_ of the system as an iteration matrix. This result corresponds to a multicomponent frozen gas mixture buf = “u+

The ultimate incremental

bU=

buf

(7)

UAf . solution for the reacting gas mixture is obtained from

+ VA,.

(8)

The determination of the species contributions and is discussed in a subsequent section. 3.2. Local time-stepping

for steady-state

due to chemical reactions

VA, is performed separately

Aow problems

The solution of a flow problem is called the steady-state remain locally unchanged with respect to time, thus

solution if the independent

field variables

(9)

UA = 0.

As a consequence following from (6) and (8)) the variation of the independent field variables due to convective transport must be compensated by appropriate production rates resulting from chemical reactions

ML1 {T [R -

s]} + VA,

= 0.

(10)

If steady-state problems are investigated, a local time-stepping procedure may improve the rate of convergence of the numerical solution. The associated finite element equation then reads MUAf

= T[R - s] .

The diagonal matrix T = I71 . . . T,

(11)

T comprises the local stability limits . . .

T,,,,] .

The local time-step is obtained as the minimum from the Courant-F’riedrich-Lewy on the finite elements connected to an individual mesh point

(1‘4 condition applied

(13)

3. Argyris et al., Physical and computational

6

aspects of reacting Aows

In (13) I,, represents a characteristic length of the finite element and o, a the element average velocity and speed of sound respectively. Again, a solution for the unknown vector ‘U is obtained via an iterative scheme similar to (6). The incremental solution for the reacting gas mixture is given by

bU= “U+ tiAf + ti& where @of and u& 3.3. Integration

(14)

refer to the flow and chemistry step in analogy to (8) respectively.

of the chemical rate equations

The system of ordinary differential equations describing the production and consumption of the participating substances due to chemical reactions may be written in the vector form ti = F(w).

(15)

The functions F(w) describe the non-linear correlation of the production rates to the concentration of the considered species [Z]. The computational procedure for the flow field requires in (4) an integration of the non-linear functions in a small interval of time 7 = ‘t - 4. An appropriate integration scheme reads

s 2t

‘w=

‘w

‘t

Fdt=

‘w+(~-C)?~F+<~;~F

where the parameter is taken as C = 0 integration procedure. A brief account investigation of the numerical stability intervals 7 for the chemical reaction and

(16)

in an explicit treatment and 0 < C 6 1 in any kind of implicit on some suitable forms for the evaluation of eq. (16) and an are outlined in [a]. Because of stability requirements, the time r for the flow evolution may differ. The ratio of the time intervals

(17) represents a subincrementation for chemical kinetics if the integer (Ic) > 1. In practical flow problems, the stability limit for the integration of the flow equations may represent a limit to the global time increment as well. If this situation arises, the integration of the chemical kinetics may be performed with constant production rates for a certain number of time-steps. If steady-state problems are investigated by local time-stepping for the flow, the corresponding scheme for the chemistry reads 2w = lw + (1 - c)k-’

‘TIF +
(18)

It has to be noted that the local stability limits for the integration of the flow and the chemistry may vary within the flow field and thus k also. As a consequence, the integration of the chemistry may be treated at each nodal point according to the individual stability requirements. In order to facilitate vectorisation, the treatment of the chemistry is adapted to the stability requirements of the critical nodal point. Following the account in [2], two suitable forms for the evaluation of (18) are considered a) Explicit

Euler-Cauchy

integration

In the special case of an explicit treatment

2w = lw + k-l lTIF

with <=O, eq. (18) is reduced to

(1%

J. Argyris et al., Physical and computational

and reproduces the Euler-Cauchy less than

integration

aspects of reacting Bows

7

scheme. For stability, the step-size F in (17) must here be

Here, p is the spectral radius of the Jacobian matrix and may be estimated by inspection diagonal. For a non-oscillating solution, this upper limit has to be reduced by a factor of two.

of the

b) Method of &m-wised direct solution A second technique to solve the implicit scheme (18), which has been applied with some successs in [4], is based on the linearisation of the implicit term. Expressing 2F by a Taylor series and neglecting terms of order higher than one, we obtain the linearised scheme

2F=

IF-

%V[2w-“w]

where N is the Jacobian find

matrix and is given by N = -(W/aut).

(21) Substituting

eq. (21) in eq. (18), we

with l& = [I+

k-l ‘!l’C 'N]

-l .

(22b)

As outlined in 121, the determination of the stability limit for this method implies a considerable effort due to the necessary differentiation of the matrix Q. c) Chemically reacting

hypersonic flow past a double ellipse

In the following, the proposed local time-stepping algorithm is studied for the two-dimensional case of the flow around a double ellipse. Particular attention is paid to the differences in the numerical results as well as in the reduced computational effort as compared to the conventional algorithm. For the numerical treatment of the problem, the flow domain around the double ellipse is idealised by a regular mesh consisting of 4898 triangular elements. The discretisation, plotted in Fig. 3.1, has been used for the computation of the inviscid frozen and the inviscid reacting flow. The flow medium is taken to be atmospheric air with the Mach number A& = 25 at an altitude of 75 km. The associated farfield conditions are obtained from the IS0 standards. The numerical simulation of the appertaining flows is performed for an angle of attack (x = 30°. The chemical reactions considered here are dissociation-recombination and exchange reactions of certain substances in atmospheric air. Neglecting ionisation, we may apply a reaction model for the components 0, N, NO, 02 and N2 which consists of 17 reaction equations [5]. As demonstrated in Fig. 3.2, the application of local time-stepping in fact results in improved convergence. A quantitative assessment of the efficiency of the local time-stepping scheme may be based on a comparison of the CPU-time which is required to reduce the density residuum by two orders of magnitude. Time measurements reveal that the application of local time-stepping reduces the computational effort by 65 % in the non-reacting case and by 60 % in the reacting case. The temperature distribution for the non-reacting and reacting case is plotted along the stagnation stream line (Fig. 3.3.a) and along the body surface (Fig. 3.3.b). The figures illustrate clearly that the steady-state solution is not significantly affected by the choice of the numerical procedure. These results were obtained with the explicit integration method for the chemical kinetics. The following studies may demonstrate that an improved adaption of the solution strategy to the characteristics of the numerical problem will further reduce the computational effort for the computation of inviscid flows with chemical reactions.

J. Argyris et al., Physical and computational

4898 triangular

aspects of reacting flows

elements

2575 nodes

Fig. 3.1. Inviscid hypersonic

3.4. Eficient

solution

strategies

flow past a double ellipse. Layout of finite element mesh.

for chemical

kinetics

The treatment of chemical reactions in the flow simulation requires additional computational time due to various effects. First, the computation of the source terms and the integration of the rate equations requires considerable computational effort. In addition, the occurence of chemical reactions increases the stiffness of the numerical procedure. Depending on the flow problem, the wave propagation may be blocked by the very small time-steps which are necessary for the integration of the chemical kinetics. According to the characteristics of chemically reacting hypersonic flows, two approaches are considered in this section in order to reduce the computational effort. The characteristics of reacting flows may be demonstrated by the study of the inviscid chemically reacting flow over a 30” wedge (Fig. 3.4.a). The flow medium is taken to be atmospheric air with the Mach number M, = 22 at an altitude of 65 km. The flow field over the wedge has been discretised by a mesh of 1512 triangular elements, as plotted in Fig. 3.4.b. The evolution of the flow field in chemical non-equilibrium is studied within a distance of 1 = 1OOOm from the wedge apex. While the how tends to be in a frozen state near the apex, it tends to achieve chemical equilibrium far downstream. The variation of certain thermodynamic state variables between those limiting states is studied at certain positions E = 0.05 m, 0.5 m, 5 m, 50 m, 500 m. Profiles of the considered variables are taken normal to the wedge surface. In order to obtain a sufficient resolution of the flow field at the considered positions, it seemed to be advantageous to scale the numerical problem. For this purpose, a number of individual problems was analysed numerically with the side length of the wedge varying from lo = 0.1 m to Ec = 1OOOm. The numerical solutions, although obtained for different wedge lengths, correspond to the same physical problem. Non-equilibrium chemistry is included by explicit integration of the chemical kinetics. The numerical solutions for non-equilibrium chemistry are compared

J. Argyris et al., PhysicaJ and computational

aspects of reacting Bows

9

10-I

to-2

10-s 2000

0

4000

6000

8000

10000

12006

14000

Time increments global time-step:

local time-step:

-

---

-

inviscid -

-

inviscid, reacting

inviscid inviscid, reacting

Fig. 3.2. Convergence history. Global vs. local time-step algorithm.

a I --_-_~~---__II I

I I

b /

4zq120 100

--_

80

--T---

T---T---T

60

-80

-60

-40

-20

T/FCC i-1

0

-400

-606

-200

4 [10F3 m]

Fig. 3.3. Temperature

0

z [10e3 m] globial time-step:

local time-step:

-

inviscid

---

inviscid, reacting

---

inviscid inviscid, reacting

distribution along stagnation stream line (a) and along body surface (b).

10

J. Argyris et ai., Physical and computationaJ aspects of reacting flows

1512 triangular elements 825 nodes

Fig. 3.4. Hypersonic flow past a wedge. Probiem description (a) and layout of finite element mesh (b).

with the two limiting cases of a frozen chemistry and a flow in chemical equilibrium. Figures 3.5 und 3.6 show profiles of temperature and density at the considered positions normal to the surface along the q-axis. As a result of the high temperature behind the shock wave, the dissociation process yields an increase of atomic species downstream of the apex. This increase of atomic species is accompanied by a lower temperature and thus by higher density. As a consequence, the normahsed distance Q/Z,-,, which represents the thickness of the shock layer, decreases, The approximation of the flow to a frozen state near the apex as well as to an equilibrium state far downstream is observed. The frozen flow as well as the flow in chemical equilibrium are independent of the wedge side length and qS/le is thus constant for any Ee. The investigation of a chemically reacting hypersonic flow over a 30°

J. Argyris et al., Physical and computational

aspects of reacting ffows

11

100 finite rate chemistry:

80 -

-

lo = O.lm

----

lo = 100m

---

ls=lm

----

lo = 1OOOm

.____._.

lo = 1Om

60 frozen chemistry

17110 [lo-s]

-------------____.______ 40 -

10

5

Fig. 3.5. Cross-sectional

15

20

25

31 D

T/Too 1-l temperature profiles along the wall normal.

100 finite rate chemistry:

80

-

lo = O.lm

----

lo = 100m

---

le=lm

----

lo = 1OOOm

____.--.

lo

=

10m

60 frozen chemistry

rlllo [lo-s]

0

5

10

15 T/Too

Fig. 3.6. Cross-sectional

20

25

1-l

density profiles along the wall normal.

half-angle cone shows comparable behaviour [6]. The dependence of the shock stand-off distance following geometric variations has also been studied for chemically reacting flows around spheres [7]. According to eq. (13), the geometric length of the problem directly affects the local stability limit of the flow whereas the local stability limit for the integration of the chemical kinetics is not sensitive to length scales. As a consequence, subincrementation of the chemical kinetics is required for the computation of the flow over the wedge with large side length in order to ensure stability. In contrast to this, in the flow over the wedge with small side length, the production rates are small and almost constant in time and thus must not be computed at every incremental step.

12

J. Argyris et d., Physical and computational aspects of reacting flows

The first approach to reduce the computational effort is based on the fact that the time scales of flow and chemical phenomena differ by orders of magnitude. An adaptive application of different solution and coupling strategies may reduce the computational effort. The second approach is based on the fact that the reacting flow area may be divided in a number of sub-areas with different chemical characteristics. While a particle proceeding along a stream line crosses a shock wave, a number of states with different chemical characteristics are traversed. In front of the shock wave, no chemical reactions occur due to the low farfield temperatures in the hypersonic regime. Behind the shock, a thin transition layer with frozen chemistry can be observed followed by an extended layer with non-equilibrium chemistry. As a consequence, the domain where non-equilibrium chemistry occurs is restricted. This domain may be determined during the course of the numerical integration by frequent analysis of the reacting flow field. In order to reduce the computational effort, the treatment of nonequilibrium chemistry may then be restricted to a certain domain with considerable chemical production rates. a) Adaptive application of different integration techniques for the chemical kinetics A non-dimensional parameter which compares characteristic time scales of flow and reaction processes is the Damkohler number [6], [2]. An alternative Damkijhler number which is related to the computational characteristics of reacting gas flows compares maximum stable time-steps for the numerical integration of the differential equations for the flow and the chemistry respectively. This number is defined as (23) The suitable solution method is indicated by an assessment of the reacting flow domain on the basis of a local estimation of the mesh Damkiihler number. The ratio of explicit integration steps for the flow and chemical kinetics is k = s DaA CFL

(24)

where DaA is the maximum Damkohler number in the flow field, C&n is the Courant number and s > 1 is a safety parameter (s = 1.25). Time measurement shows that, for k > 4, the implicit approach (eq. 23) is more efficient than the explicit approach (eq. 20) including subincrementation, whereas for k < 4, the explicit approach with subincrementation is most efficient. If the stability limit for the integration of the flow equations is the one associated with a restriction of the global time increment (k < l), the chemical kinetics may be performed with constant production rates for k,, = l/k - 1 incremental steps where k = l/3,1/4,. . .. b) Definition of the domain of chemical non-equilibrium In order to separate the domain where non-equilibrium chemistry occurs from the domain with frozen chemistry, the local convective mass fluxes are compared with the local production rates

(25) The factor sw represents a threshold value. A mesh point is incorporated rate chemistry if this threshold value is exceeded.

in the domain with finite

c) Application of adaptive techniques The aforementioned adaptive techniques are now studied for the inviscid flow past a wedge. The characteristics of this problem have been discussed before. As demonstrated in Fig. 3.5, the flow over

J. Argyris

et al., Physical and computational

13

aspects of reacting Aows

a wedge with a side length lc = 0.1 m tends to be in a frozen state concerning the chemical reactions. Typical mesh Damkijhler numbers in the flow field are of the order lo- 2. Thus, the stability requirements for the flow computation limit the global time-step. As a consequence the chemical production rates may be considered to be constant for a certain number of time-steps. In Fig 3.7.b, the specific time increments are marked where the production rates had to-be updated. It should be noticed that the stability requirements at initial state necessitate frozen chemistry for the first 100 time-steps. In the following, the number of time-steps with constant production rates drops from initially k,, = 100 to 73 at steady state. This decrease arises from the dropping stability limit for the explicit integration of the rate equations. This effect is typical for dissociation reactions at a high temperature [8]. Application of this technique reduces the computational effort by almost 10% compared to the unmodified method. As demonstrated in Fig 3.7.a, the convergence behaviour is not affected by this method. b

a

10-O

10-l

E i: B

E x

10-Z

.z

0”

6

10-S

10-d 400

600

Time increments

800

1000

r 0

5oc

81

Time increments

Fig. 3.7. History of convergence (a) and production rate updates (b)

As demonstrated in Fig. 3.5 and Fig. 3.6, the influence of the chemical reactions on the species distribution in the flow field grows with increasing length scales. This is explained with the increasing time that a flow particle requires to pass the flow field. During that time, chemical reactions will change the gas composition in the particle. According to the Courant-Friedrich-Lewy condition (13), growing length scales result in higher local stability limits for the flow computation. As a consequence, the explicit integration of the non-equilibrium chemistry commences to limit the time-step for the coupled process. Consequently, the stiffness of the numerical problem increases. For a side length Zc = 1 m, the number of time-steps with frozen production rates is diminished from initially 15 to 9 at steady state. As a consequence, the reduction in computational time diminishes slightly as compared to the previous case. For Zc = lOm, the mesh Damkohler numbers reach the order of one. The local stability limits for the flow computation and for the explicit integration of the chemical kinetics are now of the same order and the production rates have to be updated at almost every incremental step. A further increase of the length scale leads to a situation whereby the steady-state solution for the reacting flow tends towards an equilibrium state for the chemistry. This equilibrium state is far from the frozen state which is assumed to apply for the first 100 time-steps. As a consequence, the chemistry is characterised by high stiffness. For 20 = 100m three implicit subincrements are required in order to overcome this initial situation. The local gas composition is now strongly affected by the source terms. This results in strong variations of the energy in the flow field with the consequence that the reaction constants, which are functions of the temperature, have to be updated before every incremental substep. In the further course of the computation, the chemical kinetics are treated by explicit integration including

J. Argyris et al., Physkal

14

and computational

aspects of reacting flows

subincrementation. Due to a reduction of the stability limits, the number of subincrements has to be increased from 2 to 4 in the course of the computational time integration process. For lo = lOOOm, the number of implicit subincrements rises to thirty at the initial state. The mesh Damklihler numbers are now typically of the order 10 2. In the further course of the integration, the application of the linearised implicit method is as a rule more efficient with respect to the computational effort. Figure 3.8 presents a survey of the applied integration techniques.

I

lri 1

-

explicit solution method -

-

implicit solution method

a.) Ee = 0.1 m b.) lo = 1 m c.) Zo=10m d.) lo = lOOm e.) lo = lOOUrn

/_--_LI~---___~---____II______ i

ahc

I

I

200

400

I

I

600

I

I

800

I

I

1000

Time increments Fig. 3.8. Adaptive

application

of different solution methods for the Row over a wedge with variable side length.

The optimum solution method is determined by a repeated analysis of the reacting flow domain. The frequency of this assessment can be reduced during the course of the integration. For situations where the mesh DamkGhler numbers reach the order of one, the source terms have to be computed at every incremental step. A reduction of the computational effort is then obtained by limiting the domain which is characterised by chemical non-equilibrium. The limitation is based on a comparison of the local convective fluxes with the local source terms, according to (25). The threshold value sW is decisive for the extent of the reacting domain. As demonstrated by Fig. 3.9, the region of non-equilibrium increases in the course of the integration. Higher threshold values obviously result in a less extended domain of chemical non-equilibrium, For a t~eshold due essentially greater than s, = 0.1, the convergence behaviour (Fig. 3.10) as well as the steady-state solution is affected by the restriction. The threshold value s, = 0.5 has been proved in a wide range of appiications. An examination of the numerical solution every 10 time-steps proved sufficient to achieve a satisfactory adaption of the non-equilbirum area to the flow condition. According to the high effort for the computation of the chemical kinetics in flows with large length scales, this method becomes more effective, the more costly the integration method is. For the case la = 1 m, where the simple explicit intergration method is applied, the restriction technique results in a reduction of the computations time by 6%. A reduction of nearly 20% is obtained by an application of this technique on the lo = 1000 m case, where the linearised implicit approach has to be applied in order to ensure stability. This correlation is demonstrated in Fig. 3.12 on the baais of CPU-time measurements which have been performed for a side length of la = 1 m and 10 = 1OOOm. Here, the CPU-time which is required for the computation of the flow field and the chemical kinetics is normalised by the CPU-time required for gather-scatter operation which are necce~ary in finite element computation. They are performed separately because they are not required for computations te~niques based on different spatial discretisation as the finite difference and the finite volume method [9].

J. Argyris et al., Physical and computational

aspects of reacting flows

15

50

200

400

600

800

1000

Time increments Fig. 3.9. Node history for different threshold values.

Variation of threshold value sw: .-- ---..

10-4

I 0

& TO.1

200

---

400

s, = 1.0

600

800

1000

Time increments Fig. 3.10.

Convergence

history for different threshold values.

Four test cases have been considered. The first test case is a frozen flow, where exclusively the convective transport of the conservation quantities is considered. The second test case represents a flow where finite rate chemistry is considered in conjunction with convective transfer. While the explicit integration of the chemical kinetics for lo = 1 m requires 12% more CPU-time, the additional effort for 2s = 1000 m rises to 30%. This is explained by the fact that, for the second case, the more costly linearised implicit method is required. Test cases 3 and 4 represent measurements done at two different times during the course of the integration if an adaptive restriction of the reacting area is considered. Test case 3 represents the situation after 101 time increments, where the reacting flow field is not fully developed. Test case 4 represents

J. Argyris et al., Physical and computational

16

aspects of reacting Aows

200 Variation of threshold value sw:

150 -

_.__..__

sw =O.l

-

SW = 0.5

---

SW = 1.0

9110

[lo-“]

25

0

T/Too L-1 Fig. 3.11. Temperature

profiles for different threshold values. b

a

‘.‘.

“““’

:::: :.:. :: :: :: ::

Iil :::: 1.:.

:: :; ::

. ...i. ::::

::::

:::: ‘.‘.

::::

‘.‘. :.:. A:

0.00

:.:.:.: :;:j:;:

,.,. ::::::: (1.1 ;:;:j:j

:

1

2

gather, scatter operations

3

4

test case

test case

0

2

q

flow field

Fig. 3.12. Time measurements

q

chemical kinetics

for adaptive techniques.

the situation near steady state after 1100 time increments. It should be noted that the application of the adaptive technique leads to a reduction of the computational effort for Ze = 1 m from initially 8.4% to 5.6% at steady state. In contrast to this, the computational effort for lo = 1OOOm is reduced from initially 17.5% to 13.3% at steady state.

4.

Chemical reactions in viscous hypersonic flows

In order to account for viscous effects, the partial differential equation system (1) for the flow has to be extended. Accordingly, a system of equations governing the hypersonic flow of a viscous gas is provided

J. Argyris et al., Physical and computational

by the compressible Navier-Stokes

au

equations expressed in the conservation

Wl + 91)+ w2 + 92)+ w3 + 93)=

at+

ax2

3x1

aspects of reacting Aows

17

form

(26)

0.

3x3

The additional vector arrays gi, g2, g3 comprise the terms arising from the viscous phenomena in the flow field. Assuming a chemically reacting multi-component gas mixture, the caloric and thermal properties have to be determined for the local composition of the gas, which is governed by the condition of conservation of species. Considering mass flux by diffusion, the equations for species conservation may be collectively written as in [2]

aw

a(211w+jl) + w2w+j2)

at+

ax,

ax2

where ji represents the diffusive fluxes. solution of 5 + N conservation equations. single matrix expression

+ ah3w+j3) = w. ax3

(27)

Assuming finite rate chemistry, the flow problem requires the In analogy to (3), the overall system may be assembled into the

A number of N, species conservation equations may be replaced by the respective equations for the conservation of elements, as outlined previously in section 2. Obviously, the extensions introduced in (28) by the diffusive fluxes ji and the viscous terms a do not affect the global solution procedure as described in section 3 for the inviscid flow. Merely the diffusive mass flux and the viscous contributions have to be determined in addition to the quantities of the inviscid Euler case. The viscous terms g comprise deviatoric stresses which are related via the dynamic viscosity ,u to the deviatoric rate of deformation, conductive heat flux which is related via the thermal conductivity k to temperature gradients and diffusive mass flux which is related via diffusion coefficients Drj to the concentration gradients. Mass diffusion caused by temperature and pressure gradients are neglected. The reader should recollect that the transport coefficients of a multicomponent mixture depend in general on the temperature, the pressure and the composition of the mixture. These coefficients may therefore be formally written as CL=

P(T,P7 sr> (29)

k = k(T,p,Er) Drj = Drj (T,P, Cr) Some appropriate physical models for the molecular transport

coefficients have been introduced in [2].

5. Improved solution methods for viscous chemically reacting flows 5.1. Local time-stepping

for steady-state

Aow problems

In principle, the technique of local time-stepping which has been introduced in section 3.2 is also applicable to viscous flow problems. In the context of the Taylor-Galerkin method, however, differences in the convergence behaviour as well as in the steady-state solutions arise. The behaviour of the local time-stepping technique applied to a viscous flow problem will be studied here for the chemically reacting flow over the 30° wedge. The characteristics of this flow problem have been discussed in section 3.4 for the corresponding inviscid case. In addition to the boundary conditions considered in the inviscid case, the velocity of the gas normal and tangential to the surface must be suppressed. Furthermore, the temperature of the gas on the surface of the body is required as a boundary

18

J. Argyris et ai., Physicsaf and computational aspects of reacting fiows

condition in the numerical procedure. In the following, the temperature of the gas at the surface is assumed to be T, = 1500K. The finite element mesh has to be adapted to the req~rements of viscous Bow computation. The flow field over the wedge has been discretised by a mesh of 3584 triangular elements, shown in Fig. 5.1. The triangular elements are mostly concentrated near the solid wall in order to achieve a sufficient resolution of the boundary layer.

1885 nodes 3584 triangular elements

Fig. 5.1. Hypersonic viscous flow past a wedge. Layout of finite element mesh.

In order to study the behaviour of the local time-stepping technique the viscous as well as the inviscid reacting flow over the wedge have been computed on the same finite element mesh. The history of the density residua presented in Fig. 5.2 clearly demonstrates that the local time-stepping procedure improves the convergence of the numerical algorithm for the inviscid flow. An improvement of the convergence The solution, however, converges to an rate is shown also in the initial stage of viscous computation. oscillatory flow field which cannot be corrected sufficiently as indicated by the high residuum. In order to compare the numerical solutions, cross-sectional profiles of the normalised temperature were taken at 50% side length normal to the wedge surface. Fig. 5.3 confirms that the solution method does not significantly affect the steady-state solution in the inviscid case. On the other hand, considerable differences are observed in the viscous case. The maximum temperature in the inviscid flow field is unaffected by the local time-stepping technique and arises at the same position and varies less than 0.5%. In the viscous solution, the m~mum temperature is observed at different positions in the flow field and differs by nearly 10%. These variations are also visible in the cross-sectional profiles plotted in Fig. 5.3. The above discussion does not support the application of the local time-stepping in viscous flow computations. 5.2. Adaptive ~~~~cat~o~ of physical operators The inclusion of viscous effects computational effort. Additional contributions in the finite element affects the flow field are restricted

in the numerical treatment of a flow problem requires an enlarged operations have to be carried out in order to establish the diffusive equations. In general, the areas where the diffusive transport mainly to near wall regions.

J. Argyris et al., Physical and computational

19

aspects of reacting flows

reacting flow,

-

local time-step

Time increments Fig. 5.2. Convergence

history for global vs. local time-step

algorithm.

inviscid, reacting flow: - global time-step

-

--

-

- - - local time-step

viscous reacting flow: -

global time-step

-

local time-step

QllO

[lo-s]

loo -

50 -

00 0

10

5

T/Tm Fig. 5.3. Temperature

20

15

25

3

L-1

profiles normal to the solid wall for global vs. local time-step

algorithm.

In order to reduce the number of numerical operations, we propose an adaptive approach. By a repeated analysis of the viscous flow field, certain regions are detected where the diffusive transport has to be considered. An analysis of the viscous flow field is carried out on the element level by comparing the convective and the diffusive flux contributions. An element has to be entered into the viscous domain if the ratio of viscous flux to the convective flux exceeds a threshold value s/, in any equation m

> SF. Ilagm,ill I IIa+1’2fm,i((

(30)

20

J. Argyris et al., Physical and computational sspects of reacting flows

Application of the adaptive technique In what follows, we investigate the effect of the above adaptive technique on the viscous flow field. Once more, the subject of our studies is the chemically reacting flow over the 30” wedge. The threshold value sP determines the extension of the domain in which diffusive flux has to be considered. The evolution of the viscous region is illustrated in Fig. 5.4. 50

_f____._~-~_~_‘---~-.~.~_-__ .J ,Y #J

variation of threshold value scL: s/, = 10-2 for monotonic growth of the viscous domain: ---. sp = 10-3

.r II i

=

10-Z

-

sp

----

SW= 10-l #

.-_--_--_--_.._-.--.--.-

0

I 0

I 2000

I

I

I

4000

I 6000

1

I

8000

I 10000

Time increments Fig. 5.4. Extension of the viscous domain with respect to different threshold values.

10-O

10-l

variation of threshold value sjL: -.-. SbJ-0 10-Z

10-s

70

2000

4000

6000

8000

1( DO

Time increments Fig. 5.5. Convergence

histogram

with respect to different threshold values.

J. Argyris et al., Physical and computational aspects of reacting flows

21

In the course of a time integration, it has been observed that in a state in which the viscous region reaches its maximum extension elements placed at the boundary between the viscous and the inviscid region tend to oscillate between both regions. In order to prevent such oscillations, elements once placed in the viscous region are not transferred back into the inviscid region. This ensures a monotonic growth of the viscous area.

The consequences of element oscillations for the convergence behaviour is demonstrated in Fig. 5.5. Beside the element oscillations the use of unsuitable high threshold values for the separation between viscous and inviscid computation also affect the convergence behaviour. The effect of an excessive reduction on the numerical solution is demonstrated for the temperature in Fig. 5.6 and for the tangential component of the velocity in Fig. 5.7.

200

_.-_

= 0

for monotonic growth of the viscous domain: -.-* sp = 10-3 --_-

SC” =

10-l

0 0

5

10

15

20

25

30

T/‘Tm f-1

Fig. 5.6. Temperature profiles normal to the solid wall.

For threshold values lower than al.L= 10W3, the convergence behaviour as well as the steady-state solution is not affected by the adaptive technique. The ratio of the viscous area rises from initially less than 5% to 43% at steady state. If the threshold value exceeds s,, = lo-“, the adaptive technique affects the convergence behaviour. Marked differences in the numerical results are obtained if the threshold value exeeds sfi = 10-l. Nevertheless, the velocity and the temperature gradients at the solid wall are not affected. Thus, the skin friction and the wall heat flux are not sensitive to the adaptive technique. An analysis of a solution incorporating viscosity every 100 time-steps proved sufhcient in order to achieve a satisfactory adaption of the viscous area to the actual flow condition. Finally, the efficiency of the adaptive application technique is discussed on the baais of CPU-time measurements. Fig. 5.8 shows an survey of the considered test cases. The CPU-time required for the computation of the convective contributions and the chemical kinetics and for the diffusive contributions is normalised by the CPU-time required for gather-scatter operations. The first test case represents a chemically reacting inviscid flow. The integration of the flow field and the chemical kinetics requires considerable effort. The introduction of diffusive transport in test case 2 leads to a growth of the computational effort by 17%. Application of the adaptive technique in test case 3 with a threshold value of So = lop2 reduces the additional effort to 5% as compared to the inviscid case.

22

J. Argyris et al., Physical and computational

aspects of reecting flows

variation of threshold value sM: -.-.

8w = 0

for monotonic growth of the viscous domain: -._. sfi zz 10-3 a..-

SP

z

10-l

50 -

0

1 5

0

I 15

I 10 +-kc

I 20

5

i-1

Fig. 5.7. Profiles of the tangential velocity normal to the solid wall.

..,..

:.:.. ::. :;.

..$j;:...

.,

.‘.‘. ::. ::_ ::. ::, .‘.‘, ::. ::. ::. .‘.‘. .‘.‘. ,‘.’ ::, ::.

~

.‘.‘.

1

~ 2

q q

convective transport and chemistry

q ,,,:,

diffusive contributions

gather, scatter operations

3

test caSe Fig. 5.8. Time measurements for adaptive techniques.

6. Physical models for chemically reacting viscous flows The characteristics of a chemically reacting boundary layer will be demonstrated by the study of the flow over a flat plate. The flow medium is taken to be atmospheric air. The farfield Aow is characterised by the Mach number A&., = 10, the temperature T, = 275.6 K and the pressure p, = 4150 Pa. The flow domain is covered by a mesh consisting of 7040 triangular elements, which are mostly concentrated near the solid wall (Fig. 6.1). In order to study the characteristics of this problem, sections through the flow field were taken at certain positions vertical to the solid wall. It proves helpful to define the position of the vertical cut in a normalised manner as

In the following, cross-sectional profiles of some representative

thermodynamic

variables are presented.

J. Argyris et al., Physical and computational

3649 nodes 7040 triangular

Fig. 6.1. Hypersonic

23

aspects of reacting ffows

elements

flow over a flat plate. Layout of finite element mesh.

In addition, the variation of these variables due to different models of the chemistry is studied. Solutions obtained by application of an equilibrium chemistry solver [2] are traced with thin lines while solutions derived by finite rate chemistry are traced with thick lines. In Fig. 6.2, profiles of the local pressure are presented, The downstream evolution of the flow displays a smearing tendency of the shock due to the coarsening of the upper part of the mesh. The assumption of chemical equilibrium leads to a slight decrease of the shock distance to the solid wall and to lower pressure values in the shock layer.

20

I

‘\ ‘1. [10y2

m]

lo

1.0

1.2

1.4

-

&fG$ = 1000

---

G=

--------

&

= 2000

----

6

= 2500

----

x

= 3000

1.6 P/Pm

1500

1.8

2.0

i-1

Fig. 6.2. Normalised pressure profiles for non-equilibrium

vs equilibrium chemistry

at different positions.

24

J. Argyris et al., Physical and computational

aspects of reacting flows

The variations obtained for different models of the chemistry will become more pronounced for thermodynamic variables which are more sensitive to the chemistry. Note, that the following figures exclusively represent results within the boundary layer. The temperature distribution in the flow field, as demonstrated in Fig. 6.3, is strongly affected by the model applied for the chemistry.

15 -

-

&

= 1000

---

&=1500

-------w &=200(-j ---&= 2500

[UC3 m]

lo -

5-

0

6

3

T/Too Fig. 6.3. Normalised temperature

12

9

15

1-I

profiles for non-equilibrium

vs equilibrium chemistry

at different positions.

20* -

JiTG = 1000

---

-=

1500

_.*.--*-

-=

2000

----

+

_---

s=

I 0

2

4

= 2500 3000

I

I

6

L

8

50 1x1 Fig. 6.4. Profiles of atomic oxygen for non-equilibrium

vs equilibrium chemistry

at different positions.

J. Argyris et al., Physical and computational

aspects of reacting flows

25

The assumption of local chemical equilibrium results in a remarkable reduction of the temperature. This is explained by the advanced dissociation of oxygen. As plotted in Fig. 6.4 the msss fraction of atomic oxygen is much reduced if non-equilibrium chemistry is considered. The evolution of the boundary layer is also affected by the modelling of the chemistry. The establishment of chemical equilibrium results in a slight decrease in the boundary layer thickness. An important difference between inviscid and viscous flows concerns the specification of particular boundary conditions for the latter. In the following, a brief survey of appropriate solid wall boundary conditions is presented. 6.1. Solid wall boundary conditions In the case of a viscous Aow over a body, the medium sticks on the surface of the body and the velocity of the gas normal and tangential to the surface must be suppressed. In the presence of chemical reactions and mass diffusion within a boundary layer, the influence of the surface on the gas composition, denoted as catalycity, becomes significant. Due to the complex physical nature of such interactions with the surface, one has to specify simplified models. In order to assess the influence of surface catalycity on the flow, a simple model which has been applied with some success in [lo], [ll] will be considered in conjunction with two surface models which are widely used as boundary conditions in hypersonic flow computations. a) Surface catalycity The catalytic behaviour of a solid wall is characterised by the rate at which atomic species are recombined at a certain rate. The consumption of atoms due to surface reactions is balanced by the diffusive mass flux to the wall. Thus, the following boundary conditions for the mass transfer may be specified

-gas, ay

h/i,&,s x0.

Here, n/i, denotes the molecular weight and w,” the production rate of species r at the surface. If an adiabatic wall is considered, the heat transported by conduction must be compensated by the energy flux due to mass diffusion. The following boundary conditions for the heat transfer may be specified kg

+ c

M,ti,Sh, = 0. N

One limiting case is characterised by a negligible influence of the wall surface on the gas composition. Thus, no recombination occurs at the surface (wf = 0). As a consequence, following eqns. (31) and (32), all species derivatives as well as the temperature gradients normal to this perfectly non-catalytic wall vanish. A perfectly catalytic wall represents the second limiting case. Recombination as well as dissociation occurs spontaneously at an infinite rate. As a consequence, local chemical equilibrium is established at the wall (cf = (
26

J. Argyris et al., Physical

and computational

aspects

of reacting

flows

.

In (33), Ai stands for the chemical symbol of the i-th participating substance, and v[~, v:; are the stoichiometric coefficients in the j-th reaction appertaining to the i-th substance. Any reaction j may occur with a probability rj if the species Ai sticks to the surface. The surface recombination rate Lj,S is then given by the expression

(34)

In (34), R is the universal gas constant and m, the total number of considered surface reactions. As a consequence, we deduce from (31) that the derivatives of the species cannot vanish normal to the partially catalytic wall. Furthermore, an adiabatic wall with non-zero temperature gradients is established in conformity with (32). Starting with the assumption that an air model comprising the 5 species 0, N, NO, 02, Ns is used to describe the chemical kinetics in the gas phase, an appropriate model for the surface chemistry is given by

The recombination coefficients yo, TN are in general complex functions specific material data. For a detailed treatment of surface reactions, be considered. At first, the adsorption of atomic species to the surface recombination at the surface reduces the concentration of atomic species. and atoms from the surface into the gas phase takes place. In order to study the effect of surface catalycity, the recombination considered to be equal and constant.

of the temperature and various a number of processes have to has to be described. Catalytic Finally, desorption of molecules coefficients yo and YN may be

c) Surface radiation In stagnation regions and shock layers of hypersonic flows, the temperatures rise to such values that radiation contributes significantly to the energy balance. As a consequence, the flow field becomes nonadiabatic and total enthalpy variations also arise in inviscid regions. In addition, one has to take into account the radiative heat transfer to a solid wall besides the ordinary convective heat transfer. The treatment of radiative gas dynamics is quite a complex problem. If the air is assumed to be transparent, the adsorption and emission of radiative energy is negligible within the flow field. Thus, the radiative heat transfer is limited to the solid wall. The radiative heat transfer influences the boundary conditions. While the boundary condition for the mass flux (31) remains unchanged, an additional term concerning the radiative heat flux has to be considered in the energy balance. The extended formulation for the energy balance at a radiating wall then reads

In (35), E denotes the emissivity of the surface material and D the Stefan-Boltzmann constant. consequence, following eq. (31), the species derivates cannot vanish normal to a partially catalytic According to (35), a radiation adiabatic wall is characterised by non-zero temperature gradients.

As a wall.

d) Hypersonic flow over a flat plate In the following the effect of the considered conditions for a solid wall boundary is discussed for a viscous boundary layer. First, some studies concerning the influence of surface catalycity on the flow field

J. Argyris et al., Physical and computational

reacting flows

aspects of

27

presented. A discussion about the effect of surface radiation follows. Finally, the combination of both boundary conditions will be taken into account. The influence of surface catalycity rises with growing recombination probabilities. Typical values for the recombination coefficients rj vary within the range of 10V3 to 10-l. Figure 6.5 demonstates the effect of surface catalycity on the temperature and mass fraction of atomic oxygen in the boundary layer. In addition, the corresponding results for the two limiting cases of a perfectly catalytic and a perfectly non-catalytic surface are presented together with the solution which is obtained if equilibrium chemistry is considered. The cross-sectional profiles were taken at a normalised position & = 2000.

are

b

a

20 I I I

I I I I

I I

I

I I I I

I I I

I I I

I I

I I I

-_--+---+---+_--+_--_-

I

I I I

I I I

I I I

I I I

2

4

6

8

I

I

15

-+---+---+_--+_-__-

lo

ilOT

m]

I

0

3

6

9

12

15

0

T/Too [-I nonequilibrium

Eo [%I

chemistry:

-

noncatalytic

---

partially catalytic

-

--

-

- --

10

equilibrium -

wall

- partially catalytic partially catalytic

- - - - - - - - perfect catalytic

chemistry:

noncatalytic

wall

wall y = 10F3 wall 7 = low2 wall 7 = 10-l wall

Fig. 6.5. Profiles of normalised temperature

(a) and atomic oxygen (b) for different models of catalytic

surfaces.

The catalytic atom recombination leads at first to a reduction of atomic species at the catalytic surface. The production of molecular species and the consumption of atomic species is compensated by mass flux to the surface. As a consequence, the species gradients do not vanish normal to the solid wall. The increase of the temperature at catalyic surfaces is explained by the exothermal nature of the recombination reactions. The higher temperature results in a flux of energy into the flow field. As stated before, the temperature gradients are non-zero normal to the solid wall. As expected, higher recombination probabilities lead to higher recombination rate and thus to increased consumption of the atomic species. For yo = TN = 10-l, nearly all atomic oxygen vanishes at the surface. The equilibrium condition for the chemistry in the case of a perfectly catalytic surface causes progressing dissocitation. As a consequence, the surface temperature drops to the value obtained by application of chemical equilibrium solver. We notice that the perfectly catalytic wall does not, however, describe the physical properties of technical surfaces. The variation of the surface temperature and atomic oxygen along the plate is plotted in Fig. 6.6. Due to dissociation, the concentration of atomic oxygen increases downstream of the leading edge. Catalytic atom recombination compensates the dissociation process to a certain extent. As a consequence, the temperature rises with increasing influence of the surface catalycity.

J. Argyris et al., Physical

28

T/Tee

computational

aspects of reacting flows

a

[-I

20 I I I

I I I

I I I

I I I

I I I

1.5

1.0

0.5

I I

I I

1.5

2.0

I I

10 0.0

0.5

1 I

1.0

1

I

I

I

I 2.5

0.0 3.0

0.0

0.5

VGG PO31 nonequilibrium

1.5

2.0

2.5

: II

VQG [lo31

chemistry:

noncatalytic

--__-_ ____.----__

1.0

equilibrium

noncatalytic

wall

partially catalytic

wall y = 10d3

partially catalytic

wall y = lo-’

partially catalytic perfect catalytic

chemistry: wall

wall y = 10-l wall

Fig. 6.6. Distribtion

of temperature

(a) and atomic oxygen (b) at the surface.

The ability of a surface to radiate a substantial amount of energy is quantified by the emissivity E. This material property may vary within a range of 0 6 E < 1 depending on the surface material. For a number of materials used for the heat shields of hypersonic vehicles E = 0.8, is a reasonable value [6]. Fig. 6.7 demonstrates the effect of surface radiation on the temperature and mass fraction of atomic oxygen in the boundary layer. Results for non-equilibrium chemistry in connection with a radiating surface are compared with the corresponding results obtained with the equilibrium chemistry solver. Again, the cross-sectional profiles were taken at a position E = 2000. The radiative heat transfer at the surface results at first in a remarkable reduction of the temperature in the boundary layer and at the solid wall. As expected, higher emmission coefficients lead to a stronger cooling of the wall. Because surface catalycity is not considered in the present case the transport of energy to the surface has to be balanced by conduction in order to keep the surface adiabatic. As a consequence of the lower temperatures in the boundary layer, less dissociation is observed. The atomic species vanish completely at the wall if local chemical equilibrium is considered within the flow field. Finally, the aforementioned boundary conditions may be combined in order to estimate their relative importance. In Fig. 6.8, the effect of surface radiation and catalycity on the temperature and mass fraction of atomic oxygen is demonstrated. All results are obtained considering non-equilibrium chemistry. The flow field in the boundary layer of a radiating wall is compared with the corresponding results for a partially catalytic and a perfectly catalytic wall. Again, the cross-sectional profiles were taken at a position & = 2000. Obviously, the radiative heat transfer has a strong effect on the energy distribution within the boundary layer. Concerning the surface temperature, the simulation of the surface catalycity is less important. Only for the equilibrium boundary condition which has to be considered for a perfectly catalytic wall results in a reduction of the temperature. On the other hand, the species distribution is more affected by the influence of the catalytic surface. Thus, the fraction of atomic oxygen is negligible at the wall if catalycity is considered in combination with surface radiation.

J. Argyris et al., Physic& and computational

aspects of reacting ffows

a

b 20

12

0

15

,

I

I

,

0

2

4

6

l-1

T/T,

to

nonequilibrium chemistry:

equilibrium chemistry:

-

nonradiating wall

-

radiating wall, E = 0.1

-

-

-

-

----

29

8

WJ,]

nonradiating wall

- radiating wall, E = 0.8

--

radiating wall, E = 0.8

Fig. 6.7. Profiles of normalised temperature

(a) and atomic oxygen (b) for radiating and noradiating surfaces.

a

I

b

/

I

I 20

20

5

15

I

__-+--_ i

lo

3

6

9

12

f--_+---+-__ !

!

!

!

0.2

0.4

0.6

0.8

[10y3 m]

0.0

15

so 1%

‘-7’/7’00f-1 nonequilibrium chemistry: -

noncatalytic,

-

-

-

nonradiating wall

noncatalytic, radiating wall, s = 0.8

----

partially catalytic, radiating wall, 7 = lo-‘?

-

perfect catalytic, radiating wail, E = 0.8

---

10

Fig. 6.8. Profiles of normal&d

temperature

E = 0.8

(a) and atomic oxygen (b) for catalytic and radiating surfaces.

-I

1.0

30

J. Argyris et al., Physical and computational

6.2. Effect of simplified description of transport a)

aspects of reacting Bows

coefficients

Similarity parameters for the transport properties

For the numerical treatment of the viscous flow of multicomponent gas mixtures with chemical reactions, the molecular transport coefficients of (29) have to be provided. A number of models describing the dependence of the transport coefficients on thermodynamic and gas kinetic properties have been developed in the past. The employment of these models for the numerical treatment of a flow problem requires additional computational effort because of the complex structure of the expressions. Therefore, it seems to be advantageous to replace the formulations by simplified expressions which make use of similarity considerations. The consequences of the simplified modelling for the flow field have to be investigated. A suitable model for the transport coefficients which has been applied in hypersonic flow computations is the Straub model [2]. It has been established especially for the modelling of the molecular transport coefficients of a 0, N, NO, 02, N2-mixture as arising in atmospheric air. It provides specific expressions for the dynamic viscosity, the thermal conductivity and polynary diffusion coefficients. In the following, the Straub model will be compared with two models which have been documented in the problem description of the European Hypersonic Data Base (EHDB) [13]. In the first model (model l), the dynamic viscosity of the particular species is expressed as a function of the temperature. The dynamic viscosity of the mixture is obtained via Wilke’s mixture law [ll]. The thermal conductivity and the polynary diffusion coefficients are obtained via similarity expressions. The second model (model 2) provides expressions for the dynamic viscosity, the thermal conductivity and an averaged polynary diffusion coefficient. As reference, the viscosity may also be approximated by the simple Sutherland law

PI.

A comparison of the different models for the the dynamic viscosity coefficient is presented in Fig. 6.9. The mixture is assumed to be in chemical equilibrium at a fixed pressure p = 5 . lo4 Pa, but varying temperature. Below T = lOOOK, where no real gas effects are apparent, all approximations agree well. With increasing temperature, the Sutherland law provides too low values as compared with the more complex models. Within the considered temperature range, the complex models display scatter within a range of 6%. 0.3

0.2 P

[low3 Pa s] 0.1

equilibrium composition

-

0

2

4

6

8

:

Straub model

- - - - - - - - Sutherland 0.0

at p = 5000 Pa

----

Model 1

----

Model 2

law

10

T [lo3 K] Fig. 6.9. Dynamic viscosity of dissociating

air in chemical equilibrium

The simplest model for the thermal conductivity can be established via the specification of a similarity parameter. In the dimensionless Navier-Stokes equations, the Prandtl number relates the energy dissipated by friction to the energy transported by thermal conduction. It is defined by Pr = ~c~lc-~.

(36)

J. Argyris et al., Physical and computational

31

aspects of reacting Aows

Assuming a constant Prandtl number, the thermal conductivity k can be computed from the dynamic viscosity CL.In Fig. 6.10.a, the temperature dependence of the thermal conductivity is presented for the considered models. The calculation of the thermal conductivity via the similarity expression (36) results in a deviation below 4% as compared to the complex model. From (36), the Prandtl number can be deduced if the dynamic viscosity and the thermal conductivity are provided. Fig. 6.10.b shows the temperature dependence of the Prandtl number for the Straub model and model 2. If the Straub model is applied, Pr = 0.72 is a reasonable assumption. k [J/msK]

a

Pr

b

1-l

0.6 I I I I

I\

I I I ,p----r--_,

0.7 -

! 0.6 -

i

#f1

i; I I I I I I 2

1

-

-

-

composition

I I I I

I

I

I I I I I

I I_ I

I I 1 ,..-p-----

I I

1I

I I

I

I

I

1

I

I

I I I I

I I I I

I I I I

I 4

’

I 6

’

I 8

’

T [lo3 K]

T [lo3 K]

equilibrium

I I I I

at p = 5000 Pa

:

Straub model

----

Model 1, Pr = 0.72

Straub model, Pr = 0.72

----

Model 2

Fig. 6.10. Thermal conductivity

(a) and Prandtl number (b) of dissociating

air in chemical equilibrium.

The specification of diffusive mass fluxes allows the introduction of a second similarity parameter which relates the convective flux of species to the diffusive flux of the species. It is known as the Schmidt number and is defined by SC = &D)_l.

(37)

In some cases, the Schmidt number may be assumed to be constant and then it is again sufficient to provide only the dynamic viscosity ,u. Fig. 6.11.a demonstrates the temperature dependence of the diffusion coefficients. For the Straub model, we demonstrate the polynary diffusion coefficients quantifying the molecule-molecule and the atom-atom diffusion. The calculation of diffusion coefficients for the mixture via the similarity expression (37) leads to mass averaged values. They agree well with the values obtained from the expression provided by model 2. From (37), the Schmidt number can be determined if the dynamic viscosity and the mixture diffusion coefficient are provided. Fig. 6.11.b shows the temperature dependence of the Schmidt number for model 2. Notice that SC = 0.70 is a reasonable assumption for problems in the temperature range currently considered. c) Effect of simplified description on a boundary layer The consequences of the simplified description of the transport coefficients for the flow field as well as the an assessment of the savings in computational time will be investigated for the example of the hypersonic flow over the Aat plate which has been introduced before.

J. Argyris et al., Physical and computational

32

aspects of reacting fiows

a

2) [m2s-l]

L-1

Sc

b

0.8

10-o

I I 1

0.7 10-l

*A---.. !’

I/’

-

A”

/‘.

0.6

I I ,I

I

I

2

0

‘:

I I I /

I I

I I

/

I

Y.. _..-_I---.__ ’

I ,I

I

4

I I I ,

1

6

I I I, 8

I 10

T [IO3 K]

10-Z -

Straub model

-

Straub model, SC = 0.72

----

Model 1, SC = 0.70

----

Model 2

10-S 2

0

4

6

10

8

T [lo3 K] Fig, 6.11. Polynary diffusion coefficients

(a) and Schmidt number (b) of dissociating

air in chemical equilibrium.

In order to investigate the effect of a simplified description of the thermal conductivity exclusively, mass diffusion is neglected at first. Fig. 6.12 presents cross-sectional profiles for the temperature and atomic oxygen at a normalised position a = 2000. The profiles represent the results for the various models which have been introduced before. For the Straub model, the solution for equilibrium chemistry is incorporated. a

b 20

I

I

I I I I I I ___-+_--+---+---+----_~,rj I I

I I

I I I I

I I I I

I I

I I

20 I I I I ~~_~~~~_+~~~~+~~~~~_~~_ I I I I [IO!3

m]

3

6 T/Too

nonequilibrium ---

9

12

15

0

5

I I I I

I I I I

/

I

I I I

I I I

10

15

Eo I%1

L-1 chemistry

I I I

+--_-+_---+---/

lo------

I I I ~-___-_+--_-+_--_+-_--

0

I I I I

:

Straub model Straub model, Pr = 0.72

equilibrium -

-

-

chemistry:

Straub model Straub model, Pr = 0.72

- - - - - - - - Sutherland law, Pr = 0.72 ----

Model 1, Pr = 0.72

----

Model 2

Fig. 6.12. Effect of simplified modelling of transport

properties on the temperature

(a) and atomic oxygen (b).

J. Argyris et al., Physical and computational

33

aspects of reacting Aows

Within the boundary layer, significant deviations of the temperature profiles are obtained for the Sutherland law in comparison with the more complex models. The temperature distributions, obtained by an application of these more complex models, vary over a range of 3%. These variations are concentrated to the area close to the solid wall, where the flow is dominated by the diffusive transport. It is remarkable that the application of similarity parameters does not affect the solution. In a second investigation, the consequences of a simplified modelling of the diffusion coefficients for the flow field have been studied. Again, cross-sectional profiles for the temperature (a) and atomic oxygen (b) are presented in Fig. 6.13. The profiles represent the results for the various models which have been introduced before. For the Straub model, the diffusion matrix has been replaced by a single diffusion coefficient which is obtained from the similarity expression (37). The corresponding solution for equilibrium chemistry is again incorporated. a

b

I

, 20

20

15

15

(

1

3

6 T/Too

nonequilibrium

9

12

15

0

chemistry:

!

~

6

8

[%I

equilibrium chemistry:

---

Straub model, Pr = 0.72, SC = 0.72

-

Fig. 6.13.

4 Eo

-

----

2

L-1

Straub model, Pr = 0.72

-

-__+_-_+_--+__-+-! !

Straub model, Pr = 0.72 -

-

Straub model, Pr = 0.72, SC = 0.72

- - - Model 1, Pr = 0.72, SC = 1.00 Model 2 Effect of simplified modelling of transport

properties

on the temperature

(a) and atomic oxygen (b).

Once more, the solutions, obtained by application of the more complex models, show variations in the temperature profiles within a range below 3% which are concentrated to the area close to the solid wall. The reduction of the polynary diffusion matrix to a single mixture averaged diffusion coefficient has some consequences for the agreement of the solution obtained by the Straub model. The consequences of the simplified modelling concerning the computational effort are discussed on the basis of CPU-time measurements. Fig. 6.14 shows the survey of the results obtained for the reduced modelling of the thermal conductivity (a) and the polynary diffusion coefficients (b). The CPU-time which is required for the computation of the flow field and the chemical kinetics is normalised by the CPU-time required for the gather-scatter operations. The first test case represents a chemically reacting inviscid flow. The incorporation of diffusive transport in test case 2 results in an increase of the vector operations of nearly 50% and 70% if mass diffusion is considered. The introduction of similarity parameters leads to no significant reduction of the CPU-time, as demonstrated in test case 3. The gain of the simplified modelling is small in comparison to the effort, which is neccessary to account for the viscous contributions. Solely the introduction of a very simple material law like Sutherland’s law may reduce the computational effort significantly. This has been demonstrated in test case 4. Corresponding CPU-time measurements

34

J. Argyris et al., Physical

and computational

aspects

of reacting

flows

have been undertaken for model 1 and 2. Again, the reduced modelling of the transport coefficients yields no significant reduction in the computational effort.

1.1.

IL :::: :::: :::: :::: ‘.

::::

:::: , ,, .., ‘.‘.

:::: ‘7. ‘.‘.

::::

:::::::

jjl$$ :::: ;;;;;j; ::::

1:::::: ::::::: :::. 1.1. ::::::: 1.:.

4

test case u

gather, scatter operations

1

3

2

test case

q

diffusive contributions

q

convective transport and chemistry

Fig. 6.14. Effect of reduced modelling of transport properties on computational time.

Acknowledgment The authors would like to express their appreciation of the support given for this research within the R&D and the RfQ Program for the Aerodynamics of the European Space Project “Hermes” according to contract ATP-H/CP No. 3914/91 with Dassault-Aviation. Also acknowledged is the imaginative support given by Dornier Luftfahrt GmbH, Aerodynamics Department BF.

References J. Argyris, I.% Doltsinis, H. Friz, Hermes Space Shuttle: Exploration of Reentry Aerodynamics, Comp. Meth. Appl. Mech. Eng. 73 (1989) 1-51. PI J. Argyris, I.St. Doltsinis, H. Friz, J.Urban, An Exploration of Chemically Reacting Viscous Hypersonic Flows, Comp. Meth. Appl. Mech. Eng. 89 (1991) 85-128. L. H. Schafer, Untersuchung aero-thermodynamischer Vorglnge beim Wiedereintritt einer Raumfahre in die [31 Erdathmosphlre unter Beriicksichtigung von Realgas-Effekten, Diploma Thesis, ICA, University of Stuttgart, 1992. [41 J.-A. Desideri, N. Glinsky, E. Hettena, Hypersonic Reactive Flow Computations, Computers and Fluids, Vo1.18, No.2 (1990) 151-182. 151 C. Park, On Convergence of Chemically Reacting Flows, AIAA Paper 85-0247, AIAA 23rd Aerospace Science Meeting, Reno, Nevada (1985). 161 J. V. Rakich, H. E. Bailey, C. Park, Computation of Nonequilibrium, Supersonic Three-Dimensional Inviscid Flow over Blunt-Nosed Bodies, AIAA Journal, Vol. 21 (1983) 834-841. 171 N. Botta, M. Pabdolfi, M. Germano, Nonequilibrium Reacting Hypersonic Flow About Blunt Bodies: Numerical Prediction, AIAA-88-0514 (1988). k31J. Urban, Berechnung von chemischen Reactionen in reagierenden Hyperschallstromungen., Diploma Thesis, ICA, University of Stuttgart, 1989. PI H. Rieger, A. Jameson, Solution of Steady Three-Dimensional Compressible Euler and Navier-Stokes Equations by an Implicit LU-Scheme, AIAA-88-0619 (1988). WI C. D. Scott, Wall Catalytic Recombination and Boundary Conditions in Nonequilibrium Hypersonic Flows, 3rd Joint Europe/US Short Course in Hypersonics, Aachen (1990).

PI

J. Argyris et al., Physical and computationill

aspects of reacting Bows

35

J. Warnatz, Reacting Viscous Flow and Gas-Surface Interaction Modelling, Proceedings of the INHIA GAMNI / SMAI Workshop on Hypersonic Flows for Reentry Problems, Part II, Antibes (1991) 150-163. [12] J. D. Anderson, Hpersonic and High Temperature Gas Dynamics (McGraw-Hill Series in Aeronautical and Aerospace Engineering, 1989) 623-626. [13] European Hypersonic Data Base: Problems For Analysis, 1992.

[ll]