**A Magneto-Fluid-Dynamic Model and Computational Solving Methodologies for Aerospace Applications**

Francesco Battista1, Tommaso Misuri2 and Mariano Andrenucci2 *1CIRA SCpA Italian Aerospace Research Centre 2ALTA SpA Italy* 

#### **1. Introduction**

Computational plasma physics is concerned primarily with the study of the evolution of plasma by means of computer simulation. The main task of this computational branch is to develop methods able to obtain a better understanding of plasma physics. Therefore a close contact to theoretical plasma physics and numerical methods is necessary. Ideally computational plasma physics acts as a pathfinder to guide scientific and technical development and to connect experiment and theory. To build a valid computer simulation program means to devise a model which is sufficiently detailed to reproduce faithfully the most important physical effects, with a computational effort sustainable by modern computers in reasonable time (Dandy, 1993).

Computational models have played an important role in the development of plasma physics since the beginning of the computer age. Advances in our understanding of many plasma phenomena like magnetohydrodynamic instabilities, micro-instabilities, transport, wave propagation, etc. have gone hand-in-hand with the increased computational power available to researchers. Several trends are evident in how computer modelling is carried out: the models are becoming increasingly complex, for example, by coupling separate computer codes together. This allows for more realistic modelling of the plasma. Presently several efforts are carried out in different countries to develop plasma numerical tools for several applications such as fusion, electric propulsion, active control over hypersonic vehicles: these efforts lead to a growing experience in CMFD field (see Park et al. 1999, Kenneth et al 1998, Taku and Atsushi 2004, Cristofolini et al 2007, Miura and Groth 2007, MacCormack 2007, Yalim 2001, Giordano and D'Ambrosio 2004, Battista 2009).

The chapter presented was carried out in the context of a research activity motivated by renewed interest in investigating the influence that electromagnetic fields can exert on the thermal and pressure loads imposed on a body invested by a high energetic flow. In this regard, spacecraft thermal protection and the opportunity to use active control surfaces during planetary (re)entry represent the driving engineering applications. The contents of the study should be considered, to a certain extent, a systematic re-examination of past work complemented with somewhat innovative ideas.

So, in this chapter, methodologies for plasma modelling have been developed and then implemented and tested into a numerical code EMC3NS, developed in the frame of this

A Magneto-Fluid-Dynamic Model and

temperature.

the i-th specie,

\_

*i k ik i k prod term M* 

2

1. the velocity is not relativistic *u c* ;

Computational Solving Methodologies for Aerospace Applications 123

with classical aero-thermo-dynamic codes (e.g. H3NS developed by CIRA) that are nearly in the same field of application. For these scopes three set of equations have been analyzed: multi-fluid equations (Wagner 1998), MHD equations (Helander 2006), and MFD equations

The first set of equations (three fluid equations), where each species (charged and not) is considered as a single fluid interacting with the other through collisions, is one of the more appropriate way to describe plasma flows. However, even if the collision process is well described, it requires a large computational cost because of the large number of equations to be solved (especially in presence of more chemical species (Giordano 2002)). The MHD set of equations describes the motion and the electrodynamics of an electrically neutral but electrically conducting fluid. These assumptions do not permit to have any information about the species present in the flowfield; moreover diffusive transport and charge separation are excluded by this model. These drawbacks are critical in the exclusion of this approach for our purposes. The MFD equations for a single multispecies fluid is the set of equations that more than the others is suitable for our scope since it does not loose information about species (charged and not), and furthermore the collisions process can be modelled by transport coefficients. Besides this model is not computationally heavy and the equation structure is easily adaptable to the one used in typical finite volume codes for aerothermodynamics application. The model that will be used and described further is able to consider a multi-temperature gas with vibrational temperature and electronic temperature in non equilibrium with the translational

for a single fluid composed by more than one species (Giordano 2002).

The assumptions used to write down the system of equations are the following:

*ii*

The equations are written in the following; mass balance equations of the components in

*i im i <sup>i</sup> V S dV dV V J ds dV <sup>t</sup>*

 

*e m m t*

 ,

*i me T*

*p T M*

*<sup>D</sup> L* , *<sup>i</sup> L* .

(2.1)

2

/ *ce*

*ne v L <sup>t</sup> <sup>m</sup>* 

;

2 are the production rate of

0

*e*

2. the phenomena under consideration are slow enough that

where *mi J* 1are the components mass diffusive fluxes and the *<sup>i</sup>*

<sup>1</sup> <sup>2</sup>

1 1 *<sup>i</sup> n*

 

*m j ij j i i j <sup>M</sup> J diffusive flux i component MD c D E V B D T*

3. the plasma is sufficiently collisional *<sup>i</sup>*

chemical non equilibrium conditions yield,

1

*r*

4. the Debye length and the ion gyro radius are small

work. This one has been used and tested primarily to support the design of experiments in problems of active flow control over bodies in an electromagnetic field in order to modify shock waves position, and consequently the aerothermal environment.

These kind of experiments are not so common in the scientific community especially using air as working gas, thus the lack of experimental data and lesson learned in this field make numerical computation play an important role. So the use of valid computer simulation is crucial for a correct design of a flow control experiment set up diagnostics and to understand main physical phenomena related to this kind of problems.

The use of magnetic fields to control external flows is not new (for instance see Bush (1958), Resler & Sears (1958)); recent computational and experimental technologies have moved this approach from mere possibility to real application, highlighting the considerable advantage that may be gained in re-entry operation.

The physical ideas are incorporated into the so-called magneto-fluid-dynamic interaction concept: global body forces can be applied to a weakly ionized plasma using electromagnetic devices embedded in the vehicle. Therefore, there is a growing interest in using weakly ionized gases (plasmas) and electric and magnetic fields in high-speed aerodynamics. Wave and viscous drag reduction, thrust vectoring, reduction of heat fluxes, sonic boom mitigation, boundary-layer and turbulent transition control, flow turning and compression, on-board power generation, and scramjet inlet control are among plasma and MHD technologies that can potentially enhance performance and significantly change the design of supersonic and hypersonic vehicles and thrusters.

Meanwhile, despite many studies devoted to these new technologies, a number of fundamental issues have not been adequately addressed. Any plasma created in gas flow and interacting with electric and magnetic fields would result in gas heating. This heating can certainly have an effect on the flow and, in some cases, can be used advantageously. However, a more challenging issue is whether significant non-thermal effects of plasma interaction with electric and magnetic fields can be used for high-speed flow control. In conventional MHD of highly conducting fluid, electric and magnetic effects give rise to ponderomotive force terms, which can be interpreted as gradients of electric and magnetic field pressures. These ponderomotive forces are successfully used for plasma containment in fusion devices and also play an important role in astrophysics. One might hope that these forces can also be used for control of high-speed flow of ionized air. However, the great importance of ponderomotive forces in fusion and astrophysical plasmas is due to the fact that those plasmas are fully, or almost fully, ionized and, therefore, are highly conductive. In contrast, high speed air encountered in aerodynamics is not naturally ionized, even in boundary layers and behind shocks if the flight Mach number is below about 8, due to the low static temperature. So in this case studies are necessary to set up technologies and methodologies to control the flow in such weak ionized regime.

#### **2. Governing equations**

In the present section the set of equations to be solved is exposed; as previously described, the choice of a model, and thus the acceptance of all the underlying hypotheses, determines the nature of the results that numerical computations provide. The analysis has been carried out in order to provide a set of equations suitable for the problematic exposed in the introduction, with a limited computational cost and that have a structure easily adaptable

work. This one has been used and tested primarily to support the design of experiments in problems of active flow control over bodies in an electromagnetic field in order to modify

These kind of experiments are not so common in the scientific community especially using air as working gas, thus the lack of experimental data and lesson learned in this field make numerical computation play an important role. So the use of valid computer simulation is crucial for a correct design of a flow control experiment set up diagnostics and to

The use of magnetic fields to control external flows is not new (for instance see Bush (1958), Resler & Sears (1958)); recent computational and experimental technologies have moved this approach from mere possibility to real application, highlighting the considerable

The physical ideas are incorporated into the so-called magneto-fluid-dynamic interaction concept: global body forces can be applied to a weakly ionized plasma using electromagnetic devices embedded in the vehicle. Therefore, there is a growing interest in using weakly ionized gases (plasmas) and electric and magnetic fields in high-speed aerodynamics. Wave and viscous drag reduction, thrust vectoring, reduction of heat fluxes, sonic boom mitigation, boundary-layer and turbulent transition control, flow turning and compression, on-board power generation, and scramjet inlet control are among plasma and MHD technologies that can potentially enhance performance and significantly change the

Meanwhile, despite many studies devoted to these new technologies, a number of fundamental issues have not been adequately addressed. Any plasma created in gas flow and interacting with electric and magnetic fields would result in gas heating. This heating can certainly have an effect on the flow and, in some cases, can be used advantageously. However, a more challenging issue is whether significant non-thermal effects of plasma interaction with electric and magnetic fields can be used for high-speed flow control. In conventional MHD of highly conducting fluid, electric and magnetic effects give rise to ponderomotive force terms, which can be interpreted as gradients of electric and magnetic field pressures. These ponderomotive forces are successfully used for plasma containment in fusion devices and also play an important role in astrophysics. One might hope that these forces can also be used for control of high-speed flow of ionized air. However, the great importance of ponderomotive forces in fusion and astrophysical plasmas is due to the fact that those plasmas are fully, or almost fully, ionized and, therefore, are highly conductive. In contrast, high speed air encountered in aerodynamics is not naturally ionized, even in boundary layers and behind shocks if the flight Mach number is below about 8, due to the low static temperature. So in this case studies are necessary to set up technologies and

In the present section the set of equations to be solved is exposed; as previously described, the choice of a model, and thus the acceptance of all the underlying hypotheses, determines the nature of the results that numerical computations provide. The analysis has been carried out in order to provide a set of equations suitable for the problematic exposed in the introduction, with a limited computational cost and that have a structure easily adaptable

shock waves position, and consequently the aerothermal environment.

understand main physical phenomena related to this kind of problems.

advantage that may be gained in re-entry operation.

design of supersonic and hypersonic vehicles and thrusters.

methodologies to control the flow in such weak ionized regime.

**2. Governing equations** 

with classical aero-thermo-dynamic codes (e.g. H3NS developed by CIRA) that are nearly in the same field of application. For these scopes three set of equations have been analyzed: multi-fluid equations (Wagner 1998), MHD equations (Helander 2006), and MFD equations for a single fluid composed by more than one species (Giordano 2002).

The first set of equations (three fluid equations), where each species (charged and not) is considered as a single fluid interacting with the other through collisions, is one of the more appropriate way to describe plasma flows. However, even if the collision process is well described, it requires a large computational cost because of the large number of equations to be solved (especially in presence of more chemical species (Giordano 2002)). The MHD set of equations describes the motion and the electrodynamics of an electrically neutral but electrically conducting fluid. These assumptions do not permit to have any information about the species present in the flowfield; moreover diffusive transport and charge separation are excluded by this model. These drawbacks are critical in the exclusion of this approach for our purposes. The MFD equations for a single multispecies fluid is the set of equations that more than the others is suitable for our scope since it does not loose information about species (charged and not), and furthermore the collisions process can be modelled by transport coefficients. Besides this model is not computationally heavy and the equation structure is easily adaptable to the one used in typical finite volume codes for aerothermodynamics application. The model that will be used and described further is able to consider a multi-temperature gas with vibrational temperature and electronic temperature in non equilibrium with the translational temperature.

The assumptions used to write down the system of equations are the following:


The equations are written in the following; mass balance equations of the components in chemical non equilibrium conditions yield,

$$\frac{\partial}{\partial t} \int\_{V} \rho\_i dV + \int\_{S} \left( \rho\_i \underline{V} + \underline{J\_{m\_i}} \right) \cdot \underline{ds} = \int\_{dV} \Omega\_i dV \tag{2.1}$$

where *mi J* 1are the components mass diffusive fluxes and the *<sup>i</sup>* 2 are the production rate of the i-th specie,

 <sup>1</sup> <sup>2</sup> 1 1 *<sup>i</sup> n i me T m j ij j i i j <sup>M</sup> J diffusive flux i component MD c D E V B D T p T M* 2 1 \_ *r i k ik i k prod term M* 

A Magneto-Fluid-Dynamic Model and

equilibrium constants *Keq,r*, i.e.:

ionization are:

where *'ir* and

Computational Solving Methodologies for Aerospace Applications 125

chemical species involved in the reactions, the production rate *Ωi* of species *i* can be written as

where *kfr* and *kbr* are the forward and backward rate constants of the reaction r, which are

*fr*

*br*

The backward rate constants are related to the forward rate constants through the

*fr*

*br k*

*r*

In some cases the equilibrium coefficients can be given as a function of the temperature, and

For argon the reaction coefficients that have been considered in order to account for

n. Reaction A [mol-cm-s-K] β Ea [J/mol] 1 AR AR+ + e- 1.52e18 0.50 1520000 2 AR + e AR+ + 2e- 2.50e34 -3.8 1510793

> 2 2 12

*k T k A T e Pd N e*

this formula is consequent from the expression of binary collision rate where m12 is the reduced mass, d12 is the mean diameter, kB is the Boltzmann constant and P is the steric

For air ionization one of the most important reactions governing the distribution of free

*N O NO e* In this work, in order to model air ionization several 7-species models and 11 species models have been used (Gupta et al. 1989, Park 1993, Kang et al.1973) and compared also with a

<sup>8</sup> *E E a a T T B A*

 

12

*m*

where the reaction rate coefficients of the first reaction have been derived from:

factor. The second reaction rate coefficient has been taken from (Gokcet 2004).

electrons present in high temperature air plasmas is (Dunn & Lordi 1969):

 *eq*

*K*

*fr fr k AT e* 

*br br k AT e* 

' '' '' '

1 1 1

*r s s*

*N N N*

( )

 

assumed to follow the Arrhenius temperature law:

energies *Ei* depend upon the adopted kinetic scheme.

*i i ir ir fr br r i i M k qk q*

where the pre-exponential factors *Ai*, the temperature exponents

the backward reaction rate can be calculated directly from the Eqn. (2.9).

*''ir* are the stoichiometric coefficients for the r-th reaction and

*ir ir*

 

*fr*

*br*

*E RT*

*E RT*

*i N* 1,..., 1 *<sup>s</sup>* (2.7)

(2.8)

(2.9)

*<sup>k</sup>* (2.10)

*<sup>r</sup>* , and the activation

*<sup>i</sup>* are the

Pointing out the charge density as 1 *<sup>n</sup> <sup>i</sup> C A iS i i eN M* (where *iS* is equal to (+1) for electrons, to (-1) for ions and to (0) for neutral components) and summing up for all the charged species, then the electric charge balance equation is obtained

$$\frac{\partial}{\partial t} \int\_{V} \rho\_c dV + \int\_{S} (\underline{\mathbf{J}\_q} + \rho\_c \underline{\mathbf{V}}) \cdot \underline{dS} = \mathbf{0} \tag{2.2}$$

where the current is the sum of convection and conduction current *J VJ c Q* 3; the momentum equation is

$$\frac{\partial}{\partial t} \int\_{V} \rho \underline{V}dV + \int\_{S} \left(\rho \underline{V} \underline{V} + p \underline{I} - \underline{\tau}\right) \cdot \underline{dS} = \int\_{V} \left(\rho\_{c} \underline{E} + \underline{I} \times \underline{B}\right)dV\tag{2.3}$$

where , <sup>2</sup> 3 *j i i j ij i j u u V x x* is the viscous stress tensor and the total energy equation reads:

$$\frac{\partial}{\partial t} \int\_{V} \rho e\_{m}dV + \iint\_{S} \left( (\rho e\_{m} + p)\underline{V} + \underline{f\_{U}} - \underline{\tau} \cdot \underline{V} \right) \cdot \underline{dS} = \int\_{V} (\underline{f} \cdot \underline{E})dV + \int\_{\mathcal{S}V} \underline{q}\_{r} \cdot \underline{dS} \tag{2.4}$$

where the total energy is 2 1 2 *<sup>B</sup> <sup>m</sup> vibr k T mv e e* , the heat flux due to radiation is

 2 2 00 0 , cos sin *<sup>r</sup> q n I ddd* .

The vibrational energy equations:

$$\frac{\partial}{\partial t} \int\_{V} \rho \mathbf{e}\_{v\bar{\jmath}}dV + \int\_{\partial V} \underline{n} \cdot \left(\rho \mathbf{e}\_{v\bar{\jmath}} \underline{V}\right) ds = \int\_{V} \rho \dot{\mathbf{e}}\_{v\bar{\jmath}}dV\tag{2.5}$$

provides the terms in the Eqn.(2.4).

The opportunity of considering electronic temperature in the vibrational energy equations is currently object of study. The production terms of the vibrational energy equation are evaluated through the classical Landau Teller non equilibrium equation (Vincenti & Krouger 1967).

#### **2.1 Chemical model**

Given a generic set of Nr chemical reactions,

$$\sum\_{i=1}^{N\_s} \nu\_{ir}^{\cdot} \mathbb{X}\_i \Leftrightarrow \sum\_{i=1}^{N\_s} \nu\_{ir}^{\cdot} \mathbb{X}\_i \tag{2.6}$$

$$\stackrel{\circ}{\quad} \underline{J\_Q} = \sum\_{j=1}^n \overline{\mathcal{A}}\_{e^j}^p \cdot \nabla c\_j + \overline{\mathcal{A}\_e} \cdot \left( \underline{E} + \underline{V} \times \underline{B} \right)$$

*eN*

electrons, to (-1) for ions and to (0) for neutral components) and summing up for all the

