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

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

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: 1. the velocity is not relativistic uc  ; 2. the phenomena under consideration are slow enough that The equations are written in the following; mass balance equations of the components in chemical non equilibrium conditions yield, where i m J 1 are the components mass diffusive fluxes and the i  2 are the production rate of the i-th specie, Fluid Dynamics, Computational Modeling and Applications 124 Pointing out the charge density as (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 where the current is the sum of convection and conduction current is the viscous stress tensor and the total energy equation reads: where the total energy is 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).

Chemical model
Given a generic set of N r chemical reactions, '' ' A Magneto-Fluid-Dynamic Model and Computational Solving Methodologies for Aerospace Applications 125 where ' ir and '' ir are the stoichiometric coefficients for the r-th reaction and  i are the chemical species involved in the reactions, the production rate Ω i of species i can be written as where k fr and k br are the forward and backward rate constants of the reaction r, which are assumed to follow the Arrhenius temperature law: where the pre-exponential factors A i , the temperature exponents  r , and the activation energies E i depend upon the adopted kinetic scheme.
The backward rate constants are related to the forward rate constants through the equilibrium constants K eq,r , i.e.: 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:

 
this formula is consequent from the expression of binary collision rate where m 12 is the reduced mass, d 12 is the mean diameter, k B 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):

 
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 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.where R and c p 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: The single species specific heats and enthalpies are computed by using the (Gordon and McBride 1971)  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).

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.
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. where is the degree of ionization.For what concerns Q eH an expression for air could be found in (Baum 1965).

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: (2.20)

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: In the previous equations the cyclotronic vector has been considered that divided by  gives the Hall parameter vector that represent the Hall parameter along all the space dimension: Decoupling the equations multiplying by ne, and reminding from previous section that From Eqn. (2.22) the tensor e  is directly derived.

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: so they can be simply discretized following a finite volume approach: The unknown vector   W , the flux vector ( F ) and the source term vector (  ) used in the finite volume discretization are reported below.

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: 1. Fluid-Dynamics equations sub-system 2. Maxwell equations sub-system 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: 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 intensity.


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 (keeping E and B constant).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.

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).


Electromagnetic terms added at the equations implemented  Thermodynamic properties formulation in terms of Gordon-Mc Bride Polynomial fits.(Gordon and McBride 1971)  Extension of thermodynamic properties (specific heat, enthalpy and entropy) to higher temperatures.


Included ionized species transport properties.


Effects of electromagnetic field on transport accounted.Models considered for these upgrades have been widely discussed in previous Sections.

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.: With: The superscript "M" reminds that this method is being applied to the Maxwell equations.

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: Here the subscript "N" addresses the cell we are referring to, while the superscript "K" indicates the time interval considered.

M IN
F and M OUT F can be written as: 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: 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):  This equation is just a system of tot n scalar equations, where tot n ( x y z nnn  ) is the total number of grid cells.It can be written in the following more compact form: Where: 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.

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: Each term is evaluated at instant "k".Rearranging the previous equation for the simple case of cubic cells: www.intechopen.comA Magneto-Fluid-Dynamic Model and 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 carefully in order to have an acceptable time step and to grant the method stability.
Considering the very high speed of signal propagation ( 8 31 0 m / s c  ), 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.

 
1 11 1212 11 32 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).

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).Initially the condition adopted on the fluxes crossing the boundary cells was simply the following one: , 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)  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.

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

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.

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 T 0 =6850 K and P 0 =0.94738MPa.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 n e behaviour.

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 uB  .

 
F q EuB   (4.1) This force exerted upon the charged particles in the shock layer tends to drive them in the direction uB  , 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: The Faraday current density generates the JB  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 uB  directions and weaken Faraday currents as well as JB  force.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.

Conclusions
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).

Latins
Debye length and the ion gyro radius are small D L   , i L   .

Fig. 3 .
Fig. 3. System solving strategySo, given the initial conditions (both for fluid-dynamic variables and electromagnetic fields) the solution procedure to solve the system is: The fluid-dynamic subsystem is solved.Maxwellequations 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 intensity.Thenewly 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 (keeping E and B constant).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.8)  : is used to indicate a spatial interval   : is used to indicate a time interval , x y nn : are the numbers of cells along x and y apxis respectively www.intechopen.comFluid Dynamics, Computational Modeling and Applications 136 [C], [L], [R], etc… are 6x6 square matrixes (remind that F and W   are six-elements vectors).

Fig. 7 .
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.

Fig. 8 .
Fig. 8. Magnetic field B z 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.

Fig. 9 .
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 .
Fig. 10.Numerical and experimental data comparison at different RAM-C stations.

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

Fig. 13 .Fig. 16 .
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 inCristofolini et al. (2006) and numerically investigated.It consists in an array of the magnets that has been assembled to form a conical

Fig. 17 .Fig. 18 .
Fig. 17.Computed pressure field without and with MHD interaction 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 N t is dependent on the fluxes in the nearby cells at instants 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 a Reusable Experimental Re-entry Vehicle Wing Leading edge, 46th AIAA Thermophysics Conference Miami AIAA 2007-4048.Baum, C. E. (1965).The calculation of Conduction Electron Parameters, Tech.rep., AFWL.
DDebye Length e 