*dV J V dS <sup>t</sup>*

where the current is the sum of convection and conduction current *J VJ*

*V S*

*V S V*

charged species, then the electric charge balance equation is obtained

*V*

 

 

*t* 

3

.

Given a generic set of Nr chemical reactions,

 

*j i i j ij i j u u*

*x x*

 

where , <sup>2</sup>

*t*

where the total energy is

, cos sin *<sup>r</sup> q n I ddd*

The vibrational energy equations:

provides the terms in the Eqn.(2.4).

 

2 2

 

Krouger 1967).

1 *<sup>n</sup> <sup>p</sup> Q ej j e*

*j*

**2.1 Chemical model** 

<sup>3</sup>

*J c EVB*

 

00 0

1 *<sup>n</sup> <sup>i</sup> C A iS i i*

( )0 *c qc*

 

*<sup>c</sup>*

 

(2.3)

*VdV VV pI dS E J B dV <sup>t</sup>*

*m mU <sup>r</sup> V S V V e dV e p V J V dS J E dV q dS*

2

*vj vj vj VV V*

The opportunity of considering electronic temperature in the vibrational energy equations is currently object of study. The production terms of the vibrational energy equation are evaluated through the classical Landau Teller non equilibrium equation (Vincenti &

> ' '' 1 1

*ir i ir i*

 

*N N s s*

*i i* 

*e dV n e V ds e dV*

1 2 *<sup>B</sup> <sup>m</sup> vibr k T mv e e*

 

 (2.4)

 

*M* 

(2.2)

is the viscous stress tensor and the total energy

, the heat flux due to radiation is

(2.5)

(2.6)

*iS* is equal to (+1) for

*c Q* 3; the

(where

Pointing out the charge density as

momentum equation is

equation reads:

where *'ir* and *''ir* are the stoichiometric coefficients for the r-th reaction and *<sup>i</sup>* are the chemical species involved in the reactions, the production rate *Ωi* of species *i* can be written as

$$\boldsymbol{\Omega}\_{i} = \rho \boldsymbol{M}\_{i} \sum\_{r=1}^{N\_{r}} (\boldsymbol{\nu}\_{ir}^{\ast} - \boldsymbol{\nu}\_{ir}^{\cdot}) \left( \boldsymbol{k}\_{fr} \prod\_{i=1}^{N\_{s}} \boldsymbol{q}^{\nu\_{ir}^{\cdot}} - \boldsymbol{k}\_{br} \prod\_{i=1}^{N\_{s}} \boldsymbol{q}^{\nu\_{ir}^{\cdot}} \right) \qquad \left( \boldsymbol{i} = 1, \dots, N\_{s} - 1 \right) \tag{2.7}$$

where *kfr* and *kbr* are the forward and backward rate constants of the reaction r, which are assumed to follow the Arrhenius temperature law:

$$k\_{fr} = A\_{fr} T^{\beta\_{fr}} \cdot e^{\frac{-E\_{fr}}{RT}} \tag{2.8}$$

$$k\_{br} = A\_{br}T^{\beta\_{br}} \cdot e^{\frac{-E\_{br}}{RT}} \tag{2.9}$$

where the pre-exponential factors *Ai*, the temperature exponents *<sup>r</sup>* , and the activation energies *Ei* depend upon the adopted kinetic scheme.

The backward rate constants are related to the forward rate constants through the equilibrium constants *Keq,r*, i.e.:

$$K\_{eq\_r} = \frac{k\_{fr}}{k\_{br}}\tag{2.10}$$

In some cases the equilibrium coefficients can be given as a function of the temperature, and the backward reaction rate can be calculated directly from the Eqn. (2.9).

For argon the reaction coefficients that have been considered in order to account for ionization are:


where the reaction rate coefficients of the first reaction have been derived from:

$$k = A \ T^{\beta} \cdot e^{\frac{-E\_s}{\Re T}} = P d^2 \, \_{12}N \, \_{A}^{2} \sqrt{\frac{8\pi k\_B T}{m\_{12}}} e^{\frac{-E\_s}{\Re T}}$$

this formula is consequent from the expression of binary collision rate where m12 is the reduced mass, d12 is the mean diameter, kB is the Boltzmann constant and P is the steric factor. The second reaction rate coefficient has been taken from (Gokcet 2004).

For air ionization one of the most important reactions governing the distribution of free electrons present in high temperature air plasmas is (Dunn & Lordi 1969):

$$N + O \leftrightarrow NO^{+} + e^{-}$$

In this work, in order to model air ionization several 7-species models and 11 species models have been used (Gupta et al. 1989, Park 1993, Kang et al.1973) and compared also with a

A Magneto-Fluid-Dynamic Model and

Hence, for a mixture of thermally perfect gases:

**2.2 Thermodynamic model** 

heats of the single species i-th:

McBride 1971) polynomial fits:

temperature.

Computational Solving Methodologies for Aerospace Applications 127

By definition a thermally perfect gas is characterized by specific heats that depend only on

*p* 

*i i i*

Fig. 1. Enthalpy fit extension for temperatures ranging from 6000K and 15000K.

*p c*

*R*

where *R* and *cp* are, respectively, the gas constant and the specific heat at constant pressure of the mixture. They are obtained as weighted average of the gas constants and specific

*R YR p i pi*

The single species specific heats and enthalpies are computed by using the (Gordon and

12 3 4 5

<sup>1</sup> 23 4 5 *ha a a aa a TT T T RT T*

*a aT aT aT aT*

*i*

2 4 3 56 234

234

(2.14)

(2.15)

*dh c T dT <sup>p</sup>* (2.11)

*RT* (2.12)

*c T Yc T* (2.13)

model generated using (Park 1993) and (Kang et al.1973) derived in order to better fit experimental data w.r.t. the original 11 species models. The scheme adopted is reported hereinafter specifying for each reaction the respective source.


Table 1. Kinetic scheme obtained using (Park 1993) and (Kang et al.1973)

#### **2.2 Thermodynamic model**

126 Fluid Dynamics, Computational Modeling and Applications

model generated using (Park 1993) and (Kang et al.1973) derived in order to better fit experimental data w.r.t. the original 11 species models. The scheme adopted is reported

n. Reaction A [mol-cm-s-K] Β Ea [J/mol] Ref

1 O2 2 O 1.0e22 -1.5 494706.8 P 2 N2 2 N 3.0e22 -1.6 941190.08 P 3 NO N + O 1.1e17 0.0 627737.2 P 4 NO + O N + O2 8.4e12 0.0 161715.08 P 5 N2 + O N + NO 6.4e17 1.0 319257.96 P 6 N + O NO+ + e- 1.6e06 1.5 265892.8 K 7 N + O NO+ + e- 6.7e21 -1.5 0 K 8 O + O O2+ + e- 1.6e17 -0.98 671803.52 K 9 O2+ + e- O + O 8.0e21 -1.5 0 K 10 N + N N2+ + e- 1.4e13 0 671803.52 K 11 N2+ + e- N + N 8.0e21 -1.5 0 K 10 O2 + N2 NO + NO+ + e- 1.38e20 -1.84 1172330.4 K 12 NO + NO+ + e- O2 + N2 8.0e21 -2.5 0 K 13 NO + N2 NO+ + N2 + e- 2.2e15 -0.35 897955.2 K 14 NO+ + N2 + e- NO + N2 2.2e26 -2.5 0.0 K 15 NO + O2 NO+ + O2 + e- 8.8e15 -0.35 897955.2 K 16 NO+ + O2 + e- NO + O2 8.8e26 -2.5 0 K

+ e- 3.9e33 -3.78 1317832.4 P

+ e- 8.8e26 -3.82 1401807.84 P

<sup>+</sup> 7.80e11 0.5 0.0 P

+ 2.02e18 0.81 108087.2 P

<sup>+</sup> O+ + O2 2.92e18 -1.11 232803.2 P

<sup>+</sup> N2 + N+ 7.80e11 0.5 0.0 P

<sup>+</sup> O2 + NO+ 1.50e13 0.0 0.0 P

23 O + NO+ NO + O+ 3.63e15 -0.6 422371.52 P 24 NO + O+ O + NO+ 1.50e13 0.0 0.0 P 25 N + NO+ NO + N+ 1.00e19 -0.93 507178.4 P 26 NO + N+ N + NO+ 7.80e11 0.0 0.0 P 27 O2 + NO+ NO + O2+ 1.80e15 -0.6 274375.2 P

29 O + NO+ N++ O2 1.34e13 0.31 642453.69 P 30 N++ O2 O + NO+ 1.00e14 0.0 0.0 P

Table 1. Kinetic scheme obtained using (Park 1993) and (Kang et al.1973)

hereinafter specifying for each reaction the respective source.

17 O + e- O+ + e-

18 N + e- N+ + e-

20 O+ + O2 O + O2

21 N2 + N+ N + N2

19 O + O2

22 N + N2

28 NO + O2

By definition a thermally perfect gas is characterized by specific heats that depend only on temperature.

Hence, for a mixture of thermally perfect gases:

$$d\mathbf{l} = c\_p \begin{pmatrix} T \\ \end{pmatrix} dT \tag{2.11}$$

$$p = \rho RT\tag{2.12}$$

where *R* and *cp* are, respectively, the gas constant and the specific heat at constant pressure of the mixture. They are obtained as weighted average of the gas constants and specific heats of the single species i-th:

$$R = \sum\_{i} Y\_i R\_i \ c\_p(T) = \sum\_{i} Y\_i c\_{pi}(T) \tag{2.13}$$

The single species specific heats and enthalpies are computed by using the (Gordon and McBride 1971) polynomial fits:

$$\frac{c\_p}{R} = a\_1 + a\_2 T + a\_3 T^2 + a\_4 T^3 + a\_5 T^4 \tag{2.14}$$

$$\frac{h}{RT} = a\_1 + \frac{a\_2}{2}T + \frac{a\_3}{3}T^2 + \frac{a\_4}{4}T^3 + \frac{a\_5}{5}T^4 + \frac{a\_6}{7} \tag{2.15}$$

Fig. 1. Enthalpy fit extension for temperatures ranging from 6000K and 15000K.

A Magneto-Fluid-Dynamic Model and

Considering collision frequency as

could be found in (Baum 1965).

full MFD system: equations are the following:

they can be rewritten in conservative form,

spatial direction, neglecting time derivatives as:

**2.5 Maxwell equations** 

**2.6 Collision modeling** 

where *He* 

**2.4 Electrical conductivity** 

Computational Solving Methodologies for Aerospace Applications 129

In order to evaluate electrical conductivity its classical expression has been written down:

<sup>0</sup> *<sup>e</sup>*

where *υeH* is the collision frequency between electrons and heavy particles measured in Hz.

<sup>8</sup> 1081794.2

velocity (under Maxwellian distribution hypothesis) the following expression yield:

<sup>12</sup> 104.535909

Maxwell equations have to be solved together with the fluid equations in order to solve the

 *V V*

0 0

*dVE*

*nn* is the degree of ionization. For what concerns *QeH* an expression for air

By substituting the constant values the expression becomes:

<sup>2</sup>

*t*

*t*

*V V*

*V V*

*V V <sup>V</sup>*

*eH*

2

*e eH n e m*

> *e n*

*eH e Q T* 2/1

*dVB*

*t*

0

(2.19)

0

(2.20)

*t*

 

0

*<sup>E</sup> c B dV J dV*

*BdV n E dS*

<sup>2</sup>

*E J dV c n B dS dV*

In order to correctly describe collision processes that in this kind of equations determine the transport properties of the fluid, a momentum equation could be written down for each

.

*eH H eH e nQ u* and substituting the electron mean

(2.16)

(2.17)

(2.18)

These fits are written for gases with the maximum temperature of 6000K, so it has been necessary to extend the enthalpy fits to values of temperature typical of ionization problem (up to 15000K). Here is briefly reported the fitting procedure for argon: from literature (Drellishak et al. 1963) thermodynamic properties are derived for high temperatures then the coefficients are found using a least squares fitting procedure. The results for Argon enthalpy are shown in Fig. 1. Thermodynamic properties for high temperature air, used to find fit coefficients, have been found in (Hans & Heims 1968).

#### **2.3 Transport model**

Another critical point in the simulation of high enthalpy flows is the determination of the transport coefficients. In fact, the widely used Sutherland law is suitable only at low temperatures, while more complex models must be used when temperature exceeds 1000 K. An interesting way to compute transport coefficients is based on the calculation of the collision integral starting from intermolecular potentials knowledge (Hirshfield et al. 1954) and it leads to the following expressions for the diffusivities between the i-th and the j-th species and also for conductivity and viscosity with different sets of constants given in Gupta et al. 1958, Hans & Heims 1968. In this way computational efficiency is maximized since transport coefficients are computed only once at the start of each calculation.

$$\ln D\_y = \sum\_{n=1}^{N} d\_{n,y} \left( \ln T \right)^{n-1} \qquad \ln \mathcal{A}\_\iota = \sum\_{n=1}^{N} b\_{n,\iota} \left( \ln T \right)^{n-1} \qquad \ln \mu\_\iota = \sum\_{n=1}^{N} a\_{n,\iota} \left( \ln T \right)^{n-1} \qquad \Box$$

Once derived the single species properties, total conductivity and viscosity are calculated using Wilke (Hirshfielder 1954) formulas. The mixture diffusion coefficient for species i-th is obtained as reported in Bird (1954).

Polynomial coefficients for ionized species have been found fitting data from (Capitelli et al. 2000); the results of the fitting procedure for NO+, N+ and O+ are shown for viscosity in Fig. 2.

Fig. 2. Viscosity in function of temperature for charged species NO+, N+ and O+

#### **2.4 Electrical conductivity**

128 Fluid Dynamics, Computational Modeling and Applications

These fits are written for gases with the maximum temperature of 6000K, so it has been necessary to extend the enthalpy fits to values of temperature typical of ionization problem (up to 15000K). Here is briefly reported the fitting procedure for argon: from literature (Drellishak et al. 1963) thermodynamic properties are derived for high temperatures then the coefficients are found using a least squares fitting procedure. The results for Argon enthalpy are shown in Fig. 1. Thermodynamic properties for high temperature air, used to

Another critical point in the simulation of high enthalpy flows is the determination of the transport coefficients. In fact, the widely used Sutherland law is suitable only at low temperatures, while more complex models must be used when temperature exceeds 1000 K. An interesting way to compute transport coefficients is based on the calculation of the collision integral starting from intermolecular potentials knowledge (Hirshfield et al. 1954) and it leads to the following expressions for the diffusivities between the i-th and the j-th species and also for conductivity and viscosity with different sets of constants given in Gupta et al. 1958, Hans & Heims 1968. In this way computational efficiency is maximized

> 1 , ln ln *<sup>n</sup> <sup>N</sup>*

Once derived the single species properties, total conductivity and viscosity are calculated using Wilke (Hirshfielder 1954) formulas. The mixture diffusion coefficient for species i-th is

Polynomial coefficients for ionized species have been found fitting data from (Capitelli et al. 2000); the results of the fitting procedure for NO+, N+ and O+ are shown for viscosity

*Tb*

*n i n i* <sup>1</sup>

<sup>1</sup>

.

1 , ln ln *<sup>n</sup> <sup>N</sup>*

*Ta*

*n i n i*

since transport coefficients are computed only once at the start of each calculation.

Fig. 2. Viscosity in function of temperature for charged species NO+, N+ and O+

find fit coefficients, have been found in (Hans & Heims 1968).

<sup>1</sup>

1 , ln ln 

*n ij <sup>n</sup> ij TdD*

obtained as reported in Bird (1954).

in Fig. 2.

*<sup>n</sup> <sup>N</sup>*

**2.3 Transport model** 

In order to evaluate electrical conductivity its classical expression has been written down:

$$
\sigma\_0 = \frac{n\_e e^2}{m\_e \nu\_{eH}}
$$

where *υeH* is the collision frequency between electrons and heavy particles measured in Hz. By substituting the constant values the expression becomes:

$$
\sigma = 2.81794 \left| 10^{-8} \frac{n\_e}{\nu\_{eH}} \right|
$$

Considering collision frequency as *eH H eH e nQ u* and substituting the electron mean velocity (under Maxwellian distribution hypothesis) the following expression yield:

$$
\sigma = 4.535909 \, 10^{-12} \, \alpha \frac{T\_\ast^{-1/2}}{\mathcal{Q}\_{\ast H}} \tag{2.16}
$$

where *He nn* is the degree of ionization. For what concerns *QeH* an expression for air could be found in (Baum 1965).

#### **2.5 Maxwell equations**

Maxwell equations have to be solved together with the fluid equations in order to solve the full MFD system: equations are the following:

$$\int\_{V} (\nabla \times \underline{E})dV = -\frac{\partial}{\partial t} \int\_{V} \underline{B}dV\tag{2.17}$$

$$\int\_{V} \left( \varepsilon\_{0} c^{2} \nabla \times \underline{B} \right) dV = \int\_{V} \left( \underline{J} + \varepsilon\_{0} \frac{\partial \underline{E}}{\partial t} \right) dV \tag{2.18}$$

they can be rewritten in conservative form,

$$\frac{\partial}{\partial t} \int\_{V} \underline{B}dV + \int\_{\mathcal{V}} \left(\underline{n} \times \underline{E}\right)dS = 0\tag{2.19}$$

$$
\int\_{V} \left(\frac{\partial \underline{E}}{\partial t}\right) dV - \int\_{\mathcal{V}} \left(c^{2} \underline{\eta} \times \underline{B}\right) dS + \int\_{\mathcal{E}\_{0}} \frac{J}{\mathcal{E}\_{0}} dV = 0 \tag{2.20}
$$

#### **2.6 Collision modeling**

In order to correctly describe collision processes that in this kind of equations determine the transport properties of the fluid, a momentum equation could be written down for each spatial direction, neglecting time derivatives as:

A Magneto-Fluid-Dynamic Model and

1

*m x y z x y z*

*u v w e B B B E*

.. ..

*E*

 

*<sup>W</sup> <sup>E</sup>*

1 ,1

*n v*

*e*

*e*

,

*v Nv*

...

Computational Solving Methodologies for Aerospace Applications 131

 

*m z*

 2

*m y x*

*w uw vw w p p e w E E*

 

2

0

1

*w*

.. ..

0

*c B*

*x*

1 1

*w e w*

*vNv*

*e w*

*n v*

...

 

*qx z c z*

 

*J B uB*

0 0 0

1

.. ..

*n v*

...

 

*Nv*

0 0 0

*qx c qy c qz c*

 

*J u J v J w*

 

 

 

*c z qx y c y qy x c x qxx c x qy y c y qzz c z*

*E J B uB J B vB J E uE J E vE J E wE*

 

1 1

, <sup>2</sup>

*<sup>y</sup> inv z*

*c B <sup>F</sup>*

2

*E*

*x*

0

0

2

*c B v*

1

*x*

.. ..

> 1 1

*v e v*

*vNv*

*e u*

*n v*

...

 

2

*v uv v p uw p e v E*

0

*E E*

*z y*

*m*

2

*c B c B u*

*z y* ,

*<sup>z</sup> inv y*

*c B <sup>F</sup>*

0

1

.. ..

> 1 1

*u e u*

*vNv*

*e u*

0 0

*c y qz x c y cy zx xz*

*E JB JB E J B wB*

*c x qy z c z qz y c z cx yz zy*

*E JB JB E J B vB J B wB*

*n v*

...

, 2

*inv x*

*x y z*

*J J J* 0 0 0

*cz xy yx xx yy zz*

*E JB JB JE JE JE*

 

> 1 1

.. ..

*n v*

*Nv*

...

 

*F*

2

 

*u u p uv uw p e u*

$$\begin{cases} u\_x = \mu E\_x - \frac{D}{n} \frac{\partial n}{\partial x} + \left( u\_y \phi\_{cz} - u\_z \phi\_{cy} \right) \Big|\_V \\ u\_y = \mu E\_y - \frac{D}{n} \frac{\partial n}{\partial y} + \left( u\_z \phi\_{cx} - u\_x \phi\_{cz} \right) \Big|\_V \\ u\_z = \mu E\_z - \frac{D}{n} \frac{\partial n}{\partial z} + \left( u\_x \phi\_{cy} - u\_y \phi\_{cc} \right) \Big|\_V \end{cases} \tag{2.21}$$

In the previous equations the cyclotronic vector *T c xyz ie <sup>q</sup> BBB m* has been considered that divided by gives the Hall parameter vector that represent the Hall parameter along all the space dimension: *T x yz* .

Decoupling the equations multiplying by *ne*, and reminding from previous section that 2 0 *ne m* , with a matrix pivoting technique it can be obtained the relation between electric field and current, i.e.

$$
\begin{pmatrix} j\_{qr} \\ j\_{qr} \\ j\_{qr} \end{pmatrix} = \sigma\_0 \begin{bmatrix} \frac{1+\beta^2\_{-z}}{1+\beta^2\_{-z}+\beta^2\_{-y}+\beta^2\_{-z}} & \frac{\beta\_{\gamma}\beta\_{z}+\beta\_{z}^2}{1+\beta^2\_{-z}+\beta^2\_{-y}+\beta^2\_{z}} & \frac{\beta\_{\gamma}\beta\_{z}-\beta\_{\gamma}}{1+\beta^2\_{-z}+\beta^2\_{-y}+\beta^2\_{-z}}\\ \frac{\beta\_{\gamma}\beta\_{x}-\beta\_{z}}{1+\beta^2\_{-z}+\beta^2\_{-y}+\beta^2\_{z}} & \frac{1+\beta^2\_{-z}}{1+\beta^2\_{-z}+\beta^2\_{-y}+\beta^2\_{z}} & \frac{\beta\_{\gamma}\beta\_{z}+\beta\_{z}}{1+\beta^2\_{-z}+\beta^2\_{-y}+\beta^2\_{-z}}\\ \frac{\beta\_{z}\beta\_{x}+\beta\_{y}}{1+\beta^2\_{-z}+\beta^2\_{-y}+\beta^2\_{-z}} & \frac{\beta\_{z}\beta\_{y}-\beta\_{x}}{1+\beta^2\_{-z}+\beta^2\_{-y}+\beta^2\_{-z}} & \frac{1+\beta^2\_{-z}}{1+\beta^2\_{-z}+\beta^2\_{-y}+\beta^2\_{-z}} \end{pmatrix} E\_{\gamma} - \frac{\text{Dm\nu}}{\text{n}^2\text{e}} \frac{\partial\mathbf{n}}{\partial\mathbf{y}} \tag{2.22}
$$

From Eqn. (2.22) the tensor *<sup>e</sup>* is directly derived.

#### **3. Numerical discretization of the full system and solving methodologies**

In the previous parts the fundamental equations governing an MFD plasma flow have been described. They constitute a system of 11 scalar equations plus one scalar equation for each species considered and one scalar equation for each vibrating species considered. In the following the numerical strategy and methods adopted to solve the systems will be exposed All the equations briefly reported in the previous section have the following conservative form:

$$\frac{\partial}{\partial t} \int\_{V} (\underline{\mathbf{V} \mathbf{\hat{n}}}) dV + \int\_{\partial V} (\underline{\mathbf{F} \cdot \mathbf{n}}) dS = \int\_{V} \underline{\mathbf{\Omega}} dV \tag{3.1}$$

so they can be simply discretized following a finite volume approach:

$$V\frac{d\underline{\mathcal{W}}}{dt} + \sum\_{\beta=1,\delta} (\underline{\mathcal{F}} \cdot \underline{\mathcal{u}} \cdot \Delta \mathcal{S})\_{\beta} = V\underline{\Omega} \tag{3.2}$$

The unknown vector *W* , the flux vector ( *F* ) and the source term vector ( ) used in the finite volume discretization are reported below.

*z cx x cz y y*

*D n u u u E n x D n u u u E n y D n u u u E n z*

*x x*

*z z*

In the previous equations the cyclotronic vector

 *x yz* .

*T*

2

1

 

 

divided by

0 *ne m*

form:

space dimension:

2

field and current, i.e.

*qx*

*j j j*

*qy*

*qz*

From Eqn. (2.22) the tensor

 ( )

Decoupling the equations multiplying by *ne*, and reminding from previous section that

222 222 222 2 0 222 222 222

*yx z y yz x*

*x yx z xz y zyx zyx zyx*

1

 

 

111

*<sup>e</sup>* is directly derived.

 

**3. Numerical discretization of the full system and solving methodologies** 

 *VV V W dV F n dS dV*

1,6 ( ) *dW V Fn S V*

The unknown vector *W* , the flux vector ( *F* ) and the source term vector ( ) used in the

*t*

so they can be simply discretized following a finite volume approach:

*dt*

finite volume discretization are reported below.

In the previous parts the fundamental equations governing an MFD plasma flow have been described. They constitute a system of 11 scalar equations plus one scalar equation for each species considered and one scalar equation for each vibrating species considered. In the following the numerical strategy and methods adopted to solve the systems will be exposed All the equations briefly reported in the previous section have the following conservative

111

 

111

222 222 222

*zx y zy x z zyx zyx zyx*

*zyx zyx zyx*

, with a matrix pivoting technique it can be obtained the relation between electric

*y cz z cy*

 

> 

 

*c xyz ie*

*<sup>q</sup> BBB*

2

(3.1)

(3.2)

1

 

 

(2.21)

*T*

has been considered that

2 2

(2.22)

*Dm n <sup>E</sup> ne x*

*Dm n <sup>E</sup> ne y*

*x*

*y*

*z*

2 2

2 2

*Dm n <sup>E</sup> ne z*

( )

( )

*m*

gives the Hall parameter vector that represent the Hall parameter along all the

*x cy y cx*



A Magneto-Fluid-Dynamic Model and

Fig. 3. System solving strategy

intensity.

implemented:

the solution procedure to solve the system is: The fluid-dynamic subsystem is solved.

**magneto-fluid-dynamic problems treatment** 

coupled with an implicit evaluation of the source terms.

gases including ionized species (Arrhenius formulation).

(keeping E and B constant).

Computational Solving Methodologies for Aerospace Applications 133

So, given the initial conditions (both for fluid-dynamic variables and electromagnetic fields)

 Maxwell equations are solved, using as input the value of current density which has been found at the end of step one. This leads to update the electromagnetic field

 The newly found values for E and B are sent back again to the fluid dynamic over to update the residual of the equations, which perform another step like the first one

These three steps make the solution march in time for an interval *dt*. They have to be

CIRA H3NS code numerical structure has been used as starting platform. H3NS is a RANS structured multi-block finite volume solver that allows for the treatment of a wide range of compressible fluid dynamics problems considering air as a working gas. The fluid can be treated considering air as a prefect gas or as a mixture of perfect gases in of thermo-chemical non equilibrium (5 species air and 3 vibrational temperatures). With respect to the numerical formulation, conservation equations are written in integral form, and discretized with a finite volume, cell centred technique. Eulerian fluxes are computed with a Flux Difference Splitting method (Borrelli and Pandolfi 1990). Second order formulation is obtained by means of an Essentially Non Oscillatory reconstruction of interface value. Viscous fluxes are computed with a classical centred scheme. Time integration is performed by employing an explicit multistage Runge-Kutta algorithm

More accurate description of the models implemented and validation tests in 2D and 3D hypersonic problems could be found in (Schettino et al. 2008), (Battista et al. 2007). In order to be able to treat magneto fluid dynamic problems the following upgrades have been

Capability in the treatment of a generic multi-component reacting mixtureof perfect

**3.2 Fluid-dynamic subsystem: Starting platform description and upgrades through** 

repeated iteratively until the required total time of simulation is reached.

#### **3.1 Numerical solving strategy description: Weak coupling**

From the vectorial form of the system of equations we have a synthetic view of what kind of system we have to solve. It is also clear that it includes terms of different physical nature, in particular it combines acoustic and electromagnetic terms. In the present analysis it has been chosen to split this system of equations into the two following sub-systems:


This strategy, named from literature (Giordano and D'Ambrosio 2005) 'weak coupling', allows us to treat separately the Maxwell equations and then to "freeze" the resulting electromagnetic field and to solve the fluid-dynamic sub-system.

The main reason behind the choice of this procedure is the highly different characteristic speed of propagation between sound and light. The fluid-dynamic phenomena (accounted for in the first sub-system) propagate at the speed of sound, while the electromagnetic phenomena (described by the Maxwell equations) propagate at the speed of light. Since light travels at least one hundred thousand times faster than sound, it is advisable to split the main system. Moreover, the differences of signals propagation in space make crucial the adoption of different numerical methods to solve each sub-system. Since the equations are solved both in space and time, one of the most striking differences between these two numerical methods will be the choice of the integration time interval, dt. Much shorter intervals are requested when solving Maxwell equations, because of the higher propagation speed for the phenomena involved. The numerical method for the Maxwell equations solver will be investigated in the next section. Another issue to be faced with is how these two subsystems communicate each other. It is apparent that they are coupled, because the electromagnetic fields appear on the right hand side of the momentum and energy conservation equations, while the fluid velocity and the electric transport properties are directly linked to the current density (which plays an important role in the second Maxwell equation). Therefore it is necessary to find a way to solve the sub-systems alternately. The algorithm to be followed is sketched in the figure below:

Fig. 3. System solving strategy

0

 

*xy yy yz n xy yy yz U y i r y i*

*uvw J q*

 

, ,

,

*visc z*

*F*

, , 1

(3.3)

,1

*z*

*J*

0

 

> *xz yz zz n zx zy zz U z i r z i*

*uvw J q*

 

, ,,

*vz N v N z N*

*q eJ*

*v vv*

, 1 ,1 1 ,1

*z n vz v z*

...

*J q eJ*

1

,1

.. ..

*y*

*J*

, ,,

*vy N v N y N*

*q eJ*

From the vectorial form of the system of equations we have a synthetic view of what kind of system we have to solve. It is also clear that it includes terms of different physical nature, in particular it combines acoustic and electromagnetic terms. In the present analysis it has been

This strategy, named from literature (Giordano and D'Ambrosio 2005) 'weak coupling', allows us to treat separately the Maxwell equations and then to "freeze" the resulting

The main reason behind the choice of this procedure is the highly different characteristic speed of propagation between sound and light. The fluid-dynamic phenomena (accounted for in the first sub-system) propagate at the speed of sound, while the electromagnetic phenomena (described by the Maxwell equations) propagate at the speed of light. Since light travels at least one hundred thousand times faster than sound, it is advisable to split the main system. Moreover, the differences of signals propagation in space make crucial the adoption of different numerical methods to solve each sub-system. Since the equations are solved both in space and time, one of the most striking differences between these two numerical methods will be the choice of the integration time interval, dt. Much shorter intervals are requested when solving Maxwell equations, because of the higher propagation speed for the phenomena involved. The numerical method for the Maxwell equations solver will be investigated in the next section. Another issue to be faced with is how these two subsystems communicate each other. It is apparent that they are coupled, because the electromagnetic fields appear on the right hand side of the momentum and energy conservation equations, while the fluid velocity and the electric transport properties are directly linked to the current density (which plays an important role in the second Maxwell equation). Therefore it is necessary to find a way to solve the sub-systems alternately. The

*v vv*

, 1 ,1 1 ,1

*y n vy v y*

...

*J q eJ*

, ,

,

*visc y*

*F*

**3.1 Numerical solving strategy description: Weak coupling** 

electromagnetic field and to solve the fluid-dynamic sub-system.

algorithm to be followed is sketched in the figure below:

chosen to split this system of equations into the two following sub-systems:

1

0

 

*xx xy xz n xx xy xz Ux i r x i*

*uvw J q*

 

,1

*x*

*J*

, ,,

*vx N v N x N*

*q eJ*

*v vv*

1. Fluid-Dynamics equations sub-system 2. Maxwell equations sub-system

, 1 ,1 1 ,1

*x n vx v x*

...

*J q eJ*

,

*visc x*

*F*

So, given the initial conditions (both for fluid-dynamic variables and electromagnetic fields) the solution procedure to solve the system is:


These three steps make the solution march in time for an interval *dt*. They have to be repeated iteratively until the required total time of simulation is reached.

#### **3.2 Fluid-dynamic subsystem: Starting platform description and upgrades through magneto-fluid-dynamic problems treatment**

CIRA H3NS code numerical structure has been used as starting platform. H3NS is a RANS structured multi-block finite volume solver that allows for the treatment of a wide range of compressible fluid dynamics problems considering air as a working gas. The fluid can be treated considering air as a prefect gas or as a mixture of perfect gases in of thermo-chemical non equilibrium (5 species air and 3 vibrational temperatures). With respect to the numerical formulation, conservation equations are written in integral form, and discretized with a finite volume, cell centred technique. Eulerian fluxes are computed with a Flux Difference Splitting method (Borrelli and Pandolfi 1990). Second order formulation is obtained by means of an Essentially Non Oscillatory reconstruction of interface value. Viscous fluxes are computed with a classical centred scheme. Time integration is performed by employing an explicit multistage Runge-Kutta algorithm coupled with an implicit evaluation of the source terms.

More accurate description of the models implemented and validation tests in 2D and 3D hypersonic problems could be found in (Schettino et al. 2008), (Battista et al. 2007). In order to be able to treat magneto fluid dynamic problems the following upgrades have been implemented:

 Capability in the treatment of a generic multi-component reacting mixtureof perfect gases including ionized species (Arrhenius formulation).

A Magneto-Fluid-Dynamic Model and

*t*

indicates the time interval considered. *MFIN* and *MFOUT* can be written as:

sides by the surface area A, yields:

1 2

equation (3.7) is evaluated at the time "k".

"x,y,z" Cartesian grid, with cubic cells):

[LLL]

*N nn*

*z*

*x y x y*

 : is used to indicate a spatial interval : is used to indicate a time interval

1 1 2 2

1 1 2 2

1 2

fashion:

Computational Solving Methodologies for Aerospace Applications 135

case. Discretizing equation (3.1), and looking to it is possible to write in a finite volume

*<sup>M</sup> KK K <sup>N</sup> MM M*

<sup>1</sup>

*N N N N*

*x Wx W*

Here the subscript "N" addresses the cell we are referring to, while the superscript "K"

*M MM F FF IN N N* <sup>1</sup>

Considering again equation (3.6), expanding to the first order all the terms with superscript different than "K" (i.e. the terms that are not evaluated at the instant "k") and dividing both

> 11 1 22 2

*W FF xF F <sup>W</sup> <sup>W</sup> t WW*

Here the superscript "k" has been dropped, since it is clear that each term in the above

For a 3D analysis this procedure leads to the following equation (assuming to work on a

1 1

*x y x x*

[R] [LLL] [RRR]

 

1 11

*y z*

*N NN*

*t tt F F WW W W Wx Wx*

*N N N*

[C] [R] [L]

*NN N*

*NN N*

*y y z N n N n N nn N N*

*F F <sup>F</sup> t t W W*

*Nn Nn N nn N nn N*

*t t FF F F <sup>t</sup>*

*x x x y x y*

, *n n <sup>x</sup> <sup>y</sup>* : are the numbers of cells along x and y apxis respectively

1 1 <sup>1</sup> 22 2

*x x*

*N nn x x*

*<sup>F</sup> t t W FF W z x*

*y y z z*

*Wy Wy*

*N n N n*

111 222

<sup>1</sup> 1 1 <sup>2</sup> ( ) ()

*IN OUT <sup>N</sup> N N <sup>W</sup> Ax F F A A x*

(3.6)

1 2

1 1

1 1

1 1

*x y*

*W z*

*N nn*

(3.7)

1

(3.8)

*N*

*<sup>t</sup> <sup>W</sup>*

*N N*

1 1 1 1

1 1

*x x*

1 1

*NN N*

*N N N N*

*M MM F FF OUT N N*


Models considered for these upgrades have been widely discussed in previous Sections.

#### **3.3 Numerical methodologies for Maxwell equations : Implicit and explicit Maxwell equation solver design**

In this section the numerical methods used to solve Maxwell equations are discussed. So far, two methods have been employed, one of them is implicit, the other one explicit. Each method has its own strengths and its drawbacks. In the two paragraphs below, they are explained and compared. The equations to be discretized are the Maxwell equations in conservative form ((2.19), (2.20)), i.e.:

$$\frac{\partial}{\partial t} \int\_{V} \left( \underline{\underline{V}^{M}} \right) dV + \int\_{\widetilde{\underline{\varepsilon}}V} \left( \underline{\underline{F}^{M}} \right) dS = \int\_{V} \underline{\underline{\Omega}}{V} \,dV \tag{3.4}$$

With:

$$\underline{\underline{\mathbf{V}}^{M}} = \begin{pmatrix} B\_{x} \\ B\_{y} \\ B\_{z} \\ E\_{x} \\ E\_{y} \\ E\_{z} \end{pmatrix} \quad \underline{\underline{F}\_{x}}^{M} = \begin{pmatrix} 0 \\ -E\_{z} \\ E\_{y} \\ 0 \\ -c^{2}B\_{z} \\ -c^{2}B\_{y} \end{pmatrix} \quad \underline{\underline{F}\_{y}}^{M} = \begin{pmatrix} E\_{x} \\ 0 \\ -E\_{x} \\ -c^{2}B\_{z} \\ 0 \\ c^{2}B\_{x} \end{pmatrix} \quad \underline{\underline{F}\_{x}}^{M} = \begin{pmatrix} -E\_{y} \\ E\_{x} \\ 0 \\ 0 \\ -c^{2}B\_{x} \\ -c^{2}B\_{x} \\ 0 \end{pmatrix} \quad \underline{\underline{\Omega}}^{M} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ I\_{qx} + \rho\_{c}u/\rho\_{\mathcal{O}} \\ I\_{qy} + \rho\_{c}v/\rho\_{\mathcal{O}} \\ I\_{qz} + \rho\_{c}w/\rho\_{\mathcal{O}} \end{pmatrix} \tag{3.5}$$

The superscript "M" reminds that this method is being applied to the Maxwell equations.

#### **3.3.1 Implicit methodology**

The great strength of all implicit numerical schemes is their great stability. Using an implicit method it is possible to choose a larger integration time interval without compromising the stability of the method itself. This is most relevant when dealing with signals which travel so fast as the electromagnetic waves.

Explicit schemes make time advancement very simple, but often suffer severe restrictions on time step due to the loss of stability. On the other side, implicit schemes generally have much better numerical properties in term of stability and so allow a larger time step, but require the solution of a system of equations to perform time integration. When analyzing a 3D problem as in the present case, the matrix associated to the system is too large to be inverted in a reasonable computational time. This means that an iterative method has to be applied to solve the system. So, the time saved using a larger time step is actually lost in solving a large linear system of equations at each iteration. To introduce the implicit method, a 1D simple case is described. Then, the result obtained will be extended to the 3D case. Discretizing equation (3.1), and looking to it is possible to write in a finite volume fashion:

$$\frac{\Delta \mathcal{W}\_N^M}{\Delta t} \text{(A\Delta x)} \quad = \left( \left( F\_{IN}^M \right)\_N^{K+1} - \left( F\_{OUT}^M \right)\_N^{K+1} \right) \text{A} \quad + \left( \Omega\_N^M \right)^{K+1} \text{\(A\Delta x\)} \tag{3.6}$$

Here the subscript "N" addresses the cell we are referring to, while the superscript "K" indicates the time interval considered.

*MFIN* and *MFOUT* can be written as:

134 Fluid Dynamics, Computational Modeling and Applications

Thermodynamic properties formulation in terms of Gordon-Mc Bride Polynomial fits.

Extension of thermodynamic properties (specific heat, enthalpy and entropy) to higher

Models considered for these upgrades have been widely discussed in previous Sections.

**3.3 Numerical methodologies for Maxwell equations : Implicit and explicit Maxwell** 

In this section the numerical methods used to solve Maxwell equations are discussed. So far, two methods have been employed, one of them is implicit, the other one explicit. Each method has its own strengths and its drawbacks. In the two paragraphs below, they are explained and compared. The equations to be discretized are the Maxwell equations in

> *M MM V VV*

*t*

2

*y*

*F*

*M x*

2

*c B*

The superscript "M" reminds that this method is being applied to the Maxwell equations.

0

*z*

2

*M z*

*F*

*x*

The great strength of all implicit numerical schemes is their great stability. Using an implicit method it is possible to choose a larger integration time interval without compromising the stability of the method itself. This is most relevant when dealing with signals which travel

Explicit schemes make time advancement very simple, but often suffer severe restrictions on time step due to the loss of stability. On the other side, implicit schemes generally have much better numerical properties in term of stability and so allow a larger time step, but require the solution of a system of equations to perform time integration. When analyzing a 3D problem as in the present case, the matrix associated to the system is too large to be inverted in a reasonable computational time. This means that an iterative method has to be applied to solve the system. So, the time saved using a larger time step is actually lost in solving a large linear system of equations at each iteration. To introduce the implicit method, a 1D simple case is described. Then, the result obtained will be extended to the 3D

0

*E*

*c B*

*E*

 

*z*

2 2

*c B c B*

*z y*

0

*<sup>y</sup> <sup>M</sup>*

0

 

*E E*

*z*

*W dV F dS dV*

(3.4)

2

*c B c B*

*y x* 0

*qx c qy c qz c*

*J u J v J w*

*M*

0 0 0

0 0

(3.5)

0

0

*y x*

*E E*

 

Electromagnetic terms added at the equations implemented

(Gordon and McBride 1971)

conservative form ((2.19), (2.20)), i.e.:

*x y*

*B B B*

 

*E E E*

**3.3.1 Implicit methodology** 

so fast as the electromagnetic waves.

*x*

*F*

*zM x y z*

*W*

Included ionized species transport properties.

Effects of electromagnetic field on transport accounted.

temperatures.

**equation solver design** 

With:

$$F\_{IN}^{M} = \frac{1}{2} \left( F^{M}{}\_{N} + F^{M}{}\_{N-1} \right) \qquad \qquad \qquad F\_{OUT}^{M} = \frac{1}{2} \left( F^{M}{}\_{N} + F^{M}{}\_{N+1} \right)$$

Considering again equation (3.6), expanding to the first order all the terms with superscript different than "K" (i.e. the terms that are not evaluated at the instant "k") and dividing both sides by the surface area A, yields:

 1 1 1 1 1 1 1 1 11 1 22 2 1 2 *N N N N N N N N N N N N N W FF xF F <sup>W</sup> <sup>W</sup> t WW x Wx W* (3.7)

Here the superscript "k" has been dropped, since it is clear that each term in the above equation (3.7) is evaluated at the time "k".

For a 3D analysis this procedure leads to the following equation (assuming to work on a "x,y,z" Cartesian grid, with cubic cells):

1 1 1 1 1 1 1 1 [C] [R] [L] 1 1 <sup>1</sup> 22 2 [R] [LLL] [RRR] 111 222 *x y x x x x x x NN N NN N NN N y y z N n N n N nn N N N n N n t tt F F WW W W Wx Wx F F <sup>F</sup> t t W W Wy Wy* 1 1 11 [LLL] 1 1 2 2 1 1 2 2 *x y x y x y x x x y x y N N nn z N nn x x N NN N nn y y z z Nn Nn N nn N nn N <sup>t</sup> <sup>W</sup> W z <sup>F</sup> t t W FF W z x t t FF F F <sup>t</sup> y z* (3.8)

: is used to indicate a spatial interval

: is used to indicate a time interval

, *n n <sup>x</sup> <sup>y</sup>* : are the numbers of cells along x and y apxis respectively

A Magneto-Fluid-Dynamic Model and

<sup>0</sup> <sup>50</sup> <sup>100</sup> <sup>150</sup> <sup>200</sup> <sup>250</sup> <sup>300</sup> <sup>350</sup> -2

<sup>0</sup> <sup>50</sup> <sup>100</sup> <sup>150</sup> <sup>200</sup> <sup>250</sup> <sup>300</sup> <sup>350</sup> -2

( *t* ) of 10-5 s.

0

0

2

4

6

8

10 12

2

4

6

8

10 12

Computational Solving Methodologies for Aerospace Applications 137

Although here the solution, *W* , can be easily obtained, the trouble with this method regards its stability. So the CFL (Courant–Friedrichs–Lewy) number must be chosen

Considering the very high speed of signal propagation ( <sup>8</sup> *c* 3 10 m/s ), this means that the time step must be extremely short. However, this second method turns out to be faster than the implicit one due to its simplicity. The second order explicit numerical method used in this work and reported in 1D formulation in Eqn.3.11 (again, N represent the space step and K the time step) has been compared with literature numerical methods like Lax-Friederichs and Lax-Wendroff methods (Chung 2010), for the solution of the classical advection equation test consisting in the propagation of a square and a smooth wave, using as conditions a wave speed of 60 m/s, a grid spacing ( *x* ) of 10-3 m and a time step

<sup>1</sup>

0

0

2

4

6

8

10 12

Fig. 4. Lax-Friederichs results after 40 timesteps (left) and 400 timesteps (right)

Fig. 5. Lax--Wendroff results after 40 timesteps (left) and 400 timesteps (right)

2

4

6

8

10 12

*K KK K K K K K N NN N N N N N t q qq q q q q q x*

1 1 3 2

11 1212

<sup>0</sup> <sup>50</sup> <sup>100</sup> <sup>150</sup> <sup>200</sup> <sup>250</sup> <sup>300</sup> <sup>350</sup> -2

<sup>0</sup> <sup>50</sup> <sup>100</sup> <sup>150</sup> <sup>200</sup> <sup>250</sup> <sup>300</sup> <sup>350</sup> -2

(3.12)

 

carefully in order to have an acceptable time step and to grant the method stability.

[C], [L], [R], etc… are 6x6 square matrixes (remind that *F* and *W* are six-elements vectors). This equation is just a system of *ntot* scalar equations, where *ntot* ( *nnn <sup>x</sup> <sup>y</sup> <sup>z</sup>* ) is the total number of grid cells. It can be written in the following more compact form:

$$\left[D\right] \cdot \left[\tilde{\Delta}\mathcal{W}\right] = -\frac{1}{2} \left\{\Delta F^{x}\right\} \frac{\Delta t}{\Delta x} - \frac{1}{2} \left\{\Delta F^{y}\right\} \frac{\Delta t}{\Delta y} - \frac{1}{2} \left\{\Delta F^{z}\right\} \frac{\Delta t}{\Delta z} + \left\{\Omega\right\} \Delta t \tag{3.9}$$

Where:

[ ] [ ] 0 0 ... [ ] 0 0 0 ... [ ] 0 [ ] [ ] [ ] 0 ... 0 [ ] 0 0 ... 0 [ ] 0 0 [ ] 0 ... 0 0 [ ] 0 ... 0 0 0 0 0 [ ] ... 0 0 0 [ ] ... 0 0 ... 0 0 ... ... ... ... ... ... ... ... 0 [ ] 0 0 0 0 [ ] 0 0 0 ... [ ] 0 [ ] 0 [ ] 0 0 0 0 [ ] 0 0 ... 0 [ ] 0 0 [ ] 0 0 0 0 [ ] 0 ... 0 0 ... .. *C R RR RRR LCR RR RRR C RR C RR LL C RR D LL C RR LL C* 0 0 [ ] 0 0 0 0 [ ] . ... ... ... ... ... ... ... ... ... ... ... [ ] 0 0 0 [ ] 0 0 0 0 0 [] 0 0 0 [ ] 0 0 0 [ ] 0 0 0 0 0 [] [] 0 0 [ ] 0 0 0 [ ] 0 0 0 0 [] [] *RRR RR LLL LL C LLL LL C R LLL LL L C* (nxny) (nxnynz) (nxnynz) (nxny)

This system has to be solved at each time step. The unknowns are the temporal increments *W* , whose knowledge is necessary to update the vector W (containing the intensities of the electromagnetic fields) all through the spatial domain of integration. [D] is a seven diagonal block matrix. Given the considerable dimensions of the matrix [D] its inversion is very expensive in terms of computational time. Therefore, an iterative method is required to solve the system, i.e. it is strongly necessary to use a proper pre-conditioner to achieve a faster convergence to the solution.

#### **3.3.2 Explicit methodology**

An explicit method allows to directly find the unknown values for electromagnetic field at each time step, with no need of solving a large system of equations. The explicit method employed is "centred in space" and "forward in time". As a consequence, the equation (3.1) can be discretized as follows:

$$\frac{\tilde{\Delta}\mathcal{V}\_{N}^{M}}{\Delta t} \Delta V = \left( \left( F\_{IN}^{M} \right)\_{N}^{K} - \left( F\_{OUT}^{M} \right)\_{N}^{K} \right) A + \left( \Omega\_{N}^{M} \right)^{K} \Delta V \tag{3.10}$$

Each term is evaluated at instant "k". Rearranging the previous equation for the simple case of cubic cells:

$$\begin{split} \tilde{\Delta}\mathcal{W}\_{\text{N}} &= \frac{\Delta t}{\Delta V} \frac{1}{2} \Big[ \Big( F\_{\text{N}+1}^{\text{x}} - F\_{\text{N}-1}^{\text{x}} \Big) A\_{\text{x}} + \Big( F\_{\text{N}+\text{n}\_{x}}^{\text{y}} - F\_{\text{N}-\text{n}\_{x}}^{\text{y}} \Big) A\_{\text{y}} + \Big( F\_{\text{N}+\text{n}\_{x}\text{n}\_{z}}^{\text{z}} - F\_{\text{N}-\text{n}\_{x}\text{n}\_{z}}^{\text{z}} \Big) A\_{z} \Big] + \\ &+ \Big( \Omega\_{\text{N}} \Big) \Delta t \end{split} \tag{3.11}$$

[C], [L], [R], etc… are 6x6 square matrixes (remind that *F* and *W* are six-elements vectors). This equation is just a system of *ntot* scalar equations, where *ntot* ( *nnn <sup>x</sup> <sup>y</sup> <sup>z</sup>* ) is the total

> <sup>111</sup> [ ] <sup>222</sup> *x z ttt <sup>y</sup> DW F F F t*

[ ] [ ] 0 0 ... [ ] 0 0 0 ... [ ] 0 [ ] [ ] [ ] 0 ... 0 [ ] 0 0 ... 0 [ ] 0 0 [ ] 0 ... 0 0 [ ] 0 ... 0 0 0 0 0 [ ] ... 0 0 0 [ ] ... 0 0 ... 0 0 ... ... ... ... ... ... ... ... 0 [ ] 0 0 0 0 [ ] 0 0 0 ... [ ] 0 [ ] 0 [ ] 0 0 0 0 [ ] 0 0 ... 0 [ ] 0 0 [ ] 0 0 0 0 [ ] 0 ... 0 0

*C RR*

*LL C RR*

[ ] 0 0 0 [ ] 0 0 0 0 0 [] 0 0 0 [ ] 0 0 0 [ ] 0 0 0 0 0 [] [] 0 0 [ ] 0 0 0 [ ] 0 0 0 0 [] []

. ... ... ... ... ... ... ... ... ... ... ...

*LLL LL C R LLL LL L C*

*C R RR RRR LCR RR RRR*

*C RR*

(nxny)

*LL C RR*

*LL C*

*LLL LL C*

This system has to be solved at each time step. The unknowns are the temporal increments *W* , whose knowledge is necessary to update the vector W (containing the intensities of the electromagnetic fields) all through the spatial domain of integration. [D] is a seven diagonal block matrix. Given the considerable dimensions of the matrix [D] its inversion is very expensive in terms of computational time. Therefore, an iterative method is required to solve the system, i.e. it is strongly necessary to use a proper pre-conditioner to achieve a

An explicit method allows to directly find the unknown values for electromagnetic field at each time step, with no need of solving a large system of equations. The explicit method employed is "centred in space" and "forward in time". As a consequence, the equation (3.1)

> *<sup>M</sup> KK K <sup>N</sup> MM M*

*<sup>W</sup> VF F A V*

Each term is evaluated at instant "k". Rearranging the previous equation for the simple case

1 2 *x x x z x z x x y y z z N N Nx Nn Nn y N nn N nn z*

*<sup>t</sup> W F F AF F AF F A*

*IN OUT <sup>N</sup> N N*

 

(3.10)

0 0 [ ] 0 0 0 0 [ ]

(nxnynz)

*RRR*

*RR*

(3.11)

*xyz*

(3.9)

number of grid cells. It can be written in the following more compact form:

Where:

faster convergence to the solution.

*D*

(nxny)

(nxnynz)

**3.3.2 Explicit methodology** 

can be discretized as follows:

of cubic cells:

*N*

*V t*

*t*

1 1

... ..

Although here the solution, *W* , can be easily obtained, the trouble with this method regards its stability. So the CFL (Courant–Friedrichs–Lewy) number must be chosen carefully in order to have an acceptable time step and to grant the method stability.

Considering the very high speed of signal propagation ( <sup>8</sup> *c* 3 10 m/s ), this means that the time step must be extremely short. However, this second method turns out to be faster than the implicit one due to its simplicity. The second order explicit numerical method used in this work and reported in 1D formulation in Eqn.3.11 (again, N represent the space step and K the time step) has been compared with literature numerical methods like Lax-Friederichs and Lax-Wendroff methods (Chung 2010), for the solution of the classical advection equation test consisting in the propagation of a square and a smooth wave, using as conditions a wave speed of 60 m/s, a grid spacing ( *x* ) of 10-3 m and a time step ( *t* ) of 10-5 s.

$$q\_N^{K+1} = \frac{1}{3} \left[ q\_N^K + q\_{N-1}^K + q\_{N+1}^K + \frac{1}{2} \frac{\Delta t}{\Delta x} (q\_{N+1}^K + q\_{N+2}^K - q\_{N-1}^K - q\_{N+2}^K) \right] \tag{3.12}$$

Fig. 4. Lax-Friederichs results after 40 timesteps (left) and 400 timesteps (right)

Fig. 5. Lax--Wendroff results after 40 timesteps (left) and 400 timesteps (right)

A Magneto-Fluid-Dynamic Model and

called ghost cell; K indicates the time step considered.

reflection has completely messed up the solution.

following one:

*t* , time step *x* , grid spacing

equal to *c*.

iterations, the solution is clean.

1

Computational Solving Methodologies for Aerospace Applications 139

Initially the condition adopted on the fluxes crossing the boundary cells was simply the

*K K F F N N* , where the subscript N indicates the last cell of the domain and N+1 is the so

So the flux associated to the ghost cell was exactly the one associated to the last cell of the domain. Of course this leads to wave reflection, as it is apparent in the following pictures where the magnetic field generated by a wire is represented. After 200 iterations the wave has already reached the boundary and starts to be reflected, after 1000 iterations the wave

To avoid wave reflections, Orlanski (1976) conditions have been adopted. These are special conditions that define the flux through the boundaries considering what happened at previous time instants. So the boundary flux at instant *tN* is dependent on the fluxes in the

1 2

*<sup>x</sup> <sup>x</sup> FFF*

Some complication can be found if the propagation velocity is unknown, but dealing with electromagnetic waves there is no such problem. The propagation velocity is constant and

Again the previous case of a magnetic field generated by a wire is represented in the next picture. Only the boundary conditions have been changed and this time there is no reflection at all. The waves are absorbed and in both cases, after 200 and after 1000

Fig. 8. Magnetic field Bz generated by a current flowing along the "x" direction, after 200 time steps (left) and after 1000 time steps (right). Applying Orlanski conditions there is no

reflection and the situation is unchanged in the two cases.

*<sup>t</sup> <sup>t</sup> <sup>c</sup> <sup>c</sup>*

*t t c c x x*

1 1 *KKK NNN*

2 1

(3.13)

nearby cells at instants *<sup>N</sup>* <sup>1</sup> *t* and *<sup>N</sup>* <sup>2</sup> *t* . The condition can be written as:

1 1

Fig. 6. Results with the method suggested in the present work results after 40 timesteps (left) and 400 timesteps (right)

The results show that the proposed numerical method is comparable with the literature ones in terms of quality of the result; moreover, due to its simplicity it is faster than other more complex methods available in literature (Chung 2010).

#### **3.3.3 Orlanski condition for free boundaries**

When dealing with electromagnetic waves in a finite domain, one of the most complicated problems is to set a proper free boundary condition in order to avoid wave reflections. An option could be to extend the computational domain far enough from the electromagnetic wave sources. This is not acceptable, because of the cost in terms of computational time and storage memory. A 3D code like the one developed in this work (HOPE) is already much demanding and it is important to keep the total number of cells as low as possible. Besides, if the time span simulated is not very short, the waves will reach the boundary anyway and they will be reflected back into the domain (moreover the grid has to be the same for the fluid dynamic part and it is impossible to consider to solve NS multi specie equation in a too much extended domain).

Fig. 7. Magnetic field Bz generated by a current flowing along the "x" direction, after 200 time steps (left after 1000 time steps (right). The electromagnetic wave is reflected at the boundary and the solution is completely corrupted after 1000 time steps.

Initially the condition adopted on the fluxes crossing the boundary cells was simply the following one:

1 *K K F F N N* , where the subscript N indicates the last cell of the domain and N+1 is the so called ghost cell; K indicates the time step considered.

So the flux associated to the ghost cell was exactly the one associated to the last cell of the domain. Of course this leads to wave reflection, as it is apparent in the following pictures where the magnetic field generated by a wire is represented. After 200 iterations the wave has already reached the boundary and starts to be reflected, after 1000 iterations the wave reflection has completely messed up the solution.

To avoid wave reflections, Orlanski (1976) conditions have been adopted. These are special conditions that define the flux through the boundaries considering what happened at previous time instants. So the boundary flux at instant *tN* is dependent on the fluxes in the nearby cells at instants *<sup>N</sup>* <sup>1</sup> *t* and *<sup>N</sup>* <sup>2</sup> *t* . The condition can be written as:

$$F\_{N+1}^{K} = \frac{\left(1 - \frac{\Delta t}{\Delta \mathbf{x}} c\right)}{\left(1 + \frac{\Delta t}{\Delta \mathbf{x}} c\right)} F\_{N+1}^{K-2} + \frac{2\frac{\Delta t}{\Delta \mathbf{x}} c}{\left(1 + \frac{\Delta t}{\Delta \mathbf{x}} c\right)} F\_{N}^{K-1} \tag{3.13}$$

*t* , time step

138 Fluid Dynamics, Computational Modeling and Applications

0

Fig. 6. Results with the method suggested in the present work results after 40 timesteps (left)

The results show that the proposed numerical method is comparable with the literature ones in terms of quality of the result; moreover, due to its simplicity it is faster than other more

When dealing with electromagnetic waves in a finite domain, one of the most complicated problems is to set a proper free boundary condition in order to avoid wave reflections. An option could be to extend the computational domain far enough from the electromagnetic wave sources. This is not acceptable, because of the cost in terms of computational time and storage memory. A 3D code like the one developed in this work (HOPE) is already much demanding and it is important to keep the total number of cells as low as possible. Besides, if the time span simulated is not very short, the waves will reach the boundary anyway and they will be reflected back into the domain (moreover the grid has to be the same for the fluid dynamic part and it is impossible to consider to solve NS multi specie equation in a too

 Fig. 7. Magnetic field Bz generated by a current flowing along the "x" direction, after 200 time steps (left after 1000 time steps (right). The electromagnetic wave is reflected at the

boundary and the solution is completely corrupted after 1000 time steps.

2

4

6

8

10 12

<sup>0</sup> <sup>50</sup> <sup>100</sup> <sup>150</sup> <sup>200</sup> <sup>250</sup> <sup>300</sup> <sup>350</sup> -2

<sup>0</sup> <sup>50</sup> <sup>100</sup> <sup>150</sup> <sup>200</sup> <sup>250</sup> <sup>300</sup> <sup>350</sup> -2

complex methods available in literature (Chung 2010).

**3.3.3 Orlanski condition for free boundaries** 

and 400 timesteps (right)

much extended domain).

0

2

4

6

8

10 12

*x* , grid spacing

Some complication can be found if the propagation velocity is unknown, but dealing with electromagnetic waves there is no such problem. The propagation velocity is constant and equal to *c*.

Again the previous case of a magnetic field generated by a wire is represented in the next picture. Only the boundary conditions have been changed and this time there is no reflection at all. The waves are absorbed and in both cases, after 200 and after 1000 iterations, the solution is clean.

Fig. 8. Magnetic field Bz generated by a current flowing along the "x" direction, after 200 time steps (left) and after 1000 time steps (right). Applying Orlanski conditions there is no reflection and the situation is unchanged in the two cases.

A Magneto-Fluid-Dynamic Model and

reproduce experimental ne behaviour.

Fig. 11. Test 2: Mach number and total enthalpy distributions

Fig. 12. Test 2 Electron number density comparison with experimental data

**nozzle** 

Computational Solving Methodologies for Aerospace Applications 141

In the test campaign described in Dunn and Lordi (1969), electron temperatures and electron densities have been measured in the Calspan nozzle; the working gas was air and the equilibrium reservoir conditions were T0=6850 K and P0=0.94738 MPa. Geometry and grid used in this test is reported in Battista (2009). In Fig. 12 are reported the comparison of numerical data obtained with different kinetic schemes with experimental data. It is evident that the scheme proposed here is the one, among the 11 species schemes, that better

**4.2 Validation test 2 - Ionization chemistry in an expansion experiment: Calspan** 

#### **4. Models validation numerical tests**

A simple validation test for different models was presented in Battista et al. (2008, 2009). Hereinafter three main validation tests will be presented.

#### **4.1 Validation test 1 - Ionization chemistry in re-entry experiment: RAM-C II**

The aim of the RAM-C program was to obtain a better understanding of the factors that influence transmission of radio waves through plasmas and to search for methods to reduce or eliminate blackout (Schexnayder et al 1977). In this frame some experimental data have been collected about the electrons number density at different flight speed and altitudes. The RAM-C test has been considered here in order to compare numerical results in terms of electron number density with the experimental data. For our scopes, the trajectory point characterized by an altitude of 70 km and a re-entry speed of 7654 m/s, corresponding to a Mach number of about 26, has been considered.

Fig. 9. Geometry and payload configuration of the RAM-C II test with highlighted electron density measurement stations considered in this work.

The results obtained with the 11 species model proposed in this work are shown in the following Fig. 10 and compared with experimental data in terms of electron number density.

Fig. 10. Numerical and experimental data comparison at different RAM-C stations.

A simple validation test for different models was presented in Battista et al. (2008, 2009).

The aim of the RAM-C program was to obtain a better understanding of the factors that influence transmission of radio waves through plasmas and to search for methods to reduce or eliminate blackout (Schexnayder et al 1977). In this frame some experimental data have been collected about the electrons number density at different flight speed and altitudes. The RAM-C test has been considered here in order to compare numerical results in terms of electron number density with the experimental data. For our scopes, the trajectory point characterized by an altitude of 70 km and a re-entry speed of 7654 m/s, corresponding to a

Fig. 9. Geometry and payload configuration of the RAM-C II test with highlighted electron

The results obtained with the 11 species model proposed in this work are shown in the following Fig. 10 and compared with experimental data in terms of electron number density.

Fig. 10. Numerical and experimental data comparison at different RAM-C stations.

**4.1 Validation test 1 - Ionization chemistry in re-entry experiment: RAM-C II** 

**4. Models validation numerical tests** 

Mach number of about 26, has been considered.

density measurement stations considered in this work.

Hereinafter three main validation tests will be presented.

#### **4.2 Validation test 2 - Ionization chemistry in an expansion experiment: Calspan nozzle**

In the test campaign described in Dunn and Lordi (1969), electron temperatures and electron densities have been measured in the Calspan nozzle; the working gas was air and the equilibrium reservoir conditions were T0=6850 K and P0=0.94738 MPa. Geometry and grid used in this test is reported in Battista (2009). In Fig. 12 are reported the comparison of numerical data obtained with different kinetic schemes with experimental data. It is evident that the scheme proposed here is the one, among the 11 species schemes, that better reproduce experimental ne behaviour.

Fig. 11. Test 2: Mach number and total enthalpy distributions

Fig. 12. Test 2 Electron number density comparison with experimental data

A Magneto-Fluid-Dynamic Model and

facility using Mach 6 nozzle .

Fig. 14. Geometry of the test body

Fig. 15. Computational domain and grid

computational grid characterized by 78x60 cells is depicted.

Computational Solving Methodologies for Aerospace Applications 143

test body with a half-vertex angle of 22.5 deg and the maximum diameter of 40 mm. A cylindrical, 40 mm diameter section is placed behind the cone. Two pressure sensors are embedded in the iron elements between the magnets, one on the cone surface and one in the cylindrical expansion region. The test geometry is reported in Fig. 14. The free stream conditions considered are directly obtained from the nozzle argon viscous test reported in Battista (2009), since the experiment considered here has been carried out in ALTA Heat

Calculations have been carried out on the domain shown in Fig. 15 where also the

Magnetic field is given as boundary condition on the wall of the test article and has been deduced from Cristofolini (2006). For the outer boundary the Orlanski boundary condition has been used except for the face on the symmetry plane where the proper boundary

condition has been applied. Results in terms magnetic field are reported in Fig. 16.

#### **4.3 Validation test 3 - Coupled simulation in argon: ALTA test**

The basic concepts of magneto fluid-dynamic interaction for plasma flow control in the shock layer of an hypersonic vehicle are here briefly described. Referring to the axissymmetric blunt body configuration shown in Fig. 13, assuming that a part of the gas is ionized, and that some device embedded within the body generates the B field (amber lines), free charges are subject to a force per charge units *u B* .

$$
\underline{F} = q \left( \underline{E} + \underline{\mu} \times \underline{B} \right) \tag{4.1}
$$

This force exerted upon the charged particles in the shock layer tends to drive them in the direction *u B* , and (due to a different collisional behaviour between ions and electrons) on a macroscopic scale generates a current density flowing orthogonally to the velocity field direction (Faraday currents).

The current density can be expressed by the generalized Ohm law:

$$
\underline{J} = \overline{\sigma\left(\underline{E} + \underline{\mu} \times \underline{B}\right)}\tag{4.2}
$$

The Faraday current density generates the *J B* body force represented with red arrows. These forces act in a direction opposite to the flow. Furthermore, since conductivity in electromagnetic field assumes a tensorial form, Hall currents (due to Hall collisions) arise orthogonally to B and *u B* directions and weaken Faraday currents as well as *J B* force.

Fig. 13. Scheme of the interaction over a blunt body

Such concepts are used for active flow control in the experiment considered here and used for models validation.

The test considered here is the one described in Cristofolini et al. (2006) and numerically investigated. It consists in an array of the magnets that has been assembled to form a conical

The basic concepts of magneto fluid-dynamic interaction for plasma flow control in the shock layer of an hypersonic vehicle are here briefly described. Referring to the axissymmetric blunt body configuration shown in Fig. 13, assuming that a part of the gas is ionized, and that some device embedded within the body generates the B field (amber

This force exerted upon the charged particles in the shock layer tends to drive them in the direction *u B* , and (due to a different collisional behaviour between ions and electrons) on a macroscopic scale generates a current density flowing orthogonally to the velocity field

> *J EuB*

The Faraday current density generates the *J B* body force represented with red arrows. These forces act in a direction opposite to the flow. Furthermore, since conductivity in electromagnetic field assumes a tensorial form, Hall currents (due to Hall collisions) arise orthogonally to B and *u B* directions and weaken Faraday currents as well as *J B* force.

Such concepts are used for active flow control in the experiment considered here and used

The test considered here is the one described in Cristofolini et al. (2006) and numerically investigated. It consists in an array of the magnets that has been assembled to form a conical

*F qEuB* (4.1)

(4.2)

**4.3 Validation test 3 - Coupled simulation in argon: ALTA test** 

lines), free charges are subject to a force per charge units *u B* .

The current density can be expressed by the generalized Ohm law:

Fig. 13. Scheme of the interaction over a blunt body

for models validation.

direction (Faraday currents).

test body with a half-vertex angle of 22.5 deg and the maximum diameter of 40 mm. A cylindrical, 40 mm diameter section is placed behind the cone. Two pressure sensors are embedded in the iron elements between the magnets, one on the cone surface and one in the cylindrical expansion region. The test geometry is reported in Fig. 14. The free stream conditions considered are directly obtained from the nozzle argon viscous test reported in Battista (2009), since the experiment considered here has been carried out in ALTA Heat facility using Mach 6 nozzle .

Fig. 14. Geometry of the test body

Calculations have been carried out on the domain shown in Fig. 15 where also the computational grid characterized by 78x60 cells is depicted.

Fig. 15. Computational domain and grid

Magnetic field is given as boundary condition on the wall of the test article and has been deduced from Cristofolini (2006). For the outer boundary the Orlanski boundary condition has been used except for the face on the symmetry plane where the proper boundary condition has been applied. Results in terms magnetic field are reported in Fig. 16.

A Magneto-Fluid-Dynamic Model and

**Pressure**

200

**5. Conclusions** 

implicit one.

applications.

Fig. 18. Pressure along the body wall

400

600

800

1000

1200

 **(Pa)**

Computational Solving Methodologies for Aerospace Applications 145

**B CFD NO B CFD B EXP no B exp**

**x(mm)**

0 20 40 60 80

A three-dimensional model for magnetofluiddynamic problems solution has been proposed considering plasma as a single fluid in thermo-chemical non equilibrium. For chemistry, transport and thermodynamic treatment a macroscopic approach has been followed that allows for the treatment of a general mixture of perfect gases including ionized species. This, in turn, permits not only to deal with non-equilibrium air problems, but gives the opportunity to consider seeded flows or extraterrestrial atmospheric flows in presence of ionization. Models for electrical conductivity have been proposed for air depending upon gas composition and temperature, accounting also the effect of charged particles collisions that could be important in some cases. For what concerns numerical solving strategy a loosely coupling technique to solve full system of equations (Navier-Stokes and Maxwell) has been extended to 3D problems. Both numerical methods for Maxwell equations have been developed and tested, in particular both time marching explicit and implicit methods have been considered and tested for the solution of Maxwell equations. Actually explicit methodology is more rapidly reaching convergence than the

Test validation results are in good agreement with available experimental data, and shows that effect of magneto-fluid-dynamic interaction could be relevant for aerospace

Currently CIRA and ALTA with Bologna University Electrical Department are working together in the design of MFD interaction experiments in air in order to understand the real opportunity of application of these technology to Earth re-entry (Cristofolini et al 2010).

Fig. 16. Computed magnetic field lines.

In Fig. 17 has been shown the comparison in terms of pressure between magnetic and nonmagnetic case, in Fig. 18 the comparison with experimental data in terms of pressure shows a quite good agreement; it has to be remarked that the MHD interaction has a strong effect on body pressure distribution, this could be used for an active control scopes.

Fig. 17. Computed pressure field without and with MHD interaction

Fig. 18. Pressure along the body wall

#### **5. Conclusions**

144 Fluid Dynamics, Computational Modeling and Applications

**y(m)**

0.01

0.02

0.03

0.04

0.05

0.06

0.07 BX (T)

Fig. 16. Computed magnetic field lines.

0.4 0.25 0.1 6.08653E-08 0 -0.15 -0.3 -0.45

**x(m)**


In Fig. 17 has been shown the comparison in terms of pressure between magnetic and nonmagnetic case, in Fig. 18 the comparison with experimental data in terms of pressure shows a quite good agreement; it has to be remarked that the MHD interaction has a strong effect

on body pressure distribution, this could be used for an active control scopes.

Fig. 17. Computed pressure field without and with MHD interaction

A three-dimensional model for magnetofluiddynamic problems solution has been proposed considering plasma as a single fluid in thermo-chemical non equilibrium. For chemistry, transport and thermodynamic treatment a macroscopic approach has been followed that allows for the treatment of a general mixture of perfect gases including ionized species. This, in turn, permits not only to deal with non-equilibrium air problems, but gives the opportunity to consider seeded flows or extraterrestrial atmospheric flows in presence of ionization. Models for electrical conductivity have been proposed for air depending upon gas composition and temperature, accounting also the effect of charged particles collisions that could be important in some cases. For what concerns numerical solving strategy a loosely coupling technique to solve full system of equations (Navier-Stokes and Maxwell) has been extended to 3D problems. Both numerical methods for Maxwell equations have been developed and tested, in particular both time marching explicit and implicit methods have been considered and tested for the solution of Maxwell equations. Actually explicit methodology is more rapidly reaching convergence than the implicit one.

Test validation results are in good agreement with available experimental data, and shows that effect of magneto-fluid-dynamic interaction could be relevant for aerospace applications.

Currently CIRA and ALTA with Bologna University Electrical Department are working together in the design of MFD interaction experiments in air in order to understand the real opportunity of application of these technology to Earth re-entry (Cristofolini et al 2010).

A Magneto-Fluid-Dynamic Model and

BC Boundary Condition

EU Inviscid computation FDM Finite Difference Method FVM Finite Volume Method FDS Flux Difference Splitting

MFD Magneto Fluid Dynamics MHD Magneto Hydro Dynamics NE Non equilibrium computation NS Navier Stokes computation

PIC Particle In Cell

**8. References** 

CFD Computational Fluid Dynamics

Stokes solver ENO Essentially Not Oscillatory

CFL Courant - Freidricks -Lewy Number

CMFD Computational Magneto Fluid Dynamics

H3NS Hypersonic high enthalpy 3D Navier-Stokes solver.

PG Perfect Gas approximation computation

Thermophysics Conference Miami AIAA 2007-4048.

ce electron oscillation

e electronic i ionic r radiative v vibrational stag stagnation point

tot total w wall

**7. Acronyms** 

Computational Solving Methodologies for Aerospace Applications 147

EMC3NS Electro Magnetic Generalized Chemistry High Enthalpy 3D Navier-

Battista F. (2009). *Magneto-Fluid-Dynamics Methods for Hypersonic Aerothermodynamics and Space Propulsion Applications*, Ph.D. thesis, Pisa University, etd-11302009-122005. Battista, F., Clemente, M. D., and Rufolo, G. C. (2007). *Aerothermal Environment Definition for* 

Baum, C. E. (1965). *The calculation of Conduction Electron Parameters*, Tech. rep., AFWL.

*a Reusable Experimental Re-entry Vehicle Wing Leading edge*, 46th AIAA

## **6. Nomenclature**

#### **Latins**



#### **Subscripts**



### **7. Acronyms**

146 Fluid Dynamics, Computational Modeling and Applications

**6. Nomenclature** 

e- electron E energy **E** Electric field

m mass

P Pressure *q* Heat flux

R Radius T Temperature

t Time S Surface

*V* Velocity V Volume x x-coordinate y y-coordinate z z-coordinate

*<sup>D</sup>* Debye Length

μ Mobility

ρ Density σ Conductivity

*<sup>e</sup>* Dielectric tensor

Production terms

Stress tensor

Frequency

ad adiabatic wall

a activation

Collision Frequency

freestream conditions

**Greeks** 

0 

**Subscripts** 

**B** Magnetic field D Diffusion coefficients e energy, electron charge

k Rate constant J Flux, Currents

*n* Number density NA Avogadro Number

h Enthalpy per unit mass

q generic system variable

ui velocity components

Dielectric Vacuum constant

M Mach number, Molecular Weight

**Latins** 


### **8. References**


A Magneto-Fluid-Dynamic Model and

101528.

Lectures notes.

John Wiley and Sons.

Vol. 154, pp. 204–309.

Conference, 2007.

New York.

Conference.

Computational Solving Methodologies for Aerospace Applications 149

Gordon, S., McBride, P. (1971). *Computer Program for Calculation of Complex Chemical* 

Gupta, R., Yos, J., and Thompson, R. A. (1989). *A Review of Reaction Rates and Thermodynamic* 

Hansen F., Heims P.S.(1968). *A review of the thermodynamic, transport and chemical reaction rate* 

Helander P. (2006). *Magnetohydrodynamics* , 43° Culham Plasma Physic Summer School,

Hirshfielder, J. O., Curtiss, C. F., and Bird, R. B. (1954). *Molecular Theory of Gases and Liquids*,

Kang, S.-W., Jones, W. L., and Dunn, M. G. (1973) *Theoretical and MeasuredElectron-Density* 

Kenneth, G. P.and Philip, L. R., Timur, J. L., Tamas, I. G., and Darren,L. D. Z. (1998). *A* 

MacCormack, R. W. (2007). *Numerical Simulation of Aerodynamic Flow within a Strong Magnetic* 

Mitchner, M. and Kruger, C. H. (1973), *Partially Ionized Gases*, chap. IV, John Wiley & Sons,

Miura, K. and Groth, C. (2007). *Development of Two-Fluid Magnetohydrodynamics Model for* 

Orlanski, I. (1976). *A Simple Boundary Condition for Unbounded Hyperbolic Flows*, Journal of

Park, W., Belova, E., Fu, G., and Tang, X. (1999). *Plasma simulation studies using multilevel* 

Park,C. (1993). *Review of Chemical-Kinetic Problems of future NASA Missions, I: Earth Entries*,

Resler, E. and Sears, W. (1958). *The prospects for magneto-aerodynamics*, Journal of the

Schettino, A., Battista, F., Ranuzzi, G., and D'Ambrosio, D. (2008). *Rebuilding of new* 

Schexnayder C.J., Evans J.S. and Huber Paul W. (1977), *Comparison of theoretical and* 

Taku, A. and Atsushi, F. (2004). *Full Wave Analysis of MHD Modes Using Multi-Fluid Dielectric Tensor in an Inhomogeneous Plasma*, J. Plasma Fusion Res., Vol. 6, pp. 164–168.

*experimental tests on a double cone at Mach 9*, The 6th European Symposium on

*experimental electron density for RAM C flights from The Entry Plasma Sheath and its Effect on Space Vehicle Electromagnetic Systems*. NASA SP-252, 622 pages, published

Computational Physic, Vol. 21, 1976, pp. 251–269.

Aeronautical Sciences, Vol. 25, pp. 235–245.

Aerothermodynamics for Space Vehicles.

by NASA, Washington, D.C., p.277.

*physics models*, Physics of Plasma, Vol. 6, pp.1796–1803.

Journal of Thermophysics and Heat Transfer, Vol.7,No.3.

*Solution-Adaptive Upwind Scheme for Ideal Magnetohydrodynamics*, Physics of Plasma,

*Field with Hall Current and Ion Slip*, 38th AIAA Plasmadynamics and Lasers

*Non-Equilibrium Anisotropic Plasma Flows*, 38th AIAA Plasmadynamics and Lasers

*Chapman-Jouguet Detonations*, NASA Report SP-273, 1971.

*properties of high temperature air*, NACA Technical note 4359.

*Distributions at High Altitudes*, AIAA Journal, Vol. 11, pp. N.2.

*Equilibrium Compositions, Rocket performance, Incident and Reflected Shocks and* 

*and Transport Properties for an 11-Species Air Model for Chemical and Thermal Nonequilibrium Calculations to 30000 K*, Tech. rep., NASA Technical Memorandum


Bird R.B., Stewart W.E., and Lightfoot E.N. (1954). *Transport Phenomena*, John Wiley and

Borghi C.A., Cristofolini A., Carraro M.R., Gorse C., Colonna G., Passaro A., Paganucci F.

Borrelli, S. and Pandolfi, M. (1990). *An Upwind Formulation for the Numerical Prediction of Non* 

Bush, W. (1958). *Magnetohydrodynamic-hypersonic flow past a blunt body*, Journal of the

Capitelli, M., Gorse, C., Longo, S., and Giordano., D. (2000). *Collision integrals of high-*

Chung, T.J. (2010). *Computational Fluid Dynamics* , Cambridge University Press, Second

Cristofolini A. et al. (2006), *Experimental Investigation on the MHD Interaction around a Sharp* 

Cristofolini A., Borghi C.A., Neretti G., Passaro A., Baccarella D., Schettino A., Battista F.,

*a Blunt Body,* 41st Plasmadynamics and Lasers Conference, AIAA 2010-4490. Cristofolini, A., Borghi, C. A., and Colonna, G. (2007). *Numerical Analysis of the Experimental* 

D'Ambrosio, D. and Giordano, D. (2005). *A Numerical Method for Two-Dimensional Hypersonic* 

Dandy, R. (1993). *Plasma Physic*, Springer-Praxis Series in Space Science and Technology,

Drellishak K. S., Knopp C.F., Bulent Cambel A. (1963), *Partition Function and Thermodynamic* 

Dunn, M.G., Lordi, J.A. (1969). *Measurement of Electron Temperature and Number Density in* 

Giordano D. (2002). *Hypersonic Flow Governing Equation with electromagnetic field*, AIAA 2002-

Giordano, D. and D'Ambrosio, D.(2004). *Electromagnetic Fluid Dynamics for Aerospace* 

Gokcet T.(2004), *N2-CH4-Ar chemical kinetic model for simulation of Atmospheric Entry to Titan*

*properties of Argon Plasma*, The Physics of fluids Volume 6, N.9.

and Technologies Conference, AIAA 2006-8050.

in Fluid Dynamics, Oxford, UK.

AeroSpace Sciences, Vol. 25, pp. 685–690.

Lasers Conference, San Francisco, California.

(2006). *Non intrusive Characterization of the Ionized Flow Produced by Nozzle of an Hypersonic Wind Tunnel*, 14th AIAA/AHI Space Planes and Hypersonic Systems

*Equilibrium Hypersonic Flows*, 12th International Conference on Numerical Methods

*temperature air species,* Journal of Thermophysics and Heat Transfer , Vol. 14, pp.

*Cone in an Ionized Argon Flow*, AIAA-2006-3075, 37th AIAA Plasmadynamics and

(2010). *Experimental Activities on the MHD Interaction in a Hypersonic Air Flow around* 

*Results of the MHD Interaction around a Sharp Cone*," 38th AIAA Plasmadynamics

*Fully Coupled Electromagnetic Fluid Dynamics*, 36th AIAA Plasmadynamics and

*Shock Tunnel Flows. Part II: Dissociative Recombination rate in air*, AIAA Journal Vol.7,

*applications. Part I*, 35rd AIAA Plasmadynamics and Lasers Conference, AIAA

Sons, New York.

259–268.

edition.

and Lasers Conference.

N.11,2099-2104.

Paper 2004-2165.

AIAA 2004-2469.

2165.

Lasers Conference AIAA 2005- 5374.

Cambridge Press, Cambridge.


**7** 

Habib Alehossein

*Australia* 

**Mechanics of Multi-Phase Frictional** 

**Visco-Plastic, Non-Newtonian, Depositing** 

**Fluid Flow in Pipes, Disks and Channels** 

*CSIRO Earth Science and Resource Engineering, University of Queensland* 

*Background***:** Fresh concrete, fly ash and mining slurries are all frictional-visco-plastic fluids. Fresh concrete flow in Tremie pipes is used to control concrete flow rate and minimise bleeding and dilution when concrete is poured into deep submerged excavations for pile foundation construction. Slurries with very fine aggregates are used to backfill underground voids and mines to prevent subsidence and surface structural damage. Backfilling and injection of granular materials into mining induced voids, separated beddings and cracks, as either diluted slurry or concrete paste, is widely used to control subsidence. As a viable environmental solution, mine waste and rejected materials from underground coal seams are used in both backfilling and injection mine operations. For example, during longwall mining the grout slurry is pumped into the separated beds of the fractured rock mass through a pipeline connected to a central vertical borehole, which is drilled deep into the inter-burden rock strata above the coal mine seam. Either as dilute slurry or thick paste or cake, the fill material normally needs to travel a significant longitudinal distance either in a channel, a tremie pipe, a long pipeline, or radially in a disk shaped crack in the rock mass. An undesirable blockage can occur in the central borehole, in the crack or in the transportation channel or pipeline system when the slurry velocity falls below a certain critical threshold velocity, indicating a material phase change from cohesive-viscous to cohesive-frictional. This chapter presents complete analytical solutions of the required pump pressure versus fluid volume rate for such multi-phase fluids, which are categorised as frictional Bingham-Herschel-Bulkley fluids. The theory derived can be applied to flow of such fluids in pipes, disks and open channels. Furthermore, general analytical solutions have been developed for such fluids in terms of velocity and pressure gradients and velocity and pressure, as a function of flow length (e.g. pipe length, disk radial distance, or channel length) from which special and familiar equations for simpler fluids are derivable by mathematical reduction of the general equations. The formulation is distinct in considering many new aspects including: variable shear parameters rather than fixed values; inclusion of total nonlinear behaviour; and, implementation of a friction function to mimic behaviour of the depositing and consolidating stiff slurry or paste, which can cause a significant

*Bingham-Herschel-Bulkley fluids*: Recent laboratory and field experiments on mine-backfill fluids, slurries, cements, pastes and concretes proved their wide range of shear resistance

pressure rise as a result of the increased shear resistance.

**1. Introduction** 

