MPSalsa A Finite Element Computer Program for Reacting Flow Problems Part 2 - User\'s Guide

Share Embed


Descrição do Produto

SAND95–2752 May, 1996

Distribution Category UC-405

MPSalsa A FINITE ELEMENT COMPUTER PROGRAM FOR REACTING FLOW PROBLEMS PART 1 - THEORETICAL DEVELOPMENT

John N. Shadid1, Harry K. Moffat2, Scott A. Hutchinson1, Gary L. Hennigan1, Karen D. Devine3, Andrew G. Salinger1 1Parallel

Computational Sciences Department

2Chemical 3Parallel

Processing Sciences Department

Computing Sciences Department

Sandia National Laboratories Albuquerque, New Mexico 87185

Abstract The theoretical background for the finite element computer program, MPSalsa, is presented in detail. MPSalsa is designed to solve laminar, low Mach number, two- or three-dimensional incompressible and variable density reacting fluid flows on massively parallel computers, using a Petrov-Galerkin finite element formulation. The code has the capability to solve coupled fluid flow, heat transport, multicomponent species transport, and finite-rate chemical reactions, and to solve coupled multiple Poisson or advection-diffusion-reaction equations. The program employs the CHEMKIN library to provide a rigorous treatment of multicomponent ideal gas kinetics and transport. Chemical reactions occurring in the gas phase and on surfaces are treated by calls to CHEMKIN and SURFACE CHEMKIN, respectively. The code employs unstructured meshes, using the EXODUS II finite element database suite of programs for its input and output files. MPSalsa solves both transient and steady flows by using fully implicit time integration, an inexact Newton method and iterative solvers based on preconditioned Krylov methods as implemented in the Aztec solver library.

1

Acknowledgments The authors wish to thank D. K. Gartling and M. L. Martinez for their helpful discussions regarding various aspects of the finite element formulation used in MPSalsa and M. A. Christon for helpful comments on the draft of this document. In addition, the work of H. Walker on an initial implementation of the inexact Newton method is gratefully acknowledged. This work was partially supported by the U. S. Department of Energy Office of Scientific Computing under contract DEAC04-94AL85000, and by a Sandia Laboratory Directed Research and Development proposal.

2

Table of Contents: 1. Introduction .......................................................................................................................5 2. Problem Types and Equations of State ..............................................................................8 2.1 Problem Types .....................................................................................................8 2.2 Material Properties ..............................................................................................8 2.3 Units Within the Program ..................................................................................12 2.4 Exact Solutions ..................................................................................................12 3. Governing Transport-Reaction Equations .......................................................................13 3.1 Momentum Transport Equation ........................................................................13 3.2 Total Mass Conservation Equation ...................................................................14 3.3 Energy Transport Equation ................................................................................14 3.3.1 Temperature Formulation of the Energy Equation .............................14 3.3.2 Enthalpy Formulation of the Energy Equation ...................................16 3.4 Species Mass Transport Equation .....................................................................16 3.4.1 Diffusion Velocities ............................................................................18 3.5 Calculation of Diffusion Velocities ...................................................................19 3.6 Implementation of Gas Phase Reactions ...........................................................21 3.7 Implementation of Surface Phase Reactions .....................................................23 3.8 Summary of Transport Equations Implemented ................................................27 4. Boundary Conditions .......................................................................................................28 4.1 Specification of Boundary Conditions ..............................................................28 4.2 Momentum Equations .......................................................................................28 4.3 Total Continuity Equation .................................................................................30 4.4 Internal Energy Continuity Equation .................................................................30 4.5 Gas-Phase Species Continuity Equations ..........................................................31 5. Finite Element Approximation of the Transport Equations ............................................33 5.1 Interpolation Functions and Quadrature Rules ..................................................33 5.2 Evaluation of Surface Integrals .........................................................................34 5.3 Matrix Equations ...............................................................................................34

3

6. Solution Procedures .........................................................................................................38 6.1 Implementation on Multiple Processors ............................................................38 6.2 Numerical Properties .........................................................................................41 6.3 Transient Solution Algorithms ..........................................................................41 6.3.1 Forward /Backward Euler Integration ................................................42 6.3.2 Adams-Bashforth/Trapezoidal Rule Integration ................................42 6.3.3 Time Integration Procedures ..............................................................43 6.3.4 Time Step Control ..............................................................................44 6.4 Inexact Newton Method with Backtracking ......................................................45 6.4.1 Nonlinear Convergence Criteria .........................................................45 6.5 Linear System Solvers .......................................................................................45 Appendix A: Semi-discrete Conservation Form of Transport Equations ...........................47 Appendix B: Petrov-Galerkin Pressure Stabilization Constant ...........................................49 Appendix C: Derivation of the Discrete Matrix Equations .................................................50 C.1 Obtaining the Newton-Kantoravich Equations .................................................50 C.2 Residual Equations ...........................................................................................51 C.3 The Continuum Newton-Kantoravich Equations .............................................52 C.3.1 Momentum Equation ..........................................................................52 C.3.2 Mixture Total Continuity Equation ....................................................54 C.3.3 Thermal Energy Transport Equation ..................................................55 C.3.4 Species Transport Equations ..............................................................56 C.4 Discrete Newton-Kantoravich Equations .........................................................57 C.5 Definition of Jacobian Block Matrices .............................................................63

4

1. Introduction

The theoretical development and numerical procedures for the finite element computer program, MPSalsa, are presented in detail in this document. A companion user’s manual provides details on using MPSalsa for specific applications along with a number of example problems [1]. Employing unstructured meshes on massively parallel (MP) computers, MPSalsa is designed to solve two- or three-dimensional problems which exhibit coupled fluid flow, heat transport, species transport, and chemical reactions. The modeling equations defined in MPSalsa for fluid flow and mass conservation are the momentum transport and the total mass continuity equation for incompressible or variable density Newtonian fluids (Navier-Stokes equations). The heat transport equation and an arbitrary number of species transport-reaction equations couple strongly with each other through chemical reaction source terms and with the fluid flow equations through property variation and body force terms. The program uses the CHEMKIN suite of library routines to provide a rigorous treatment of ideal-gas multicomponent transport, including the effects of thermal diffusion [2]. The mixture-averaged diffusion approximation is available in addition to the computationally-expensive DixonLewis formulation. Chemical reactions occurring in the gas phase and on surfaces are also treated by calls to CHEMKIN [3] and SURFACE CHEMKIN [4], respectively. Because of this, MPSalsa can handle varying numbers and types of chemical reactions and species in a robust manner. For example, the code can handle the complex temperature and pressure dependence predicted for unimolecular reactions (using the Troe parameterization), important for chemical vapor deposition (CVD) systems, which typically run at sub-atmospheric pressures. Surface site fractions and bulkphase mole fractions are defined on all reacting surfaces using the SURFACE CHEMKIN package. Through this method, complex Langmuir-Hinshelwood-type and precursor adsorption surface mechanisms, characteristic of many real CVD and catalysis surface systems, can be incorporated into the reacting flow analysis code. The capability of modeling simple dilute species transport and reaction, without the need of linking to CHEMKIN, is also included in MPSalsa. The user can extend the models past what has been pre-defined within MPSalsa [1]. Functions can be written to represent additional source terms, special boundary conditions, and variations in physical properties, any of which can be dependent on the current solution, position, or time. The discretization method is a Petrov-Galerkin finite element method (PGFEM) with pressure stabilization. Both steady and transient flows may be analyzed. The time integration methods include true transient, pseudo-transient, and steady implicit solvers. The overall solution is obtained by fully-coupled, implicit, parallel iterative solvers based on preconditioned nonsymmetric Krylov subspace methods. Presently, MPSalsa can simulate low Mach number (< 0.3) flows, where an algorithm employing an implicit coupling between the pressure and velocity field is required. MPSalsa employs unstructured grids, using the EXODUS II finite element database suite of programs for its input and output files [5, 6, 7]. Therefore, it can be used in conjunction with the CUBIT mesh generation package [5], as well as other mesh generation packages that support the EXODUS II standard. A number of pre- and post-processing routines for the EXODUS II database can be used. Currently, two- and three-dimensional grids with Cartesian coordinates are supported.

5

MPSalsa includes both first- and second-order predictor-corrector time integration schemes; these methods use explicit predictors and fully implicit corrector methods based on forward/backward Euler and Adams Bashforth/Trapezoidal rule methods, respectively. At each time step, a prediction of the solution and its time derivative are generated from the appropriate time integration scheme. This prediction is used as the initial guess for the fully coupled non-linear problem generated at each time step. The non-linear problem is solved using an inexact Newton method. At each step of the non-linear problem, a “residual vector” and a “Jacobian matrix” are generated, based on the current solution approximation. The resulting linear problem is solved using iterative methods based on preconditioned Krylov-subspace techniques. The accuracy or convergence criteria for solving the linear subproblem is controlled by the inexact Newton algorithm. This algorithm selects the convergence criteria based on how well the linear subproblems are approximating the underlying nonlinear problem. As is the case with most adaptive ODE integration codes, the accuracy to which the non-linear problem is solved is based on a time-step truncation error estimate. The adaptive time integration method uses a user-specified error tolerance and a time-truncation error estimate from the compatible-order predictor/corrector methods to automatically select time step sizes to control time step truncation error at a user-specified tolerance. From its inception, MPSalsa has been designed for distributed memory MIMD computers with hundreds to thousands of processors. It also runs on traditional serial workstations and networks of serial workstations. Interprocessor data communication and global synchronization are accomplished by a small number of message passing routines. These routines have been ported to many different message passing protocols, including the MPI standard and the native nCUBE and Intel Paragon protocols. To achieve efficient parallel execution, the unstructured finite element mesh is partitioned or load-balanced in a preprocessing step. Here, each processor is assigned nodes from the mesh such that the computational load is balanced and the total amount of information communicated between neighboring processors is minimized. Each processor is then responsible for calculating updates for all the unknowns at each of its assigned FE nodes. Each processor also stores and performs operations on the rows in the fully-summed, distributed matrix associated with these unknowns. Along processor subdomain boundaries, replicated FE unknowns, called “ghost unknowns,” are stored and updated through interprocessor communication. These ghost unknowns, assigned to neighboring processors, are quantities needed for the local residual calculation and matrix-vector multiplication on a processor. Interprocessor communication occurs for each step of the iterative solution of the linear system as well as for each outer step in the non-linear and time-transient algorithms. This communication constitutes the major unstructured interprocessor communication cost in the program, and its algorithm has been extensively optimized within MPSalsa [8]. Solution output from the program is achieved through several means. Output can be written to either a standard serial EXODUS file format [6, 7] or a “parallel extension” of the EXODUS file format [9]. This extension consists of writing an individual standard serial EXODUS file for each processor with an extra array that maps the local node numbering scheme on an individual processor to the global node numbering scheme. The format can be used on both MP computers, such as the Intel Paragon, and distributed computing systems, such as groups of workstations. This parallel I/O capability can be used with today’s primitive parallel I/O facilities with nearly linear speedup. This report serves as an introduction to MPSalsa. A companion user’s manual contains a detailed description of the input and solution options, as well as several example problems that have been solved by MPSalsa [1]. The target problem classes of MPSalsa are discussed in Section 2,

6

along with the currently supported material types and equations of state. Section 3 introduces the governing transport-reaction equations. Special sections on the calculation of the multicomponent diffusional fluxes and gas-phase reactions are included as well. The treatment of surface species and surface reaction source terms is also discussed. Subsection 3.8 contains a summary of the bulk transport equations solved within the code. Section 4 contains a general discussion of the implementation of boundary conditions within MPSalsa where boundary conditions specific to each equation are introduced. In Section 5, the finite element implementation of the transport-reaction equations, the supported interpolation functions, quadrature rules, and methodology for calculating surface integrals are introduced. The matrix equations are also presented to display the essential form of the system of coupled equations. Terms included and excluded from the Jacobian matrix are delineated in Appendix C. Section 6 contains the solution methodology at the algorithm level. The parallel implementation of the code is described, and the nonlinear solver and the linear system solvers, along with their respective convergence criteria, are discussed. The algorithmic details of the Aztec library of Krylov solvers and preconditioners are left to companion documents [10, 11].

7

2. Problem Types and Equations of State

2.1

Problem Types

MPSalsa is designed to solve the governing transport-reaction equations for momentum, total mass, thermal energy, and species. In addition, MPSalsa allows the user to solve a reasonably general set of coupled transport-reaction equations by specification of general transport coefficients and source terms. The scope of the problem types that a program can handle is determined, in part, by the discretization scheme and solution method. MPSalsa employs a highly coupled approach to the solution of its equation set, by storing all cross terms in the Jacobian. The fully-summed distributed Jacobian is stored so that highly effective general algebraic preconditioners such as ILU with partial fill-in and block ILU factorizations may be used to reduce the total number of iterations in the linear solver. Thus, MPSalsa is most effective on highly coupled problems that require an implicit solution technique. It is less efficient on problems that can be solved with explicit or semiimplicit solution techniques, such as high Mach number flows or weakly coupled systems. MPSalsa is currently designed to solve low Reynolds number laminar flow problems, and no stabilization terms have yet been added to avoid oscillatory behavior of the solution for high cell-Reynolds numbers. Additionally, the filtering of the density by eliminating the hydrodynamic pressure dependence limits the problem classes MPSalsa can currently handle to low Mach number flows. However, within these bounds, the transport-reaction systems and geometric complexity that MPSalsa can handle are quite general. The determination of which equations are solved, as well as which operators are included, is done by specifying the “problem type.” This also determines what types of unknowns are included in the solution vector. Table 2-1 shows the available options for the problem type. As the table points out, diffusion operators are always included, while inclusion of the convection operator depends on the particular problem type. Single, general PDEs that don’t fall into any of the categories in Table 2-1 may be handled either with the energy equation/temperature unknown or the species conservation/mass fraction unknown. Systems of general PDEs are handled with the mass species transport equations and can optionally be coupled to the momentum, thermal energy and total mass equations. Each problem type has a default setting for whether the equations are linear or non-linear. MPSalsa contains logic for the efficient handling of both cases. The default linearity setting can be overridden as well. For heterogeneous or multi-physics problem types, different domains with different material types, such as a solid and an ideal gas, are used. A varying number of transport equations are then solved on each domain. While this type of problem has not been fully implemented in MPSalsa, the underlying data structures are in place. In particular, the matrix storage format, Variable Block Row (VBR) sparse matrix format [12], allows for a different number of equations to be solved for per node.

2.2

Material Properties

The assignment of material properties starts with designating each region, specified by EXODUS II element blocks, with a “material model.” Material models are broadly classified within

8

Table 2-1: Problem Types

Diffusion

energy_diff

x

x

energy_conv_diff

x

x

mass_diff

Convection

Operators Included

Mass/ Pressure

Momentum/ Velocity

Species/ Mass Fraction

Problem Types

Energy/ Temperature

Transport Equation/ Unknowns

x

x

energy_mass_diff

x

x

x

energy_mass_conv_diff

x

x

x

x

x

stokes_flow

x

x

x

fluid_flow

x

x

x

x

x

x

x

x

x

x

x

x

x

fluid_flow_energy

x

fluid_flow_mass whole_banana

x

x

x

x

x

x

advection_diff

input

x

input

input

x

x

9

MPSalsa as belonging to a “material type” which are listed in Table 2-2. The material type is used Table 2-2: Material Types Material Type

Description

CHEMKIN

Ideal Gas - Use the ideal gas mixture equation of state, and calculate transport properties and reaction rates via CHEMKIN.

NEWTONIAN

Newtonian fluid, i.e., has a Newtonian stress tensor formulation. The default is to use constant fluid and transport properties.

BOUSSINESQ

Boussinesq fluid, i.e., a Newtonian fluid with a constant thermal expansion coefficient. Density varies only in the body force term.

SOLID

Bulk solid with isotropic transport properties

NNEWTONIAN

Non-newtonian fluid (not yet implemented)

ANISOTROPIC_SOLID

Material that has an anisotropic thermal conductivity and species diffusivities (not yet fully implemented).

extensively within the code for conditional evaluation of equations of state, transport property computations and source terms. When a CHEMKIN material type is defined in a problem, MPSalsa reads the CHEMKIN binary work arrays produced by CHEMKIN preprocessors. Details of this process can be found in the MPSalsa User’s Guide [1]. From these work arrays, MPSalsa obtains the number of gas-phase species, the number of surface phases and surface-phase site fractions, and the number of bulk mole fractions. All gas-phase transport properties are obtained from the TRANLIB library [24], which evaluates gas-phase multicomponent transport properties. The ideal gas equation of state given by Eqn. 1 is used to yield expressions for the density, ρ . Ng

P o = ρRT

Yj

∑ -----Wj

ρRT = ----------------------N

j=1

(1)

g

∑ W jX j

j=1

N g is the number of gas phase species, Y j is the mass fraction of the jth species, X j is the mole fraction of the jth species, and W j is the molecular weight of the jth species. P o is the thermodynamic pressure. The CHEMKIN material type assigns a “special species label” to one of the species. The conservation equation for that species is replaced by the condition that the sum of the mass fractions must equal one: Ng

∑ Yk k=1

10

= 1.

(2)

The caloric equation of state for an ideal gas mixture is used for CHEMKIN materials. In this model, h, the specific enthalpy of the mixture, does not depend on the total pressure. Eqn. 3 provides the expression for the specific enthalpy in terms of the partial specific enthalpies for each species and the mass fractions. Since an ideal solution is assumed, the partial specific enthalpies are equal to the pure specific enthalpies of each species in its reference state. T

Ng

h =



o hˆ j ( T ) = W j ∆H f , j(T o) +

hˆ j ( T )Y j

j=1

∫ Cˆ p, j dT

(3)

To

o

∆H f , j(T o) is the heat of formation of the jth species in its standard state and at the common reference temperature (which for CHEMKIN is T o = 298.15K ). The standard state for gases corresponds to an ideal, pure gas state at 1 atm. Thermodynamic information for the CHEMKIN material type is obtained from the CHEMKIN thermodynamics data base or the CHEMKIN input file. The calculations in Eqn. 3 are carried out within CHEMKIN. Cˆ p, j , the specific heat at constant pressure for species j, is a polynomial function of temperature. In the NEWTONIAN material type, all transport properties, as well as the density, are assumed constant. This assumption can be overridden by specification of variable properties for a number of the transport properties. In the BOUSSINESQ material type, the default is for the density to be constant in all equations, except for the body force term in the momentum equations. In this term, the density is assumed to be a linear function of the temperature. The density can be expressed in terms of the coefficient of volumetric expansion, β. ρ = ρ(T 0) [ 1 – β ( T – T 0 ) ]

,

where

1 ∂ρ β = – --- -----ρ ∂T

(4) T = T0 p

Note that for an ideal gas, β = 1 ⁄ T , and, thus, it is not a constant. For the BOUSSINESQ material type, β is supplied by the user. The SOLID material type is a placeholder set aside for the future anticipated capability to do conjugate heat transfer problems in domains with both solid and fluid regions. These problems have regions where the momentum equations are not solved. Currently, this capability is not available in MPSalsa. In MPSalsa, both constant and variable thermal transport properties can be used. The NNEWTONIAN material type is defined for the specification of non-Newtonian constituitive equations (as well as the required additional Jacobian entries) for viscosity. The NEWTONIAN, BOUSSINESQ, or SOLID material types can be used if species equations are desired but the CHEMKIN subroutine library for mixtures of ideal gases is not to be used. The default for these non-CHEMKIN materials is to NOT enforce Eqn. 2. However, this default can be overridden. The lack of Eqn. 2 represents the situation where all species transport equations represent only dilute components of phases. The majority component of a phase is not represented by a species equation. For all equation types, there is a capability in MPSalsa for including both volumetric and surface source terms in the residuals and, just as importantly for stiff terms, their Jacobian contributions in the matrix used to relax the equations. Volumetric source terms are specified as part of the materials model using either built-in or user-specified functions. In contrast, surface source terms

11

are specified as surface boundary conditions. They are applied by integrating over surfaces defined in the finite element model. These boundary conditions can also be user-specified functions or built-in functions representing well-known cases, such as those that correspond to convective or radiative heat transfer and sticking coefficient reactions. For boundary conditions at surfaces where deposition or etching of bulk phases occurs, SURFACE CHEMKIN is used to describe the process’ kinetics and yield values for surface fluxes of gas-phase species. The capability for solving Stefan flow problems, i.e., problems that have a net normal mass flux at the surface that depends on the surface reaction rate, is built into this “reacting surface” boundary condition.

2.3

Units Within the Program

Non-dimensionalization of the equations is not done within MPSalsa. Except when CHEMKIN is used, no units are a priori specified within the program. CHEMKIN produces quantities such as transport properties, densities, pressure, energy, and species rates of production in terms of the CGS units system, i.e., gm, cm, sec, mole, and Kelvin. Therefore, whenever the CHEMKIN material type is used, the user inputs to the program --including boundary condition values -- should also be in CGS units. The specification of the thermodynamic pressure is in atmospheres and the default units for activation energies for gas and surface reaction rates are in cal mole-1 for the CHEMKIN material type. When a material type other than CHEMKIN is being used, the user must specify a constant set of units. Understanding the behavior of a system as a function of non-dimensional numbers, such as the Reynolds number or Grashof number, is a powerful tool. However, this must be carried out by the user indirectly. One way is through use of the continuation routine, where the user can often associate the continuation parameter with a dimensionless group. Another way, which can be seen by comparing the dimensional and non-dimensional formulations of the equations and boundary conditions, is to choose the physical properties such that a single property will represent a dimensionless group; e.g., by setting all other properties to one and using the appropriate domain size and boundary conditions, the gravity unknown will be equivalent to the Rayleigh number. An example of a non-dimensionalization of the equations is provided in the MPSalsa user’s manual [1].

2.4

Exact Solutions

MPSalsa is a large code. The use of test problems with known, exact solutions was found to be essential in verifying the code. Much of the code can be checked by comparing numerically derived solutions against exact solutions, and analyzing mesh convergence of numerical solutions. This includes all of the parallelization aspects of the code as well as the implementation of the EXODUS finite element database on multiple processors. For instance, an exact solution to the timedependent Navier-Stokes equations has been implemented[15]. There are, however, cases where exact solutions are not available to check the validity of the code. Real gases with complicated transport properties are one instance. For these situations, the code was checked against other numerical codes. Two such case studies are included in the user’s manual. One case is a comparison of a rotating disk CVD problem to the 1-D numerical code SPIN [13]; the other case is a comparison of a homogeneous, isotropic gas-phase pyrolysis study to the 0-D code SENKIN [14].

12

3. Governing Transport-Reaction Equations

The equations solved by MPSalsa are based on the governing transport equations for total mass, momentum, energy, and individual gas-phase species. Constitutive relations for the momentum, heat, and species fluxes are based on one of three models: (a) the non-equilibrium statistical mechanical theory of multicomponent, dilute polyatomic gases [16, 17, 18, 19, 20]; (b) a constant property, Boussinesq fluid model; and (c) constituitive equations supplied and linked in by the user through a set of user subroutines. The Boussinesq fluid approximation is suited to the study of convection in liquids, including liquid metals, while the multicomponent gas model is suitable for a mixture of ideal gases at atmospheric pressures or lower. The governing transport equations listed below are given in “conservative form” rather than “advective form.” In the actual numerical implementation, both the conservative and the nonconservative forms of the equations can be solved. Experience indicates that while greater accuracy is not guaranteed by the conservative formulation, long-time numerical integration stability is enhanced. For this reason, both formulations have been included in the numerical solution procedure, as described in Appendix A. An acoustically-filtered formulation of the momentum and mass conservation equations is used within MPSalsa [22, 23]. Thus, a distinction between the hydrodynamic and thermodynamic pressure values is employed in the equation set. Variations in the hydrodynamic pressure, which are assumed small compared to the thermodynamic pressure, are not included in the calculation of the density that appears in the conservation of mass, species, and momentum equations. This assumption has been shown to be valid for Mach numbers lower than 0.3 [22] and has the benefit of filtering out shock formation.

3.1

Momentum Transport Equation

The conservation of momentum is expressed by Eqn. 5 and 6. Assuming a Newtonian stress constituitive equation, there are as many scalar components of the momentum equation as there are spatial dimensions in the problem. ∂ ( ρu ) --------------- + ∇•( ρuu ) – ∇•T – ∂t where

Ng

∑ ρk gk

= 0

(5)

k=1

2 T T = – PI + ϒ = – PI – --- µ ( ∇•u )I + µ [ ∇u + ∇u ] 3

(6)

Here, T is the stress tensor for a Newtonian fluid, I is the unity tensor, ϒ is the viscous stress tensor, and P is the isotropic hydrodynamic pressure. In the pressure-filtered formulation, there is a distinction between the hydrodynamic pressure (used in the transport equations) and the pressure level used in the equation of state, P o . This distinction allows the nearly constant thermodynamic pressure level to be set independently of the relatively small pressure fluctuations due to the hydrodynamic flow. Unlike the treatment in Paolucci [22], there is no global equation for the thermo-

13

dynamic pressure, P o , in the equation set. MPSalsa assumes a Newtonian fluid mixture with zero bulk viscosity where µ is the mixture viscosity and is a function of the temperature and fluid composition. For a multicomponent ideal gas mixture it is a complex function of the temperature and 0.7 the species mole fractions with roughly a T dependence on temperature; µ is obtained from a subroutine call to the TRANLIB package [24]. The last term in Eqn. 5 is the body force term where gk is the sum of all body forces acting on species k, and N g is the total number of species. In most cases not involving charged particles and/ or electromagnetic fields, the body force on each species is the same for all species and reduces to the gravitational force, g. In that case, the last term in Eqn. 5 reduces to ρg . Currently, the only body force considered in the code is gravity which is constant for all molecular species. Additional functionality for this term will be application driven.

3.2

Total Mass Conservation Equation The conservation of total mass within MPSalsa is expressed by Eqn. 7. ∂ρ ------ + ∇•( ρu ) = 0 ∂t

(7)

In this equation, ρ is the mass density of the mixture. Two alternate equations of state are allowed for ρ. Either ρ is considered to be a constant (i.e., the incompressible case or the Boussinesq fluid case where ρ is considered to be a constant, except in the body force term), or ρ is calculated from the ideal gas mixture equation of state Eqn. 1. Thus, for an ideal gas, ρ is not a function of the variable hydrodynamic pressure; it is a function of the constant thermodynamic pressure only. Additionally, a user-defined subroutine can be employed to incorporate an alternate equation of state that is dependent on the local temperature and species compositions as well as the thermodynamic pressure.

3.3

Energy Transport Equation

For high speed flows, the conservation equation in the total energy (i.e., the internal energy plus the kinetic energy) form is normally used. This form is particularly useful for inviscid flows. However, the difficulty with this representation is that for flows in which the molecular transport of thermal energy is large, the implicit coupling of the internal energy or enthalpy to the temperature is weak. Given this, for low speed, incompressible flows this equation is generally translated into either the enthalpy or temperature form. These choices work well for the initial class of problems to be addressed by this code - low Mach number CVD problems. In the code, the specific heat/ temperature form is implemented. However, future versions of the code may include the enthalpy form as it is natural for control volume formulations in which a local conservation of energy property can be attained. 3.3.1 Temperature Formulation of the Energy Equation Eqn. 8 is the internal energy equation in terms of the temperature.

14

Ng

DP ( ρT ) Cˆ p ∂-------------- + ∇•( ρuT ) = – ∇•q + φ + Q˙ + ∑ j k • g k + ------Dt ∂t k=1 +

Ng

Ng

k=1

k=1

(8)

∑ hˆ k ( ∇• jk ) – ∑ hˆ k W k ω˙ k

In this equation, Cˆ p is the specific heat of the mixture at constant pressure. The first term on right hand side is the diffusive heat flux, q , given by Eqn. 9. The second term is the volumetric heat source term from viscous dissipation, φ , given by Eqn. 10. The volumetric energy source term Q˙ is specified by a user function, and j is the diffusive flux of the kth species relative to the massk averaged velocity, u . The net change of potential energy from body force terms into heat energy, ∑ jk • gk , is zero for the single body force term implemented so far, gravity, because gk is equal for all k. The total derivative of pressure, DP ⁄ Dt , represents the reversible exchange of mechanical energy into internal energy. The first term on the second line of Eqn. 8 is the production of internal energy due to diffusion, where hˆ k is the partial specific enthalpy of the kth species. The ˙ k as the last term in Eqn. 8 is the volumetric production of heat due to chemical reactions using ω th net production rate of the k species due to homogeneous chemical reaction and W k as the molecular weight of the kth species. Ng

q = – λ ∇T +



Ng

hˆ k j k –

k=1

RT

-------------- D k d k + q r ∑W kXk T

(9)

k=1

2 2 1 T φ = – ( τ : ∇u ) = – --- µ ( ∇•u ) + --- µ ∇u + ∇u 3 2

2

(10)

The first term in Eqn. 9 is the diffusive flux of energy due to heat conduction. λ is the heat conductivity of the mixture. For gases, it is a complicated isotropic function of the temperature and mass fractions. The second term in Eqn. 9 is the diffusive flux of energy due to species diffusion. The third term is the Dufour effect, the diffusive flux of energy due to thermal diffusion. This term is usually very small and is neglected in the implementation of the code. The last term is the flux of energy due to radiative transport, qr. It is almost always ignored when solving the gas-phase energy continuity equation; i.e., the gas is assumed to be transparent to radiant energy. However, it is included here, because it appears naturally in the specification of the boundary conditions for the energy flux on solid surfaces. For a simple thermodynamic material, the heat flux term from the species diffusive flux and the heat source term originating from the divergence of the species diffusive flux term may be combined to yield a single heat source term due to the diffusive flux, Eqn. 11. This modification is incorporated into Eqn. 8 and 9. Ng

– ∇• ∑ hˆ k j k + k=1

Ng

∑ hˆ k ∇• jk k=1

Ng

= –

∑ jk • Cˆ p, k ∇T

(11)

k=1

15

In the initial implementation of the code, some of the terms in Eqn. 8 are not included because of their relatively small contributions. The body-force source term is omitted since the gravity vector, gk, is equal for all k. The viscous dissipation term and reversible change of mechanical energy into internal energy term ( DP ⁄ Dt ) are dropped since they are small for low Mach number applications. Also, the energy flux terms due to species diffusion, as presented in Eqn. 11, have not yet been included but will be in the near future. 3.3.2 Enthalpy Formulation of the Energy Equation Eqn. 12 is the conservation of energy equation expressed in terms of the mixture enthalpy, h. ∂ ( ρh ) DP -------------- + ∇•( ρuh ) = – ∇•q + φ + Q˙ + -------- + ∂t Dt

Ng

∑ jk • gk

(12)

k=1

Eqn. 9 and Eqn. 10 are used for q and φ , respectively. The terms in Eqn. 8 due to the volumetric production of heat caused by diffusion and chemical reaction do not appear in Eqn. 12. Therefore, Eqn. 11 is not used to simplify Eqn. 12. The flux of enthalpy due to diffusion in Eqn. 8 must be explicitly evaluated and added to the heat flux caused by conduction in order to determine the total diffusional heat flux. The mixture specific enthalpy can be related to the partial specific enthalpies by Eqn. 13. For ideal gases, the partial specific enthalpy is equal to the pure component enthalpies, which are not functions of the total pressure. Ng

h(T , P) =

∑ hˆ k(T )Y k

(13)

k=1

The dependent variable most easily used with Eqn. 12 is the temperature. If the mixture enthalpy itself were used as the dependent variable, Eqn. 13 would have to be inverted to obtain the temperature. Also, the temperature appears explicitly in Eqn. 9. Because the total derivative appears on the left hand side of Eqn. 12, the enthalpy can be considered a conserved quantity. Note, this does not occur for Eqn. 8 since Cˆ p , a complicated function of the temperature and composition, appears outside of the time and convective derivatives. For discretization schemes that employ integral balances over control volumes, such as the control volume finite element methods, local as well as global conservation of ρh can be proven. For the Galerkin finite element method, conservation exists only on a global basis (see Appendix A:).

3.4

Species Mass Transport Equation The governing transport-reaction equation for each molecular species, k, is expressed by Eqn.

14. ∂ ( ρY k ) ˙k ----------------+ ∇•( ρuY k ) = – ∇• j k + W k ω ∂t

16

k = 1, ..., N g – 1

(14)

˙ k is the molar production rate of species k from gas-phase reactions, jk is the flux of species Here, ω k due to diffusion relative to the mass-averaged velocity, u . As described above, for a CHEMKIN material type, there are N g – 1 continuity equations for the molecular species; the continuity equation for the special species is replaced by Eqn. 2, the requirement that the mass fractions Yk sum to unity. Therefore, that single species in the mechanism employs a different equation to calculate its mass fraction. For the Dixon-Lewis multicomponent diffusion algorithm, this substitution does not cause any loss of accuracy. However, when the mixture-averaged diffusion coefficients are used, the effective continuity equation for the special species may have a different type and generally larger discretization error than other species in the mechanism. An “effective continuity equation” for this species, Ng, can be derived by taking the sum of all species continuity equations (Eqn. 14), k = 1, …, N g – 1 , subtracting it from the total continuity equation (Eqn. 7), and then invoking Eqn. 2. To minimize the errors in this “effective continuity equation” for the special species, the special species should be chosen to be the species with the largest mass fraction. Eqn. 2 doesn’t have to be used to ensure that the sum of the mass fraction equals one; it is implied by the continuum equations and by the property that the sum of the diffusion velocities and species mass production rates is zero. This can be seen by summing Eqn. 14 over all species and subtracting the total continuity equation, Eqn. 7. The resulting equation is Eqn. 15. ∂ ρ ∂t

Ng

Ng

∑ Y k + ρu • ∇ ∑ Y k k=1

= 0

(15)

k=1

If Eqn. 2 holds rigorously as an initial condition, Eqn. 15 ensures that the sum remains equal to one everywhere for all time. The presence of reacting surfaces, roundoff error, discretization error, and time-step truncation error, however, changes this result in the numerical problem, necessitating the use of Eqn. 2. The mass fractions Y k are the dependent variables solved for in the species conservation equation. However, mole fractions are used for specification of boundary conditions and source terms, ˙ k . The conversions between mass and mole fractions are shown in Eqn. 16. ω Wk Y k = ------- X k W

Ng

where

W =

∑ XkW k k=1

1 = --------------------------N

(16)

g

∑ Yk ⁄ Wk k=1

Other material types also use Eqn. 14 for the mass transport-reaction equation. However, they default to a different formula for the conversion of mass fraction to mole fraction. For the NEWTONIAN, BOUSSINESQ, and SOLID material types, W is assumed to be constant. Then, W k becomes a constant multiplicative factor in Eqn. 14, which can be factored out after some substitutions of definitions. The assumption of constant W is appropriate for dilute advection-diffusion of trace species in liquids and solids. When the values of W and W k are defined to be unity, the mole and mass fractions of a species become identical, and the dependent variable in Eqn. 14 can be considered to be the mole fraction.

17

3.4.1 Diffusion Velocities In Eqn. 14, jk can be written in terms of the diffusion velocity for species k, Vk. j k = ρY k V k

(17)

Several different approximations for Vk are used within MPSalsa depending upon the material type. For the NEWTONIAN, BOUSSINESQ, and SOLID material types, jk is expressed by Eqn. 18. Dk j k = – ρD k ∇Y k orV k = – ------ ∇Y k Yk

(18)

The default for these material types is to assume that D k is constant but the user can override the default and make it a user-specified function of the solution. For the CHEMKIN material type, two different approximations for the diffusion velocity are used in the code: the mixture-averaged diffusion approximation and the Dixon-Lewis formulation [20]. In the full Dixon-Lewis formulation, Vk is expressed in terms of the ordinary multicomponent T diffusion coefficients, Dkj, and the thermal diffusion coefficient, D k . 1 V k = -----------XkW

Ng

T

D k ∇T -------- -------∑ W j Dkj d j – ρY k T

(19)

j≠k

In this equation, Xk is the mole fraction for the kth species, and dj is the diffusional driving force for the jth species given by Eqn. 20 [25, 25]. Note that dj is expressed in terms of the gradient of the mole fractions instead of the mass fractions. Ng

∇P ρ d j = ∇X j + ( X j – Y j ) ------- + --- ∑ Y j Y i ( g i – g j ) P P

(20)

i=1

The second term in Eqn. 20 is the pressure diffusion term. Pressure gradients can create driving forces for separation of species with different molecular weights. However, except for applications designed to specifically use this driving force to effect a separation of isotopes, this term is usually negligible compared to other terms. The last term in Eqn. 20 is the driving force for diffusion due to differences in the body forces between species. For neutral gas transport where the only body force is gravity, this term is identically zero. In the initial implementation of the code, only the first term in Eqn. 20 is included. Other terms will be added when warranted by an application. Values for the ordinary multicomponent diffusion coefficients Dkj and the multicomponent T thermal diffusivities D k are obtained from library calls to the CHEMKIN transport parameters package [24]. Details concerning their formulation may be obtained from [24]. However, it should T be noted here that Dkj and D k have the property that the sum of the diffusive fluxes is zero. The full, multicomponent diffusion formulation is extremely expensive and possibly too expensive to be carried out in the two- and three-dimensional applications for which this code is designed. Therefore, the solution strategy concentrates on implementing approximations to the rigorous mul-

18

ticomponent diffusion formulation. A user flag is set to indicate the level of approximation to be used. The full formulation is available, however, to check the accuracy of other approximations with respect to the full multicomponent formulation. The mixture-averaged diffusion velocity formula, Eqn. 21, does not have the property that the sum of the diffusive fluxes is zero. For two- and three-dimensional applications it is, however, much less expensive. Additionally, it reduces the coupling between species equations, leading to a more efficient iterative solution of the global linear equations. T

D k ∇T D km V k = – --------- d k – --------- -------ρY k T Xk

3.5

(21)

Calculation of Diffusion Velocities

As mentioned, the cost of undertaking a full multicomponent diffusion formulation is prohibitive for two- and three-dimensional reacting flow problems. Therefore, several levels of approximation are used by the code which are similar to those used in the 1-D code, SPIN. Each of these approximations calculates the diffusion velocities, Vk, in a different manner by expressing the conservation of species mass density equation for species k in terms of a pseudo-Fickian diffusion coˆ k , and the thermal diffusion coefficient, D T , as shown in Eqn. 22 and 23. efficient, D k ∂ ( ρY k ) ˙k ----------------- + ∇•( ρuY k ) = – ∇•( j k ) + W k ω ∂t

(22)

where T ˆ k ∇Y – D T ∇ j k = ρY k V k = – ρD k k ------T

(23)

Eqn.’s 22 and 23 assume that the pressure and body-force diffusion terms are negligible. In the limˆ k is equal to the binary diffusion coefficient. The comit of a binary mixture or a dilute mixture, D bination of Eqn. 22 and Eqn. 23 has great utility as an approximate form for the Jacobian because it does not require the expensive calculation of all the cross-coupling terms. The Jacobian entries for the row corresponding to an unknown for the mass fraction of species k will be non-zero only ˆ k is treated as for unknowns corresponding to the mass fraction of species k. This assumes that D a constant in the calculation of the Jacobian, and that the dependence of ρ on Y k is also not taken ˙ k , will tend to fill into account. Of course, Jacobian entries corresponding to the reaction term, ω in those same entries.Various approximations to the multicomponent diffusion formulation, as well as the rigorous multicomponent diffusion formulation, can now be put in the pseudo-Fickian ˆ k can be obtained from the full multicomponent diffusion form. For example, an expression for D diffusion form, Eqn. 19, when forced diffusion and body-force diffusion are negligible. N

Wk g –1 ˆ D k = -------------------------- -------2 ∑ W j D kj ∇X j • ∇Y k ∇Y k • ∇Y k W j≠k

(24)

19

When the full multicomponent diffusion formulation is used, it is expected that Eqn. 19 will be used to calculate the diffusional velocities in the residuals. However, since the multicomponent difˆ k can be efficiently calculated fusion velocity has been calculated for evaluation of the residual, D ˜ , for use in the Jacobian as follows. If the multicomponent diffusion velocity is represented as V k Eqn. 25 defines the pseudo-Fickian diffusion coefficient. T

Dk  ˜  –1 ˆ k = ------------------------- Y k V D k + ------- ∇T  • ∇Y k ρT ∇Y k • ∇Y k   ˆ1 = D ˆ2 = D = In the binary limit, it can be shown from Eqn. 24 that D 12 nent diffusion coefficient reduces to the binary diffusion coefficient.

(25)

D 12; the multicompo-

Two simplified approximations to the full multicomponent diffusion formulation that are more computationally economical will now be described. The first approximation is the mixtureaveraged diffusion approximation [26, 27]. The second approximation is a more computationally intensive approximate solution of the Stefan-Maxwell equations introduced by Oran and Boris [28]. For steady-state problems, it is expected that the user first obtain a solution to the equations employing the mixture-averaged diffusion approximation. Then, if more accuracy is desired, the full multicomponent diffusion equations may be used. The Stefan-Maxwell equations have not been implemented in the code. It is expected that the mixture-averaged diffusion coefficient formulation will get the most use in the code. In the mixture-averaged diffusion formulation, Eqn. 19 for Vk is replaced by Eqn. 26. T

D km D k ∇T V k = – ---------- d k – --------- -------Xk ρY k T

(26)

The mixture-averaged diffusion coefficient, D km , can be obtained directly from a call to the CHEMKIN transport library. It is a simple function of the composition and the binary diffusion coefficients, Eqn. 27. 1 – Yk D km = ----------------------------X ⁄ D ∑ j jk

(27)

j≠k

In the above equation, Djk is the binary diffusion coefficient between species j and k. Dkm can be ˆ k by equating expressions for Vk. Assuming that forced diffusion and bodyformally related to D force diffusion are negligible, Eqn. 28 results. W ∇ X k • ∇Y k ˆ k = D ------k- -------------------------D km W ∇Y k • ∇Y k

(28)

ˆ k in formulating the Jacobian needed to relax the residuals when the mixEqn. 28 is used for D ˆ 1 is not equal to ture-averaged diffusion coefficient is used in the residuals. In the binary limit, D ˆ D1m ( D 1 = ( W D 1m ) ⁄ W 2 ), because Dkm is not equal to the binary diffusion coefficient.

20

The mixture-averaged diffusion coefficient, Dkm, has the unfortunate property that it doesn’t ensure that the diffusion fluxes sum to zero. Thus, a correction velocity is needed to ensure that this fundamental condition holds [25]. In this approach, the diffusion velocity vector is redefined to be ˆ k + Vc . Vk = V

(29)

ˆ k is the ordinary diffusion velocity computed by the various methods given above, and Vc is a V constant correction factor (independent of molecular species), defined by Eqn. 30. Ng

c

V = –

∑ Y k Vˆ k

(30)

k=1

The addition of the correction velocity to the diffusive flux expressions either requires additional terms in the Jacobian or the calculation of the entire diffusion term in the Jacobian by numerical differentiation. The current implementation of the code chooses the latter.

3.6

Implementation of Gas Phase Reactions

The gas-phase reaction mechanism enters into MPSalsa through the volumetric production ˙ k , in the species conservation rate for the kth species due to homogeneous chemical reaction, ω ˙k equations and in the temperature representation of the internal energy conservation equation. ω is calculated using the CHEMKIN package [3]. This modular approach to programming complex chemical mechanisms has found a great deal of use in the combustion and CVD community [2, 13, 29, 30] because it allows separation of the specification of a complex reaction mechanism from the programming of the numerical representation of the continuity equations. Additionally, different types of reactions (e.g., reversible and irreversible reactions, unimolecular reactions whose rate constant is parameterized by a Troe form, bimolecular reactions, third body reactions with enhanced third body collision efficiencies, and/or lumped kinetics expressions appropriate for the description of overall combustion processes) may be integrated into the numerical code without having to include complex reaction mechanisms. Moreover, changes to the mechanism do not induce changes in the numerical code, and correspondingly, mechanisms developed for one numerical code may be applied in any other numerical code conforming to the CHEMKIN interface. ˙ k used by CHEMKIN [3]. Not all of The following is a brief review of the formulation of ω the complexity possible with CHEMKIN will be discussed here. Consider N R elementary reversible or irreversible reactions involving Ng chemical species that can be represented by Eqn. 31. Ng

Ng

∑ ν'ki χk ↔ ∑ ν''ki χk k=1

i = 1, …, N R

(31)

k=1

ν' ki is the stoichiometric coefficient of the kth species for the forward direction of the ith gas-phase reaction; it is defined as a non-positive number. ν'' ki is the stoichiometric coefficient of the kth species for the reverse direction of the ith gas-phase reaction; it is defined as a non-negative number.

21

The possibility of non-integer stoichiometric coefficients is allowed as long as the reaction satisfies charge and elemental balances. The χ k represents the chemical symbol for the kth species. The pro˙ k for the kth species can be written as a summation of the rate-of-progress variables, duction rate ω q i , for all reactions involving the kth species, Eqn. 32. NR

˙k = ω

∑ νki qi

v ki = ν'' ki + ν' ki

where

(32)

i=1

The default in CHEMKIN is to assume mass action kinetic rate constants. For this case, the rate of progress variable, q i , for the ith reaction is given by the difference of the forward rates and the reverse rates, expressed by Eqn. 33 where q i has units of mol cm -3 s-1. Ng

qi =

f ki

∏ [Xk]

ν' ki

Ng



k=1

r ki

∏ [Xk]

ν'' ki

(33)

k=1 f

r

Here, [ X k ] is the molar concentration of the kth species, and k i and k i are the forward and reverse rate constants for the ith reaction, respectively. The forward rate constants for the N R reactions default to having the following extended Arrhenius temperature dependence: –E i β f k i = A i T i exp  --------  .  RT 

(34)

Other expressions for the reaction rate constants, Eqn. 34, are also allowed, such as fall-off behavior parameterized by a Troe form, Landau-Teller reaction rate forms, and third body reactions. The r reverse rate constants k i are generally (but not necessarily) related to the forward rate constants c through the concentration-based equilibrium constant for the ith reaction, K i , according to Eqn. 35. f

r ki

k = -----icKi

(35)

c

K i is in turn related to the temperature, the net molar production rate of gas production during the reaction, and the Gibbs free energy of reaction. Thermodynamic information for the equilibrium constant is calculated from CHEMKIN’s species thermodynamic information. Thermodynamic information is in a format [31] similar to that used by Gordon and McBride [32] for the thermodynamic database used in the NASA chemical equilibrium program. Chemical reaction mechanisms usually consist of stiff modes, i.e., reactions which are fast compared to other time scales in the problem. Therefore, it is imperative that the Jacobian terms ˙ k be available. The current procedure is to calculate the ω ˙ k source term only at nodes, in order for ω ˙ k is interpolated throughout the element using elemental to reduce the expense of this step. Then, ω basis functions. The Jacobian contributions for the source terms due to reaction are currently calculated at the nodes in the element via numerical differencing. Details of their implementation are discussed in Appendix A.

22

3.7

Implementation of Surface Phase Reactions

Surfaces where reactions take place create additional source and sink terms for gas-phase species. The boundary conditions for the gas-phase continuity equations for species must specify the total flux of the species at the domain interface. For the case where the interface is stationary and the growth or etching due to surface reactions can be considered not to move the interface, this boundary condition for species k can be expressed by Eqn. 36. n • ( ρuY k + j k ) = – s˙k W k

(36)

The left side of Eqn. 36 represents the total flux of species k, both convective and diffusive. The first term on the left-hand side is the Stefan flux term where n is the outward facing normal to the domain and j k represents the net diffusive flux of k from all diffusive processes, including thermal diffusion. The right-hand side represents the net destruction of gas-phase species k due to chemical reaction. Therefore, s˙k represents the net molar production rate of gas-phase species k due to chemical reaction. Integration by parts, carried out in the Galerkin formulation (discussed in Appendix A), leads naturally to surface integrals of the normal component of j k multiplied by the nodal basis functions. Thus, in applying boundary conditions to the kth gas species continuity equation, the normal component of the diffusive flux for species k is replaced by the right hand side of Eqn. 37. n • j k = – s˙k W k – n • ρuY k

(37)

Eqn. 37 can be further simplified by summing Eqn. 36 over all gas-phase species and using the property that diffusive fluxes must sum to zero to yield an expression for the Stefan flow, Eqn. 38. Eqn. 38 can then be used in Eqn. 37 to yield Eqn. 39. Ng

n • ρu = –

∑ s˙k W k

(38)

k=1 Ng

n • j k = – s˙k W k +

∑ s˙k W k

Yk

(39)

k=1

Eqn. 39 is used within MPSalsa for specification of boundary conditions for gas-phase species equations for the case of a reacting surface. Additionally, Eqn. 38 is used for specification of the normal boundary condition for the momentum equation. The tangential boundary condition for the momentum equation for reacting surfaces is set to the no-slip condition. Thus, the problem is reduced to the calculation of s˙k , k = 1, 2, …, N g . For CHEMKIN material types, s˙k is supplied by the SURFACE CHEMKIN package [4]. However, they are functions of additional unknowns corresponding to surface site fractions of surface phases and bulk mole fractions of bulk phases where each surface phase represents a different type of surface site and each bulk phase represents a different type of bulk mixture. (The reader is referred to the manual for the SURFACE CHEMKIN package [4] and to the manual for the Surface PSR program [33] for a more complete description.) Thus, the calculation of s˙k demands the solution of a subproblem at each node on the reacting surface to calculate the values of the extra unknowns corresponding to the state of the surface. The resulting non-linear system of equations is solved using Newton iteration. Since the subproblem is 23

solved at each node, it is completely local to a processor and, thus, requires no additional communication when run on parallel computers. We now describe the equations that comprise this subproblem. Let Z k(n) be the surface site fraction of the kth surface species in the nth surface phase. Let Γ n be site density for the nth surface phase (e.g. mol cm-2). Let ck(n) be the concentration of the kth surface species in the nth surface phase (e.g. mol cm-2). Then, Eqn. 40 is the conservation equation expressing the continuity balance for the kth surface species in the nth surface phase. d( AW k c k(n) ) f l surf -------------------------------- = AW k s˙k , k = K s (n) ,...,K s(n) , n = 1, …, N phase dt

(40)

Here, s˙k is the production rate from surface reactions for the kth surface species, A is the surface f l area, and W k is the molecular weight of the kth surface species. K s (n) and K s(n) are the indices th for the first and last surface species in the n surface phase, respectively. Also, ck(n) can be related to Z k(n) by Eqn. 41. Γ n Z k ( n) c k(n) = -----------------σk

(41)

Here, σ k is the number of surface sites the kth species covers. Substituting Eqn. 41 into Eqn. 40 and assuming A is not a function of time yields the equation for Z k(n) as a function of time. In general, Γ n can also be a function of time and this must also be taken into account. d Z k ( n) dΓ n f l surf Γ n --------------- = σ k s˙k – σ k --------- , k = K s (n) ,...,K s(n) , n = 1, …, N phase dt dt

(42)

For any valid surface mechanism, the following equation also holds true for each surface phase n, regardless of the Z k(n) used. l

K s(n)



k=

f K s ( n)

dΓ σ k s˙k = --------ndt

(43)

Eqn. 43 is called the surface site conservation equation. For most reaction mechanisms, the righthand side of Eqn. 43 is identically zero. If this is not the case, MPSalsa expands the solution vector at each surface to include Γ n and uses Eqn. 43 to solve for the concentration of surface sites for phase n as a function of time. On each surface, the sum of the surface site fractions must equal one. l

K s(n)



surf

Z k(n) = 1 , n = 1, ..., N phase

(44)

f

k = K s (n)

This implies that the use of Eqn. 40 leads to a singular Jacobian for the steady state case, if used for all surface species site fractions in a surface phase. Thus, one of the surface species balance 24

equations, Eqn. 40, is replaced with Eqn. 44 for each surface phase. This has the disadvantage that all the numerical round-off error is assigned to that one equation. Therefore, the equation corresponding to the species with the largest site fraction in the surface phase is replaced by Eqn. 44. The amount of material in bulk phases within the domain may not be in steady state; i.e., the bulk phases may be growing or etching (although their growth/etch rate is not assumed to affect either the volume or surface area within the domain). MPSalsa treats the mole fractions of bulkphase species as well as their growth/etch rates as unknowns to be solved for. The format of these equations depend on whether the bulk phase is growing or being etched. The following equations apply to a growing phase. In this case, the growth rate of the nth bulk phase, G(n), can be expressed by the following equation: d[ AL n C b(n) ] -------------------------------- = A dt

l

K b(n)

∑ G k(n)

k=

where

G k(n)

= AG (n)

(45)

f K b ( n)

= Max(s˙k, 0)

In this equation, Ln is the film thickness for the nth bulk phase. C b(n) is the average molar concentration of the nth bulk phase; it has units of mol cm-3, G k ( n ) is the growth rate of the kth species in the nth bulk phase, and s˙k is the production rate of the kth species returned from SURFACE CHEMKIN. It is a function of the gas phase concentrations, pressure, temperature, surface site concentrations, and the bulk phase activities. Having s˙k less than zero for some species, while it is greater than zero for other species is not appropriate for a growing bulk phase. One positive value of s˙k for a bulk phase signals that particular phase is growing. b

For a growing phase, X k (n) , the instantaneous mole fraction of the kth bulk-phase species in the nth bulk phase, is determined from the relative growth rates of all species in that phase, Eqn. 46. 0 =

G k ( n ) – X bk(n) G ( n )

(46)

The condition Max ( s˙k, 0 ) may violate the overall elemental balance condition. However, in b practice, this does not occur because X k (n) for such a species is set to zero by Eqn. 46. Then, only nonphysical mechanisms involving zeroth-order destruction of a bulk species could possibly create b the situation where s˙k < 0 and X k (n) = 0 . If all s˙k for a particular bulk phase are less than zero, that bulk phase is undergoing etching. The user can specify whether a particular phase is expected to be etched and MPSalsa solves a different set of equations for the bulk-phase components for that bulk phase. In this case, the user must also supply the initial composition of the bulk phase to be etched. The time-dependent equations used for the bulk-phase mole fractions and etch rates in the nth bulk phase undergoing etching are then given by Eqns. 47 and 48. b

b

0 = X k (n) INITIAL – X k (n)

G k(n)

(47)

= s˙k

(48)

25

b

Here, X k (n) INITIAL is the user-supplied initial estimate for the mole fraction of species k in bulk phase n, assumed to be normalized so that the sum over all bulk-phase species is one. The idea is that the initial phase is being etched away congruently. Incongruent etching, within the context of a single phase, is not allowed, at least at the level where it affects the concentrations of bulk species. b

In order to specify the thermodynamic information needed for bulk phases, the activities a k of the bulk-phase components must be determined. These are the quantities in SURFACE CHEMKIN that appear in the rate expressions for surface reactions. This is done within the code by calling a subroutine that users can modify to specify their own relationships between the bulk activities and the bulk mole fractions, temperature, and pressure. The default subroutine assumes a perfect solution relationship for all bulk phases, Eqn. 49, that almost never occurs in practice. b

b

b

a k (T , P, X k (n)) = X k (n) b

(49)

In summary, the extra unknowns, Z k(n) , Γ n(n) , X k (n) , and G k ( n ) are not included in the formal solution vector. Instead, a separate subproblem is solved for these unknowns as part of the calculation of the residual and Jacobian entries for the gas-phase problem. The two problems are coupled at the gas-species flux level, Eqn. 39. The surface subproblem depends on the gas-phase concentrations at the surface, while the main gas-phase species problem depends on the fluxes calculated from the surface subproblem. An advantage of this approach is that the surface subproblem calculation can be protected from nonphysical occurrences, such as negative gas-phase mole fractions, and made more robust than it would be if lumped in with the main problem. Also, advanced surface profile simulators may be incorporated into MPSalsa at a later date. These simulators model behavior at the micron feature size, and couple into “reactor simulators” such as MPSalsa, which model behavior at the centimeter or meter feature size, through the gas-phase flux b boundary condition described above. Solving a separate subproblem for Z k(n) and X k (n) , however, can create some concerns. For time-dependent reacting flow problems, difficulties typically associated with operator splitting techniques arise if the surface unknowns are allowed to have true time dependence (i.e., if they are not assumed to have a faster transient than the bulk and, thus, are assumed to be in pseudo-steady state at each time step of the gas-phase problem).

26

3.8

Summary of Transport Equations Implemented

Mixture Momentum: ∂ ( ρu ) --------------- + ∇•( ρuu ) – ∇•T – ρg = 0 ∂t

(50)

2 T T = – PI – --- µ ( ∇•u )I + µ [ ∇u + ∇u ] 3

(51)

∂ρ ------ + ∇•( ρu ) = 0 ∂t

(52)

Mixture Continuity:

Thermal Energy: Ng

( ρT ) Cˆ p ∂-------------- + ∇•( ρuT ) = – ∇•q c + φ + Q˙ – ∑ j k • Cˆ p, k ∇T ∂t k=1

(53)

Ng



∑ hk W k ω˙ k – ∇•qr k=1

2 2 1 T φ = – --- µ ( ∇•u ) + --- µ ∇u + ∇u 3 2

2

(54)

q c = – λ ∇T

(55)

Species Continuity: ∂ ( ρY k ) ˙k ----------------- + ∇•( ρuY k ) = – ∇• j k + W k ω ∂t

(56)

j k = ρY k V k

(57)

T ˆk D k ∇T D V k = – ------ ∇Y k – --------- -------ρY k T Yk

(58)

27

4. Boundary Conditions

4.1

Specification of Boundary Conditions

In general, the second-order transport-reaction equations in MPSalsa need either their dependent variables or their normal derivatives specified at all domain boundaries in order to define a well-posed problem. EXODUS II defines the concept of node sets and side sets on which these boundary conditions are applied. A node set is an arbitrary group of nodes in the domain. A side set is an arbitrary group of element sides in the domain. Only side sets establish the concept of a surface. Dirichlet boundary conditions specify the value of dependent variables. The usual conservation equation for the dependent variable identified with an element node, where a Dirichlet boundary condition is specified, is discarded and replaced with another equation for that variable. The new equation may be a function of the other independent or dependent variables in the problem. Dirichlet conditions that don’t need the concept of a surface may be applied on node sets as well as side sets. MPSalsa also allows for Dirichlet conditions to be applied as surface integrals of functions weighted by the elemental basis functions, i.e. Galerkin’s method. These surface integral Dirichlet conditions may be applied only on side sets. For example, the concept of a surface is needed to define normal and tangential vectors for normal and tangential velocity boundary conditions. Neumann and Robin (or mixed) boundary conditions impose conditions on the normal derivative of the dependent variable. This term is specified by replacing the normal derivative in the surface integral that arises from the integration by parts during the Galerkin finite element formulation (see, e.g., Eqn. C.25) with the boundary condition. Surface integral conditions may be applied only on side sets and are generally defined as being satisfied in a “weak sense”. In other words, they are satisfied only in the limit of no discretization error. The following is a discussion of the types of boundary conditions permissible in MPSalsa for each of the conservation equations.

4.2

Momentum Equations

For the fluid dynamical part of the problem, either the velocity components or the normal component of the total stress tensor must be specified on the boundary of the domain for each component of the vector momentum equation. On both side and node sets, Dirichlet boundary conditions of the form u m = f um(x, u, P, T , Y, t) , m = 1, 2, 3

(59)

may be applied to the velocity in the x -, y - and z -directions, respectively. In Eqn. 59, f um is a user-specified function of the dependent and independent variables. For these boundary conditions, the corresponding momentum equations are replaced by Eqn. 59 at all nodes of the designated node or side sets.

28

Surface integrals involving the components of the surface traction vector, τ = n • T = Tn , on a surface with normal, n , arise naturally in the Galerkin form of the momentum equations (see Eqn. C.14) and are added to the volumetric contributions of the Jacobian and residuals of all nodes on the surface. The components of the normal stress may be replaced in the surface integrals by user-specified functions f τ of the dependent and independent variables, as shown in Eqn. 60, where Φ i is the elemental shape function for node i on the surface.

∫ τm Φi dΓ

∫ f τ, m Φi dΓ

=

Γ

, m = 1, 2, 3

(60)

Γ

Boundary conditions may also be applied to the normal and tangential components of the velocity and normal stress. Each region of the boundary is associated with a unit normal to the boundary, n, and two orthonormal tangential components to the boundary, t1 and t2. Specification of the boundary condition for the momentum equations then involves specification of the velocity component or normal tensor component in each of the directions n, t1, and t2; that is, the user must τ•n τ • t1 specify either n • u or , and either t 1 • u and t 2 • u , or and τ • t2 . Normal and tangential Dirichlet boundary conditions on velocity are enforced using surface integrals along sides of elements. The surface integral form of a Dirichlet boundary condition on the normal velocity is given by Eqn. 61. For each elemental node on a surface, i , the boundary condition is multiplied by the elemental shape function Φ i . The integral over the surface of the resulting expression is the residual contribution for the corresponding component of the momentum equation for node i . Similar expressions enforce tangential velocity boundary conditions.

∫Γ ( n • u – f u ( x, u, P, T , Y, t ) )Φi dΓ n

= 0

(61)

Conditions on the normal stress in the normal and tangential directions are enforced by replacτ•n , τ • t 1 , and τ • t2 ing in the surface integrals with user-supplied funcτ•i , τ • j , and τ•k , tions, which are then rotated to derive expressions for which are needed in the surface integral terms in the x, y, and z momentum equations, respectively. For example, Eqn. 62 specifies traction boundary conditions in a 2-D geometry with n = i and t1 = j. In this examples, f τ is the user-supplied function specifying the traction vector.. 2 ∂u ∂v ∂u τ • n = τ • i = – P – --- µ  ----- + -----  + 2µ ----- = n • f τ(x, u, P, T , Y, t) 3  ∂x ∂y  ∂x

(62)

∂u ∂v τ • t 1 = τ • j = µ  ----- + -----  = t 1 • f τ(x, u, P, T , Y, t)  ∂y ∂x 

In Eqn. 62, f τ is shown as a function of all of the independent and dependent variables in the problem ( x, u, P, T , Y, t ) . A common outflow boundary condition is setting the normal stress, τ • n , to zero. This is the so called natural B.C. on the momentum equation. For the particular case of a reacting, impermeable wall, the Dirichlet boundary conditions in Eqn. 63 are applied using Eqn. 61.

29

Ng

n • ρu = –

∑ s˙k W k

(63)

k=1

t1 • u = t2 • u = 0

4.3

Total Continuity Equation

The incompressible Navier-Stokes equations are unchanged when the hydrodynamic pressure is changed by a constant. They are affected only by gradients in the hydrodynamic pressure and MPSalsa’s discrete equation set shares this property. Therefore, the pressure scale must be set eiτ•n ther implicitly or explicitly somewhere in the domain. This is achieved by specifying somewhere on the boundary since P appears in this expression (see Eqn. 62), or by setting a Dirichlet condition for P on one node in the domain.

4.4

Internal Energy Continuity Equation

Either the temperature or the normal heat flux must be specified on all boundaries of the domain. That is, either Dirichlet boundary conditions in the form of a user-supplied function or value must be specified for the temperature, or surface integral boundary conditions involving the heat conduction must be used. The expression in the surface integral resulting from the Galerkin integration by parts is the normal component of the heat flux vector, n • q c , where q c = – λ∇T (see Eqn. C.6). The user supplies a function that is substituted for n • q c in the surface integral. Inflow boundary conditions for the energy equation are usually specified by a Dirichlet condition on the temperature. For cases where the energy balance at a surface must be calculated, Eqn. 64 is a useful starting point in the derivation of energy boundary conditions based on heat balances. –

FLUX + PRODUCTION Γ = FLUX

+

(64) –

The heat flux to the boundary from within the solution domain is defined as FLUX . This, plus the energy stored at the interface, PRODUCTION Γ , should be equated to the heat flux exiting + the domain, FLUX . For the convection of enthalpy inlet boundary condition, PRODUCTION Γ is zero but the flux terms are defined by Eqn. 65. An extra convective heat transfer term, h c ( T – T o ) , is added to the inflow heat flux, on the outer side of the domain. In MPSalsa, the user supplies a function returning the value of n • q c as determined by Eqn. 64 and 65.

30

N

N

g  g  ˆ FLUX = – n •  ∑ ρY k h k u + λ ∇T + ∑ h k j k + q r k = 1  k=1





(65)

Ng

  + FLUX =  ∑ ρ o Y k, o hˆk(T o)u o + q r k = 1 

+ hc ( T – T o ) +

For boundary conditions corresponding to outflow areas, where neither the energy flux nor the temperature is known before hand, the specification of a zero normal temperature derivative is used (natural b.c.). n • ∇T = 0

(66)

For boundary conditions corresponding to solid walls where reactions may be occurring, Eqn. – 64 may be used to obtain a heat balance. FLUX is given by Eqn. 65 and PRODUCTION Γ is nonzero due to the growing or etching film at the interface. l

l

Kb

Ks



PRODUCTION Γ =

k=

s˙k W k h k +

f Ks

∑ k=

s˙k W k h k + Q˙ Γ

(67)

f Kb

PRODUCTION Γ includes terms due to the storage of energy due to surface and bulk-phase species and Q˙ Γ is the heat input to the boundary from external sources (e.g., resistive heating). + Typically, FLUX is specified by a heat transfer coefficient combined with radiative heat input from a black body at a known temperature. However, its exact specification is left undefined at this – point. The enthalpy terms in FLUX and PRODUCTION Γ may be combined with reacting wall boundary conditions on the species conservation equations (Eqn. 37) to yield Eqn. 68. l

Kb

– n • ( λ ∇T + q r ) – –

∑ s˙k W k hk

+ = Q˙ Γ + FLUX

(68)

k=1

The sum in Eqn. 68 is over all species defined in the problem: gas, surface, and bulk. For phase change-type reactions, the second term in Eqn. 68 can be identified with the latent heat of the phase change. Radiation contributions, q r , appear naturally in surface integral expressions for the heat flux. Currently, an MP gray body radiation treatment is under development and will be presented at a later time.

4.5

Gas-Phase Species Continuity Equations

Several types of boundary conditions may be specified on Y k , k = 1, …, N g . Theoretically, either the value of Yk or its normal derivative must be specified on a boundary. However, for low pressure systems where diffusive transport dominates, Dirichlet conditions on the species equations are discouraged as a means of specifying the flow rate of species k into the system. The actual

31

flux of species k into the domain, which consists of both convective and diffusive contributions, will be quite different than the intended flux into the domain. Therefore, flux-based conditions should be used on all boundaries of the domain for these systems. For boundary conditions corresponding to inflow areas where the flow rates of the gas-phase species are known, the flux of species k is specified by what is known as Danckwerts’ boundary condition: n • ( ρY k u + j k ) = ρ o u o Y k, o

( k = 1, ..., N g ) ,

(69)

where ρ o , u o and Y k, o are user-specified values. For boundary conditions corresponding to solid walls where reactions may be occurring, the flux of species k to the wall should be equated with the negative of the net production rate of species k at the wall. n • ( ρY k u + j k ) = – s˙k W k

( k = 1, ..., N g )

(70)

For boundary conditions corresponding to solid walls, where no reactions are occurring, the net flux of species k should be set to zero. n • ( ρY k u + j k ) = 0

( k = 1, ..., N g )

(71)

For boundary conditions corresponding to outflow areas, where neither the flux nor the concentration of species k is known, the specification of a zero normal diffusion velocity may be employed. n • jk = 0

( k = 1, ..., N g )

(72)

The boundary conditions in Eqn. 70-72 are incorporated into the finite element equations representing the continuity equation for species k via the boundary integral involving ( n • j k ) that appears from the integration by parts of the diffusive flux term (see Eqn. C.7). Specifically, ( n • j k ) is replaced with the appropriate terms from Eqn. 70-72 expressed via a user-supplied Y function f k as in Eqn. Y

n • jk = f k

( k = 1, ..., N g )

(73)

As with any Neumann or Robin boundary conditions in the finite element method, these boundary conditions are satisfied only in the limit that the discretization error goes to zero. Also, if a determination of the flux of species k is required at a reacting solid wall where Eqn. 70 is used, the flux should be evaluated using the right hand side of Eqn. 70 instead of the left hand side. The accuracy in Yk is one order of the mesh discretization size greater than the accuracy in the derivatives of Yk. Y

For non-CHEMKIN material types, Dirichlet boundary conditions of the form Y k = f k , k = 1, …, N g , and flux boundary conditions of the form in Eqn. are supported.

32

5. Finite Element Approximation of the Transport Equations

The governing transport Eqn.’s 50-58 are approximated by a Petrov-Galerkin finite element method (PGFEM). The summary presented here is intended to provide a sufficiently detailed discussion of the FE development and a practical formulation background to discuss the numerical algorithms that are used to solve the resulting linear algebra problems. The finite element procedure begins by dividing the physical domain of interest, Ω , into N e simply shaped regions Ω e called finite elements. Within each of these elements, the dependent variables ( u 1, u 2, u 3, P, T , Y k ) , k = 1, …, N g , are interpolated by continuous functions of compatible order, in terms of values to be determined at a set of global node points. To develop the FE equations for these nodal unknowns, we present the finite element expansion in terms of global interpolation functions. This development differs from an elemental basis approach only in the interpretation of the summation scope and the resulting domain of integration of the inner product. Using this approach simplifies the resulting discussion of the node-based matrix-fill algorithms in the parallel implementation of the code.

5.1

Interpolation Functions and Quadrature Rules

Within each element the mixture velocity, temperature, species mass fractions, and hydrodynamic pressure are approximated by the expansions in Eqn. 74. N

u l ( x, t ) =

∑ ( uI )J ( t )ΦJ ( x )

l = 1, 2, 3

(74)

J=1 N

P ( x, t ) =

∑ PJ ( t )ΦJ ( x ) J=1 N

T ( x, t ) =

∑ T J ( t )ΦJ ( x ) J=1 N

Y k ( x, t ) =

∑ ( Y k )J ( t )ΦJ ( x )

k = 1, …, N g

J=1

Here, Φ J ( x ) is the standard polynomial finite element basis function associated with the Jth global node, N is the total number of global nodes in the domain, and N g is the number of gas-phase species. The u 1 , u 2 , and u 3 components of velocity correspond to velocity in the x -, y -, and z directions, respectively. Equal order interpolation of all variables is used. In these and the following expansions, global interpolation functions are denoted with uppercase indices as in the expansions above. The only exception to this convention is the use of a lower case index k to denote the species number.

33

Thermodynamic and transport properties, as well as volumetric source terms, are interpolated from their nodal values using the finite element shape functions. For example, Eqn. 75 represents the computation of density at a point x . The density is not evaluated from the equation of state with values of the dependent variables at x but instead is computed at global nodes J = 1, …, N with values of the dependent variables at the global nodes, and the elemental shape functions are used to interpolate the density at x . N

ρ ( x, t ) =

∑ ρJ ( t )ΦJ ( x )

(75)

J=1

Evaluation of volumetric integrals is performed by standard Gaussian quadrature. For quadrilateral and hexahedral elements, two-point quadrature (in each dimension) is used with linear basis functions, while three-point quadrature is used for quadratic interpolated elements. For example, for tri-linear hexahedral elements, eight Gaussian quadrature points within an element are used to evaluate its volumetric integrals.

5.2

Evaluation of Surface Integrals

Evaluation of surface integrals is performed by standard Gaussian quadrature on the side of the element. As with the volumetric integrals, two-point quadrature (in each direction) is used with linear shape functions, while three-point quadrature is used with quadratic shape functions. For example, for a three-dimensional problem with linear shape functions, four Gaussian quadrature points located on the side of an element are used to evaluate its surface integrals.

5.3

Matrix Equations

Substitution of the dependent variable expansions, Eqn. 74, into the governing transport equations (Eqn. 50-58) yields a set of residual equations, Eqn. 76. Momentum: (u)

Rl

with

= f ul ( Φ J , u J , P J , T J , Y J ) , l = 1, 2, 3 , J = 1, …, N R

(u)

(u)

(u)

(u)

≡ R1 i + R2 j + R3 k

Mixture Continuity: (P)

= f P ( ΦJ, uJ, PJ, T J, YJ )

(T )

= f T ( ΦJ, uJ, PJ, T J, YJ )

R

, J = 1, …, N

Thermal Energy: R

34

, J = 1, …, N

(76)

Species Continuity: (Y )

Rk

= f Y k ( ΦJ, uJ, PJ, T J, YJ )

, k = 1, …, N g , J = 1, …, N

where the R’s denote the resulting errors or residuals from using the finite element approximation in the continuous governing equations. The Galerkin form of the method of weighted residuals [34] reduces these errors in a weighted sense, by making the residuals orthogonal to the interpolation functions. The Petrov-Galerkin pressure stabilization formulation used in MPSalsa slightly modifies the Galerkin residual by employing an additional term in the momentum equation weighting function. This weighting vector includes the standard Galerkin weighting function along with a term that is proportional to the gradient of the basis functions. This projection method allows the use of equal order interpolation for velocity and pressure without producing spurious pressure modes in the solution of the incompressible flow problem. The PGFEM used here follows the work of Hughes et al. [34] and Tezduyar et al. [35]. The resulting orthogonality relations produce the PGFEM residuals at the Ith global finite element node, Eqn. Fl

(u)

F

F

( U , P, T , Y )

(P)

(T )

I

( U , P, T , Y ) ( U , P, T , Y )

∫ Rl

Ω I

I

(Y )

F k ( U , P, T , Y )

(u)

=

I

=

∫R

(P)

∫R

(T )

Φ I dΩ = 0

(Y )

Φ I dΩ = 0



=

∫ Rk



Φ I dΩ + ρτ ∫ ∇Φ I • R

(u)

dΩ = 0





=

l = 1, 2, 3

Φ I dΩ = 0

k = 1, …, N g

For Eqn. , the vectors of global unknowns are defined as U

T

= ( ( u1 )J, ( u2 )J, ( u3 )J )

T

P = ( PJ ) T

T

Y

T

J = 1, …, N ,

(77)

J = 1, …, N ,

= (T J )

J = 1, …, N ,

= ( ( Y 1 ) J , ( Y 2 ) J , ( Y 3 ) J , …, ( Y N g ) J )

J = 1, …, N ,

where N is the total number of global nodes and N g is the total number of gas-phase chemical species. The constant ρτ in Eqn. is a stability constant defined in [34, 35] and a detailed description of the choice for ρτ is given in Appendix B. The manipulations of the integral equations (Eqn. ) to produce the system of discrete matrix equations are presented in Appendix C. Below, we summarize the result of applying a Newton linearization to solve the system of nonlinear governing transport equations. We have chosen to leave the explicit representation of the time derivative terms unaltered for purposes of the following general discussion. As shown in Appendix C, this term can easily be approximated with various levels of accuracy. The results of applying a Newton

35

linearization to the governing transport equations can be represented in matrix form. In these equations the overbar variables refer to the approximate solution generated by the previous step of the inexact Newton iteration.

Momentum transport:

M (ρ)U˙ + C(U)U + D(U)U + K



(U ) GUl m U



(U ) GP m P



(Um)

N

g   U – Q P – B T T + ∑ B Y k Y k  g   k=1

(U ) GT m T

T

Ng



(Um)

∑ GY

k

Y k = –F

(Um)

(78)

( U , P, T , Y )

k=1

Mixture continuity: Q(ρ)U + RU

(79)

T (P) (P) (P) (P) + ρτ Q(ρ) U˙ + C (U)U + D (U)U + L P – B T T –

Ng

(P)

∑ BY

k

Yk

k=1

= –F

(P)

( U , P, T , Y )

Thermal energy transport: M (ρCˆ p)T˙ + C (T )

(T )

(U)T + D

(T )

(T )

(T )

(T )U + K

(T )

(T )

T+E

(T )

T

(80)

(T )

+ G U U + G P P + ( G T – ℑ T + ℜ T )T Ng

N

g  (T ) (Y )  (T ) (T ) + ∑ G Y i – ℑ Y i + ∑ ℜ Y i k hˆ k  Y i = – F ( U, P, T , Y )  i = 1 k=1

Species transport of the kth species: (Y ) (Y ) (Y ) ˙ M (ρ)Y k + C(U)Y k + D k (Y )U + K k Y k + S k (T )T (Y )

(Y )

(Y k)

+ GU k U + G P k P + ( GT

36

(Y )

– ℜ T k )T

(81)

Ng

+

(Y k)

∑ ( GY

l

(Y )

– ℜ Y l k )Y i = – F

(Y k)

( U , P, T , Y )

l=1

In Eqn. 78-81, U , P , T , and Y are given by Eqn. 77. The details of the discretization scheme can be found in Appendix C. Finally, we present a single matrix equation for the discrete transportreaction equations. This fully reduced system is useful for our description of the time-stepping algorithms that follows. The matrices in this system can be easily obtained by comparison of Eqn. with Eqn. 78-81. R ( V˙ , V ) = M u V˙ + K ( V )V + F ( V ) = 0

(82)

U V = P . T Y

(83)

where

Here, R ( V˙ , V ) is the residual for the equation system, which is a nonlinear implicit function of the dependent variables and their time derivatives. Substituting the appropriate approximation for the time derivative produces the final system of nonlinear residual equations to be solved at each time step, t n + 1 , Eqn. 84. R ( V˙ n + 1, V n + 1 ) = 0

(84)

This system is solved using the inexact Newton scheme described in Section 6 and Appendix C.

37

6. Solution Procedures

In this section, we present the general procedures used in MPSalsa for the steady state and the time dependent solution of equations that describe the discrete problem. The choice of numerical methods in MPSalsa has been made from the standpoint of robustness, efficiency of implementation on parallel architectures, and the ease of including new solution kernels. The major solution kernels used in MPSalsa are the first- and second-order implicit time integration routines, an inexact Newton procedure and the linear system solvers of the Aztec [10,11] parallel Krylov solver library, developed in conjunction with MPSalsa. Below we summarize the properties of the discrete matrix problem Eqn. 82 and consider the details of the major solution kernels in MPSalsa. First, we give a brief overview of the implementation of the unstructured finite element method on multiple processors, since this aspect underlies much of the discussion and implementation of the solution algorithms for the linear system.

6.1

Implementation on Multiple Processors

MPSalsa is designed to solve problems on massively parallel (MP) multiple instruction multiple data (MIMD) computers with distributed memory. For this reason the basic parallelization of the finite element problem is accomplished by a domain partitioning approach. The initial task on an MP computer is to partition the domain among the available processors, where each processor is assigned a subdomain of the original domain. It communicates with its neighboring processors along the boundaries of each subdomain. There are two fundamental ways to partition the FE domain among processors: either element or node assignment. Each method has its own advantages and fundamentally affects the solution strategies and interprocessor communications. Dividing the mesh according to elements quite naturally can lead to an element-by-element (EBE) solution scheme, whereas dividing the mesh according to nodes leads most naturally to a fully-summed distributed matrix solutions. In the EBE case, each element’s matrix is stored separately and is not summed with its contributions from neighboring elements. All matrix-vector operations are performed with these dense elemental block matrices and the vector result is obtained only after summing over all elements. This scheme substantially increases the matrix storage requirements and the amount of computation needed relative to fully-summed distributed matrix solution strategies. For example, for 3-D linear hexahedral elements, this method requires approximately 60% more storage and greater than three times as many floating point computations are required for the EBE approach. Although the larger block sizes associated with the EBE approach may yield an increase in the number of operations performed per second, this improved performance is unlikely to compensate for the increased operation count. Because of this, nodal decomposition was chosen in MPSalsa to allow the implementation of computationally efficient, minimum flop algorithms for the matrix-vector multiply kernel. Also, storing the fully summed equations allows the use of robust general preconditioners, such as domain decomposition incomplete factorizations and direct sparse subdomain solvers. The parallel solution of a particular FE problem proceeds as follows. At the start of the problem, each processor is “assigned” a set of finite element nodes that it “owns.” A processor is responsible for forming the residual and the corresponding row in the fully summed distributed ma-

38

trix for the unknowns at each of its assigned FE nodes. To calculate the residual for unknowns at each assigned node, the processor must perform element integrations over all elements for which it owns at least one element node. To do this the processor requires 1) the local geometry of the element and 2) the value of all unknowns at each of the FE nodes in each element for which it owns at least one node. The required elemental geometry is made available to the processor through the initial partitioning and database distribution part of the algorithm. Here, a broadcast of all information in a serial EXODUS data base to all processors is used in MPSalsa. Then, each processor extracts its geometry information form the FE database. In addition to the broadcast algorithm, MPSalsa has the capability to use a parallel FE database [9] for geometry input as well as all parallel I/O. The unstructured interprocessor communication of FE unknowns is handled by an Aztec routine that exchanges the necessary interprocessor information [11]. Figure 6-1, which depicts a partitioning scheme of an unstructured mesh, graphically represents the above concepts. An unstructured mesh is divided into four regions by assigning ownership of the nodes. Nodes in each processor are classified as “border” and “internal” nodes, at which border and internal unknowns, respectively, are defined. Border unknowns are those unknowns whose values must be communicated to neighboring processors so they may complete their element integrations; the remaining “owned” unknowns on a processor are designated as internal unknowns. Those unknowns required for a processor’s element integrations but assigned to a neighboring processor are stored in the local solution vector and designated as “external” unknowns. Interprocessor communication occurs when an owning processor communicates the values of its border unknowns to a neighboring processor to update the value of the neighboring processor’s corresponding external unknowns. Figure 6-1 demonstrates how Processor 0 would classify the nodes in the internal, border, and external categories. Processor 0 has three neighboring processors. During the interprocessor communication phase, it sends each neighboring processor a message containing the values of each border unknown that the neighboring processor needs. The value of each border unknown may be needed by more than one processor, as it may appear in the external node lists of more than one of the neighboring processors. Processor 0 also receives a message from each of its surrounding processors containing the values of its external unknowns. Processor 0 doesn’t have to know about unknowns defined at elemental nodes which don’t have the , , or symbols attached to them. On each processor, a solution vector is stored which corresponds to the internal, border, and external unknowns defined on that processor. The solution vector is reordered locally so that local internal unknowns appear first, border unknowns appear second, and external unknowns, grouped by the owning neighboring processor, appear last. A local-to-global mapping vector is maintained, so that the global solution vector may be regenerated using “fan-in” operations. This local reordering scheme minimizes the gather/scatter operations involved in the interprocessor communication step. Only a gather operation at the originating processor to gather all of the border unknowns needed by a single neighboring processor into a contiguous space in memory is required. This message can then be directly sent to the contiguous space in the destination processor’s solution vector corresponding to the external unknowns owned by the originating processor. No scatter operations are needed on the destination processor. Moreover, the communications stencil required for this operation may be calculated once and used over and over again for a static mesh discretization. The communications stencil refers to the content of the message that each processor needs to send to each of its neighboring processors and the length of the return message containing the external unknown values from each neighboring processor.

39

Proc 2

Proc 3

Proc 1 Proc 0

Figure 6-1: Division of the nodes of an element amongst the processors, and the further differentiation of the nodes into interior ( ), border ( ), and external ( ) categories on Processor 0.

MPSalsa stores the Jacobian matrix in a distributed version of the Variable Block Row (VBR) sparse matrix format [12]. Each processor is responsible for storing rows of the Jacobian corresponding to its unknowns. Once a specific partition and assignment of the unknowns to internal, border, and external sets has been defined and the local solution vector has been reordered a distributed VBR sparse matrix is constructed. Each row of the Jacobian may include column entries corresponding to internal, border, and external unknowns defined on that processor. During the matrix-vector multiply kernel of the Krylov subspace iterative methods, each processor is responsible for carrying this out for its rows. This necessitates an interprocessor communication step wherein all external entries in the vector are updated with values from the neighboring unknowns, before the start of the operation. Calculation of matrix-vector products on rows corresponding to the internal unknowns requires no external node values and can therefore proceed simultaneously with the communication step. Much of MPSalsa’s parallel implementation is designed with the goal of maximizing the speed of this matrix-vector multiplication, which essentially requires minimizing the time needed to perform the communications. This subsection has described several strategies employed by MPSalsa to achieve rapid interprocessor communications: reordering of the solution vector to minimize work involved with the communications step, the pre-setup of the communications stencil,

40

and the ability to do calculations during the communications step. The other basic algorithmic aspect of highly efficient unstructured communication is the partitioning of the FE mesh in a way that reduces the total communication volume and message start-ups while achieving load balance over all of the processors. To do this, MPSalsa currently uses a static partitioning generated by Chaco [42], a general graph partitioning code that was developed in conjunction with MPSalsa. Chaco supports a variety of new and established graph partitioning heuristics, such as spectral techniques, geometric methods, multilevel algorithms and the Kernighan-Lin method. All of these approaches may be applied in bisection, quadrisection, or octasection mode to recursively partition general graphs for mapping onto hypercube and mesh architectures of arbitrary size. Using these techniques, a problem mapping with low communications volume, good load balancing, minimum message start-ups and small amounts of congestion can be generated.

6.2

Numerical Properties

The system of transport-reaction equations, Eqn.’s 50-58, is a system of nonlinear non-selfadjoint PDEs. After applying the Petrov-Galerkin approximation to these equations and using the Newton-Kantoravich linearization, the final matrix problem is given by Eqn. 82. In this form, we have kept the explicit time derivative term to make the following discussion of the time integration methods clearer. As described by Eqns. 79 - 82 and the definitions of the block matrices in Appendix C, these discrete equations form a nonsymmetric system of stiff Differential Algebraic Equations (DAEs). The nonsymmetric global matrix operator is a result of the convection operators in the transport part of the equations, and the stiffness in the equations is the result of the disparate time scales for the fast chemical kinetics terms and the relatively slow transport processes of diffusion and convection. The stiffness and the strongly coupled nature of the reaction operators, combined with the elliptic behavior of the pressure for incompressible flows, lead to a natural choice of fully implicit time integration techniques to provide stable time integration. The nonsymmetric character of these equations requires the use of nonsymmetric iterative methods.

6.3

Transient Solution Algorithms

The transient time integration methods used in MPSalsa follow closely the development of Gartling [21] in the NACHOS II code and the work of Gresho [36]. When appropriate, we have used the discussion from [21] with the author’s permission. Two types of implicit predictor/corrector integrators are used in MPSalsa: Forward/Backward Euler and Adams-Bashforth/Trapezoidal Rule. As discussed above, implicit solution methods are preferred for transport-reaction equations. Explicit methods suffer from a number of difficulties, including a) the strong elliptic nature of pressure in incompressible flows, b) severe time step limitations needed to maintain stability, c) fully integrated and consistent mass matrices require the inversion - defeating the efficiency of the explicit method, d) the reduction of accuracy due to diagonalizing M (ρ) to avoid (c). Effective explicit time integration demands 1-pt-quadrature and the associated stable lumped mass matrix. Though computationally expensive, implicit methods

41

are desirable because of their stability and ability to integrate efficiently to steady state solutions for problems where the diffusion operator is important. The implicit time integrators in MPSalsa are based on predictor/corrector methods to improve their accuracy and efficiency. Both integrators may be used with either a constant or dynamic time step selection algorithm. A solution of the resulting nonlinear, algebraic system for each time plane is obtained by the inexact Newton method described in Section 6.4. 6.3.1 Forward /Backward Euler Integration The first-order integration method in MPSalsa employs a forward Euler scheme as a predictor, with the backward Euler method as a corrector. The scheme uses the forward Euler predictor, V

p

n+1

= V n + ∆t n V˙ n , ∆t n = t n + 1 – t n .

(85)

The implicit backward Euler corrector uses the following approximation for the time derivative of the solution vector 1 V˙ n + 1 = -------- ( V n + 1 – V n ) , ∆t n

(86)

to solve the residual equations at t n + 1 .

R ( V˙ n + 1, V n + 1 ) = 0

(87)

In Eqn. 85 and 86, the subscript indicates the time plane index, the superscript p denotes the predicted value at time t n + 1 . The solution of the implicit corrector, Eqn. 87, at t n + 1 is obtained by the inexact Newton scheme outlined in the Section 6.4. The rate of convergence of Newton’s method is greatly increased if the initial solution estimate is “close” to the true solution. The solution predicted from Eqn. 85 provides this initial guess for the inexact Newton scheme. Appendix C provides the details of developing the discrete Newton equations for the governing transport-reaction equations. 6.3.2 Adams-Bashforth/Trapezoidal Rule Integration An explicit integration method that is the second-order analogue to the forward Euler method is the variable step, Adams-Bashforth predictor given by V

p

n+1

∆t n ∆t n ∆t n  - V˙ n – --------------- V˙ n – 1 . = V n + -------- 2 + -------------2  ∆t n – 1 ∆t n – 1 

(88)

This formula can be used to predict the solution vector, given the time derivatives at the previous two time steps, V˙ n and V˙ n – 1 . A compatible corrector equation is available in the form of the trapezoidal rule. This corrector uses an approximation to the time derivative as

42

2 V˙ n + 1 = -------- ( V n + 1 – V n ) – V˙ n . ∆t n

(89)

Eqn. 89 is then used in Eqn. 87 to find the solution at t n + 1 . Eqn. 89 is also used to calculate the time derivative at t n + 1 for later use in the predictor equation, Eqn. 88, in later time steps. 6.3.3 Time Integration Procedures The integration formulas above form the basis for the solution of time-dependent problems in MPSalsa. The similarity of the first- and second-order methods makes it possible to include both procedures in a single algorithm. The major steps in the time integration procedure are outlined here. At the beginning of each time step, it is assumed that all of the required solution and time derivative vectors are known and the time increment for the next step has been selected. To advance the solution from time t n to time t n + 1 requires the following steps: p

1) A tentative solution vector, V n + 1 , is computed using the predictor equation (either Eqn. 85 or Eqn. 88). 2) The implicit corrector equation, Eqn. 87, using Eqn. 86 or Eqn. 89 for the time derivative, is solved for the actual solution, V n + 1 . This involves the iterative solution of the linear matrix p equation arising from Newton’s method. The predicted values V n + 1 are used to initialize the FE residuals and the Jacobian matrix for the Newton iterations. 3) The time derivative vectors are updated using the new solution V n + 1 and Eqn. 86 or Eqn. 89. These equations can be conveniently described by the following relationship for the time derivative, V˙ n + 1 = CJ ( V n + 1 – V n ) – ( α – 1 )V˙ n ,

(90)

where  1  ------- ∆t n CJ =  2  ------ ∆t n 

and

1   α =   2 

order = 1 .

(91)

order = 2

The relation, dV˙ n + 1 ⁄ dV n + 1 = CJ , can be used in the formulation of the Jacobian. 4) A new integration time step is computed. The time step selection process is based on the analysis of the time truncation errors in the predictor and corrector formulas as described in the Section 6.3.4. If a constant time step is being used, this step is omitted.

43

6.3.4 Time Step Control The time integration procedures above can be used with either a user-defined constant time step or a dynamically controlled time step that is initialized with the user-defined time step size. In general, the a priori selection of a time step size can be a very difficult task, especially for stiff reacting flow equations with complex fluid flows. One of the benefits of using the predictor/corrector algorithms is that they provide a rational basis for dynamically selecting the time step size. The details of time step control algorithm can be found in Gresho et al. [36]. The general formulation of the time step selection process comes from well-established procedures for solving ordinary differential equations. By comparing the time truncation errors for two time integration methods of comparable order, a formula can be developed for predicting the next time step, based on a user-specified error tolerance. In the present case, the time truncation errors for the explicit predictor and the implicit corrector steps are analyzed and provide the required formulas. The time step estimation formula is given by [36] as m

∆t n + 1 = ∆t n ( br ε ) ,

(92)

where m = 1/2, b = 2 for the first-order method and m = 1/3, b = 3 ( 1 + ∆t n – 1 /∆t n ) for the second-order scheme. Also, r ε is a ratio of the desired time integration error to an estimate of the time integration error. Clearly, when r ε is large, a larger time step can be taken and when r ε is small, a shorter time step must be used. In practice, we have selected a measure of the time integration error that works well for the combined fluid flow and reaction kinetics problem. In MPSalsa, this ratio is computed as N

unk εˆ i 1 r ε = ----------- ∑ ----------------------, N unk p i = 1 (V i – V i)

(93)

where the subscript i refers to the component of the solution vector, N unk is the total number of unknowns and εˆ i is the desired integration accuracy for this component. For the fluid velocity unknowns, εˆ i = ε r u ∞ , where ε r is the relative accuracy desired; for temperature, εˆ i = ε r T ∞ . These measures enforce a minimum relative accuracy of time integration for the computed value locally, compared with a measure of the maximum value of the variable in the domain and are very similar to the values used in NACHOS II [21]. The hydrodynamic pressure, P , does not influence the step size control norm since there is no time derivative of the pressure in the governing transport equations. However, the determination of convergence at each time step does involve the pressure unknown. For the mass fraction unknowns, MPSalsa requires that the local time truncation error be small relative to the magnitude of the local variable and to an absolute measure of accuracy since even trace amounts of a specific chemical species can produce significant changes in the kinetics. To accomplish this, MPSalsa uses εˆ i = ε r Y k, i + ε a , where ε a is the desired absolute accuracy.

44

6.4

Inexact Newton Method with Backtracking

In this section, we briefly discuss an implementation of Newton’s method that uses approximate iterative solution techniques to solve the sequence of linear problems produce by the Newton linearization scheme. The particular implementation we use follows the work of Eisenstat and Walker [37,38,39]. This method differs from standard Newton implementations as follows. First the inexact Newton scheme uses iterative solution techniques rather than direct matrix inversion methods. Second, at each stage of the Newton iteration, the algorithm selects an appropriate level of convergence required for the iterative linear solver. This strategy is used to increase robustness of the nonlinear algorithm and to ensure that the linear equations are not over-solved at early stages of the Newton iteration when the Jacobian matrix is not very accurate. Third, this algorithm requires that at each step of the Newton iteration, the nonlinear residual must decrease. If this condition is not satisfied, a backtracking algorithm decreases the Newton step size and re-evaluates the residual at this new proposed solution. The backtracking algorithm is called recursively until the residual reduction criteria is satisfied and a new approximate solution is obtained. 6.4.1 Nonlinear Convergence Criteria Two separate convergence requirements are enforced for the Newton scheme. The first requires that the ratio of the norm of the current nonlinear residual to the norm of the initial residual –2 be reduced by a preset factor (default: 10 ). The second criterion requires that the Newton correction for any variable be suitably “small” compared to the magnitude of the variable. This criterion is very similar to the ratio used to dynamically control the time step size and is standard in general purpose ODE packages such as LSODE [44]. This convergence criterion is given by N

unk ∆V i 1 ----------- ∑ ------------------------ < 1. N unk εr V i + εa

(94)

i=1

This criterion requires the ratio of the Newton correction ∆V i be small relative to the variable V i with constant ε r , and to be small in absolute terms compared to ε a . This assures that all variables, even variables with small magnitude (e.g., trace species), are considered in determining when to halt the Newton iteration.

6.5

Linear System Solvers

The linear systems generated by the Newton iteration are iteratively solved using preconditioned Krylov methods. The methods are among the fastest and most robust iterative methods currently available. Our implementation of MPSalsa uses a parallel preconditioned Krylov solver library called Aztec [10, 11]. The Aztec library provides an efficient and well-defined interface to a number of advanced parallel iterative solution methods. These include the well-known conjugate gradient (CG) method for symmetric positive definite systems and a number of closely related algorithms for the solution of nonsymmetric systems (e.g. generalized minimum residual method (GMRES) and transpose free quasi-minimum residual method (TFQMR)) as well as various algebraic and domain decomposition preconditioners.

45

For robust and efficient solution procedures, MPSalsa and Aztec use a sparse block storage scheme called the variable block row (VBR) format [12]. Storing the matrix in a sparse format allows very efficient iterative computational kernels to be used [36, 40] and allows for the use of robust general preconditioning methods [10]. These robust schemes are critical to the solution of the strongly-coupled physics solved in MPSalsa. In the VBR format, the nonlinear dense coupling of the Jacobian at each FE node is stored intact as a small dense block. Details of the Aztec solver library can be obtained from [10, 11].

46

Appendix A: Semi-discrete Conservation Form of Transport Equations

In this section, we present a short discussion on the implementation of the governing conservation equations and how the discrete form of the conservation properties (implied by the continuous PDEs) are enforced. The heuristic explanation that follows is a generalization of the analysis from [41] for incompressible flows. For this discussion we consider the continuity equation along with a generic transport equation devoid of all dissipation and source terms: ∂ρ ------ + ∇•( ρu ) = 0 , ∂t

(A.1)

∂( ρΘ ) --------------- + ∇•( ρuΘ ) = 0 . ∂t

(A.2)

A conservative formulation will conserve the quantity ρΘ ; i.e., the time derivative of this quantity will be zero for this ideal case. Rewriting Eqn. A.2, we obtain Eqn. A.3 with the parameter, β , introduced. ∂Θ ∂ρ ρ  ------- + u • ∇Θ + βΘ  ------ + ∇•( ρu ) = 0  ∂t   ∂t 

(A.3)

Since Eqn. A.1 holds only in the continuous case and not in the weak sense, the second term of Eqn. A.3 is not zero. It can be seen that if β is taken to be zero then the advective form of the transport equation is obtained. Similarly, if β is taken to be one then the conservative form of Eqn. A.2 is obtained. Next, the Galerkin form of Eqn. A.3 is produced: ∂Θ

∂( ρ )

- Φ dΩ + ∫ ρu • ∇Θ Φ I dΩ + β ∫ Θ  ----------- + ∇•( ρu ) Φ I dΩ ∫ ρ ----- ∂t  ∂t I





= 0.

(A.4)



In this equation, we will use the following two identities from the chain rule: ∂Θ

- Φ dΩ ∫ ρ -----∂t I



=

∂ρ d ρΘΦ I dΩ – ∫ Θ ------ Φ I dΩ , and ∫ ∂t dt

∫ ρu • ∇Θ ΦI dΩ







=

∫ ∇•( ρuΘ )ΦI dΩ – ∫ Θ ∇•( ρu )ΦI dΩ . Ω



Integrating the last identity by parts, taking the mass flux to be zero on the boundary of Ω , and substituting the two resulting identities into Eqn. A.4, we obtain Eqn. A.5. d ∂ρ ρΘΦ I dΩ = ( 1 – β ) ∫ Θ  ------ + ∇•ρu Φ I dΩ + ∫ ρΘu • ∇ Φ I dΩ  ∂t  dt ∫ Ω



(A.5)



Summing this equation over all of the basis functions and using the relations

47

N

∑ ΦI I=1

N

= 1 and

∑ ∇Φ I

= 0,

I=1

we obtain the global equation, Eqn. A.6. d ∂ρ ρΘ dΩ = ( 1 – β ) ∫ Θ  ------ + ∇•ρu dΩ ∫   dt ∂t Ω

(A.6)



This implies that ρΘ is conserved globally for all time if and only if β = 1 . In our implementation, we represent the material derivative in Eqn. A.2 by the relationship in Eqn. A.3. In this way, we can easily compare the advective formulation with the conservative formulation. In addition, since we use an iterative solution method the second term of Eqn. A.3 (the continuity residual) can be moved to the right hand side and evaluated from the previous iteration. Since the magnitude of this term is small, this procedure will converge rapidly and the time step difficulties associated with the density variation in Eqn. A.2 are avoided. If the mass flux of ρΘ is not zero on the boundary, additional surface integral terms arise from surface integral boundary conditions imposed on the diffusive flux of Θ . Expressions analogous to Eqn. A.6 may be derived to yield the global balance equation.

48

Appendix B: Petrov-Galerkin Pressure Stabilization Constant

In this section we present the specific procedure for calculating the pressure stabilization parameter, ( ρτ ) , used in the pressure stabilization scheme. MPSalsa uses the Petrov-Galerkin pressure stabilization of Hughes et al. [34] and Tezduyar [35] to allow the use of equal order interpolation of velocity and pressure for incompressible flows. Below we present the pressure stabilization term proposed by Tezduyar. This formulation is a generalization of the work of Hughes et al. to finite Reynolds number flows. The stability constant for the eth global element, ( ρτ ) e , e = 1, …, N e is given by Eqn. B.1. 2

 1 ( h *e )  ρ e ------ ------------ 12 µ e ( ρτ ) e =  *  he  ρ e -----------------2 V ρe 

*

0 ≤ Re e ≤ 6 (B.1) *

6 < Re e *

Here, V is a global scaling velocity, Re e is a modified element Reynolds number based on the * global scaling velocity and the effective element length of the eth element, h e . Also, ρ e and µ e are element based quantities, equal to the average value of the density and viscosity, respectively, at the elemental nodes of the eth element. The modified element Reynolds number for the eth element is defined according to Eqn. B.2. Note that it does not conform to the normal definition of an element Reynolds number because the value of the velocity in Eqn. B.2 is globally based. Finally, V is evaluated as the L2 norm of the velocity vector evaluated at all nodes in the domain. *

* Re e

V he ρe = --------------------µe

(B.2)

The element length is defined as a length scale based on the area (volume) of the element and the equivalent area (volume) of a circle (sphere) in two (three) dimensions, as shown below. 1 * h e = 2 --- ∫ dΩ π

in 2-D

(B.3)

in 3-D

(B.4)

Ωe

6 * h e = --- ∫ dΩ π

1 --3

Ωe

The integration is over element e only.

49

Appendix C: Derivation of the Discrete Matrix Equations

In this section, the definition of the PGFEM approximation to the governing nonlinear transport equations is presented. Before a numerical solution to the nonlinear system can be attempted the system is linearized. In our development, we apply Newton’s method directly to the system of nonlinear PDEs before discretization. When Newton’s method is applied directly to the PDE system rather than to its discrete approximation, it is often called the Newton-Kantoravich (NK) method [43]. This method produces a sequence of linear systems whose solution converges to the nonlinear solution. Discretization of the Newton-Kantoravich equations then produces the discrete Newton iteration system.

C.1 Obtaining the Newton-Kantoravich Equations The Newton-Kantoravich (NK) iteration for a system of equations can be derived as follows. Given a system of equations defined as F ( v∗ ) = 0 ,

(C.1)

where v* is the actual solution to the system, solve the following sequence of equations iteratively * until v = v . Given v , the present approximation to v* , solve the NK equation, F v ( v )v' = – F ( v ) ,

(C.2)

for the update v' . The term F v ( v ) is the Frechet derivative of F , which is the operator that linearizes the equations around the solution v with respect to perturbations in the direction v' : ∂ [ F ( v + εv' ) ] . ∂ ε→0 ε

F v ( v )v' ≡ lim

(C.3)

The new approximation to the solution v* is v , and is calculated from v = v + v' . The Newton-Kantoravich equations are developed by first defining the system of residual equations, Eqn. C.1, and then carrying out the computation of the left hand side of Eqn. C.2 as described by Eqn. C.3. The residual equations are derived by using the governing transport equations (Section 3), the conservative formulation (Appendix A), the Petrov-Galerkin pressure stabilization (Appendix B), and integration by parts. These equations are shown in Section C.2. Section C.3 describes the procedure for calculating the continuum NK system for each equation. Section C.4 contains the details of each matrix term in the discretized NK system.

50

C.2 Residual Equations The residuals in Eqn. C.1 are formulated as a variational problem with weighting function Φ I , which will later be associated with the finite element basis function at global node I when we discretized the equations. The second order terms have been reduced using integration by parts (Green’s Theorem) resulting in boundary integral terms.

Residual for the three components of the velocity unknown, based on the momentum equation: F

(u)

( u, P, T , Y ) =

∂u

- Φ dΩ + ∫ ρu • ∇u Φ I dΩ + ∫ ∇Φ I • T dΩ – ∫ n • T Φ I dΓ ∫ ρ ----∂t I Γ







(C.4)

∂ρ – ∫ ρgΦ I dΩ + β ∫ u ------ + ∇•( ρu ) Φ I dΩ ∂t Ω



Residual for the hydrodynamic pressure unknown based on the total mixture continuity equation modified by a pressure stabilization term: F

(P)

( u, P, T , Y ) =





∂ρ ------ + ∇•( ρu ) Φ I dΩ ∂t

(C.5)

∂u + ρτ ∫ ∇Φ I • ρ ------ + ρ ( u • ∇u ) – ( ∇•T ) – ρg dΩ ∂t Ω

Residual for the temperature unknown based on the Cˆ p -T formulation of the energy equation: F

(T )

( u, P, T , Y ) =

∂T

-Φ dΩ + ∫ Cˆ P ρu • ∇T Φ I dΩ + ∫ λ ( ∇T • ∇Φ I ) dΩ ∫ Cˆ P ρ ----∂t I Ω





(C.6)

Ng

– ∫ ( φ + Q˙ )Φ I dΩ + Ω

∑ ∫ ( W k hˆ k ω˙k + jk • Cˆ P, k ∇T )ΦI dΩ

k = 1Ω

∂ρ + ∫ n • q c Φ I dΓ + β ∫ ( Cˆ P T ) ------ + ∇•( ρu ) Φ I dΩ ∂t Γ



51

Residual equation for the kth mass fraction based on the species continuity equation for the kth mass fraction: (Y )

F k ( u, P, T , Y ) =

∂Y k ˆ k ( ∇Y • ∇Φ ) dΩ - Φ dΩ + ∫ ρu • ∇Y k Φ I dΩ + ∫ ρD ρ k I ∫ -------∂t I





(C.7)



T

Dk + ∫ ------- ( ∇T • ∇Φ I ) dΩ + ∫ ( j k • n )Φ I dΓ T Ω

Γ

∂ρ – ∫ W k ω˙k Φ I dΩ + β ∫ Y k ------ + ∇•( ρu ) Φ I dΩ ∂t Ω



C.3 The Continuum Newton-Kantoravich Equations The Newton-Kantoravich equations are derived by applying Eqn. C.3 to the system of residuals, Eqn. C.4-C.7. The solution v is, for our system, the set of unknowns: v = ( u, P, T , Y ) In this development, the mass conservation terms, i.e., those terms with the β in them, are not represented since they are assumed to be negligible in comparison with the other terms. The momentum equation derivation is done in significant detail, whereas the remaining equations are only summarized. After each equation is presented, a summary of the assumptions used in its derivation is presented.

C.3.1 Momentum Equation Eqn. C.4 is the starting point for the linearization of the momentum equation. In the following derivation of the NK equations we have often assumed that the physical properties are independent of the solution. (Notable exceptions to this are the inclusion of the dependence of the density on the temperature and mass fractions in the body force term of the momentum balance and the dependence of the reaction kinetics on temperature and mass fractions.) This assumption affects only the accuracy of the linearization. This approximation will slow the convergence of the NK scheme but does not decrease the accuracy of the converged solution generated by the NK iteration since the dependence of the physical properties on the solution is included in the residual calculation. If it is not severe, the effect of fluid property variation is handled by a successive substitution approach and through pseudo-transient time integration methods. The time derivative in Eqn. C.4 may be expanded according to Section 6.3.3’s definition. ∂u

- Φ dΩ ∫ ρ ----∂t I



= CJ ∫ ρ ( u – u old )Φ I dΩ – ( α – 1 ) ∫ ρu˙ old Φ I dΩ Ω

(C.8)



Here, CJ and α are defined in Section 6.3.3, and u old and u˙ old denote the solution vector and its time derivative at the previous time step. Eqn. C.3 is applied to the right hand side of Eqn. C.8 to

52

obtain Eqn. C.9. Lastly, u old and u˙ old are constants in this procedure, and thus, do not produce additional terms in the linearization. Note that u' is the update difference for the velocity at the current time. The density in Eqn. C.9 is also evaluated at the current time. Note that we have neglected to include the property variations of the density in this term, ∂ ∂( u + εu' ) ρ ------------------------Φ dΩ = CJ ∫ ρu'Φ I dΩ ∫ ∂ ε ∂t I ε→0 lim

(C.9)





The contribution to the linearized momentum equation from the nonlinear convection term is computed in Eqn. C.10. Property variations of the density are also ignored here. Two linear terms arise out of the nonlinear convection operator. ∂ ρ ( u + εu' ) • ∇( u + εu' ) Φ I dΩ = ε → 0 ∂ε ∫ lim



∫ ρu • ∇u' ΦI dΩ + ∫ ρu' • ∇u ΦI dΩ



(C.10)



The stress tensor volume integral term, already a linear operator, is computed as Eqn. C.11; the property variation of the viscosity is neglected. 2 ∂ T ( P + εP ', u + εu' ) • ∇Φ I dΩ = – ∫ P ' ∇Φ I dΩ – --- ∫ µ ( ∇•u' ) ∇Φ I dΩ ∫ 3 ∂ ε ε→0 lim







(C.11)

+ ∫ µ [ ∇u' + ∇u' ] • ∇Φ I dΩ T



In contrast to the time derivative and convection operator, property variations of the density are often taken into account in the body force term in accordance with the Boussinesq approximation or the ideal gas law. Variations of the density with temperature and mass fraction are taken into account in Eqn. C.12. ∂ lim ∫ ρ ( T + εT ', Y k + εY ' k )gΦ I dΩ = ∂ ε→0 ε Ω

∂ρ -T 'gΦ I dΩ + ∫ ----∂T



Ng

∂ρ

- Y ' gΦ dΩ (C.12) ∑ ∫ -------∂Y k k I

k = 1Ω

The surface traction term, n • T , is assumed to be linear and is represented in terms of the stress tensor in Eqn. C.13 to make the natural boundary conditions more apparent. ∂ [ n • T ]Φ I dΓ = ε → 0 ∂ε ∫ lim

Γ

∫ [ n • T' ]ΦI dΓ

(C.13)

Γ

2 T = – ∫ P ' – --- ( ∇•u' ) nΦ I dΓ + ∫ µn • [ ∇u' + ∇u' ] Φ I dΓ 3 Γ

Γ

Combination of Eqn. C.9-C.13 produces the left side of the continuous NK equations for the momentum transport equation, Eqn. C.14. Eqn. C.14 is a linear, continuous system, whose succes-

53

sive solution converges to the solution of the non-linear residual equation. The right side of Eqn. C.2 is given by the momentum residual equation, Eqn. C.4. CJ ∫ ρu'Φ I dΩ + ∫ ρu • ∇u' Φ I dΩ + ∫ ρu' • ∇u Φ I dΩ – ∫ P ' ∇Φ I dΩ Ω





(C.14)



∂ρ 2 T – --- ∫ µ ( ∇•u' ) ∇Φ I dΩ + ∫ µ [ ∇u' + ∇u' ] • ∇Φ I dΩ – ∫ ------T 'gΦ I dΩ ∂T 3 Ng









∂ρ

- Y ' gΦ dΩ – ∫ [ n • T ' ]Φ I dΓ ∑ ∫ -------∂Y k k I

k=1 Ω

Γ

C.3.2 Mixture Total Continuity Equation Eqn. C.5 is the starting point for this section. The development of the individual terms in the mixture total continuity NK equation is reasonably straightforward. The only exceptions to this are the time derivative term for the density and the divergence term for the stress tensor, which appears in the pressure stabilization term. These terms are discussed in detail below. The time variation of density term will be developed in significant generality; however, in the current MPSalsa implementation, this term is not included in the NK equations. ∂ρ - Φ dΩ = ∫ ----∂t I







∂ρ -----∂P

T, Y k

∂P ∂ρ ------ + -----∂t ∂T

P, Y k

∂T ------ + ∂t

Ng

∑ k=1

∂ρ --------∂Y k

T , P, Y j ≠ k

∂Y k --------- Φ I dΩ ∂t

(C.15)

In Equation C.15, the pressure dependent term is neglected due to the assumption of low Mach number flow. The remaining temperature and mass fraction terms are omitted from the NK equations, as they also are in the time derivative term and convection operator in the momentum equations. The divergence of the stress tensor is calculated as follows. For convenience, we first represent the stress tensor as a sum of the pressure and the shear-stress tensor, T = – PI + ϒ , where the shear-stress tensor is defined as Eqn. C.16. 2 T ϒ = – --- µ ( ∇•u )I + µ [ ∇u + ∇u ] 3

(C.16)

Then, the divergence of the stress tensor can be represented as ∇•T = – ∇P + ∇•ϒ .

(C.17)

Presently, with the assumption of constant properties in the mixture continuity equation and the consistent assumptions of only density variation in the momentum residual term, the NK left hand side of the mixture continuity equation is Eqn. C.18. All but the first term arise from the pressure stabilization term in the residual equation, Eqn. C.5.

54

∫ ∇•ρu' ΦI dΩ + ρτ ∫ ∇ΦI • [ CJρu' + ρ ( u • ∇u' ) + ρ ( u' • ∇u ) ] dΩ



(C.18)



∂ρ – ρτ ∫ ∇Φ I • ------T ' + ∂T Ω

Ng

∂ρ

-Y ' ∑ -------∂Y k k

k=1

g dΩ + ρτ ∫ ∇Φ I • ∇P ' dΩ – ρτ ∫ ∇Φ I • ( ∇•ϒ' ) dΩ Ω



As described above, we have assumed that a number of fluid properties are constant for the development of the NK equations. Additionally, the current version of MPSalsa assumes that the last term involving the divergence of the shear-stress tensor is negligible relative to the other terms. Eqn. C.19 expands the divergence of the shear-stress tensor. 1 T 2 ∇•ϒ ' = ∇µ • – --- ( ∇•u' )I + [ ∇u' + ∇u' ] + --- µ ∇( ∇•u' ) + µ ∇2u' 3 3

(C.19)

As can be seen, this term is composed of second-order derivatives of the velocity vector and a term that involves the spatial variation of the viscosity. In dropping the ∇•ϒ' term from Eqn. C.18, we have assumed that the property variation contribution and the magnitude of the second derivative terms are small. Clearly, these second derivative terms will be zero for a linear interpolation of the velocity. In the case of quadratic interpolation, we follow the assumption of Hughes et al. [34] and neglect these terms for the present time.

C.3.3 Thermal Energy Transport Equation The starting point for the NK thermal transport equation is Eqn. C.6. The major assumptions made in this development relate to the dropping the viscous dissipation term, φ . In addition, the enthalpy diffusional transport term is also omitted since, in most applications, it is expected to be small relative to the thermal energy source term due to chemical reactions or the enthalpy convective transport term. The first term of the third line of Eqn. 20, which represents the heat flux due to species diffusion, has not, to date, been included in MPSalsa but will be in the future since this term can be significant in certain applications. CJ ∫ ρCˆ p T 'Φ I dΩ + ∫ ρCˆ p u • ∇T ' Φ I dΩ + ∫ ρCˆ p u' • ∇T Φ I dΩ + ∫ λ ( ∇T ' • ∇Φ I ) dΩ Ng

+



˙k ˆ ∂ω -------W h ∫ k k ∂T-ΦI T ' dΩ +

k = 1Ω Ng

+











Ng

Ng

∑ ∑

˙k ∂ω ˆ -------W h I Y ' i dΩ ∫ k k ∂Y -Φ i

k = 1 i = 1Ω

∂Q˙ ˆ P, k ∇T ' Φ dΩ – -----j • C k I ∫ ∫ ∂T-ΦI T ' dΩ –

i = 1Ω

(C.20)



Ng

∂Q˙

--------- Y ' k Φ I dΩ ∑ ∫ ∂Y k

k = 1Ω

55

+∫

Γ

∂q c -------- • n T 'Φ I dΓ + ∂T

Ng

∑∫

k = 1Γ

∂q --------c- • n Y ' k Φ I dΓ ∂Y k

In Eqn. C.20, the parametric dependencies on temperature and the mass fractions of the volu˙ k , are included. metric heat source term, Q˙ , and the heat source term due to chemical reactions, ω The last line of Eqn. C.20 indicates that the parametric dependence of the surface integral boundary conditions on temperature and mass fractions is included. However, the parametric dependencies of the density, ρ , and the thermal conductivity, λ , are neglected. Also, the temperature dependence of the partial specific heat, Cˆ p , and the partial specific enthalpy, hˆk , is neglected. C.3.4 Species Transport Equations The development of the species transport NK equations follows from the residual equation, Eqn. C.7. The linearization of the reaction rate source terms are produced by numerical differentiation. The NK left hand side of Eqn. C.2 for the species transport equations is shown in Eqn. C.21. The right hand side is Eqn. C.7. ˆ k ( ∇Y ' • ∇Φ ) dΩ CJ ∫ ρY ' k Φ I dΩ + ∫ ρu • ∇Y ' k Φ I dΩ + ∫ ρu' • ∇Y k Φ I dΩ + ∫ ρD k I

(C.21)









T

T

Dk Dk + ∫ ------- ( ∇T ' • ∇Φ I ) dΩ – ∫ ------2- ( ∇T • ∇Φ I )T ' dΩ ΩT Ω T N

g ˙k ∂ω – W k ∫ ---------Φ I T ' dΩ – W k ∑ ∂T

+∫

Γ





i = 1Ω

∂j -------k • n T 'Φ I dΓ + ∂T

Ng

∑∫

i = 1Γ

˙ ∂ω --------k-Φ I Y ' i dΩ ∂Y i ∂j -------k- • n Y ' i Φ I dΓ ∂Y i

Currently, contributions from property variations of the density and variations in the regular ˆ k and D T and thermal diffusion coefficients are neglected. This lack of property variations for D k has been shown to have a deleterious effect on the convergence of MPSalsa for non-dilute multicomponent gas-phase diffusion problems, especially when the correction velocity formalism is used. These additional terms may soon be added to improve the robustness of the solution procedure for non-dilute systems. The parametric dependence of the volumetric chemical reaction rate, ˙ k , on the temperature and mass fractions is included as is the parametric dependence of the surω face integral term on the temperature and the mass fractions. Therefore, strong coupling between the mass conservation equations and the temperature equation through volumetric and surface chemical reaction source terms may be handled in a robust fashion.

56

C.4 Discrete Newton-Kantoravich Equations Using the expansions derived in Section C.3, the following discrete equations are obtained for the Newton-Kantoravich equations at the I th finite element node. These equations represent the discrete form of the Newton-Kantoravich equations, i.e., the Jacobian formulation of Newton’s method (also known as the Newton-Raphson method). The unknowns are discretized using the expansions in Eqn. 74. The discrete set of residual equations is formed by forcing the variational equations presented in Section C.2 to be satisfied for each basis function in the domain. In the equations that follow, we use lower case indices to denote the spatial index of vectors and upper case indices to indicate the index of the finite element expansion. We also use a “group finite element representation” for some of the combined physical-property/solution-unknown products. In this representation, groups of terms are evaluated only at global nodes, not quadrature points, and their spatial dependence is handled with the elemental basis functions. For example, the density is always evaluated at the finite element nodes. Then, its value at a point x is determined through interpolation with the finite element basis functions as in Eqn. C.22. The “~” oversymbol denotes this representation. N

ρ˜ (x) =

∑ ρ J Φ J ( x) .

(C.22)

J=1

After the discrete Jacobian equations are described, the discrete residual equations are presented. The two sets of equations may then be compared to show what has been left out of the Jacobian representations. Discrete Terms in the Newton Formulation of the Momentum Transport Equation (mth component) at the Ith global node The starting point for development of the discrete momentum Jacobian is Eqn. C.14, the linearized expression for the momentum equations. The index m refers to the velocity unknown in the x m -direction. Thus, u' m, J is the update of the velocity unknown corresponding to the x m -direction at the Jth global finite element node. Though the primed variables in the expressions below are not actually part of the Jacobian expression, they denote the column indices where the Jacobian entry is located. Dependent variables denoted by overbars represent variables defined in the previous Newton iterations. Since they are dependent variables, their values within an element are interpolated via the basis functions as in Eqn. C.23. N

um =



u m , J Φ J ( x)

and

J=1

∂u --------m- = ∂x l

N

∑ J=1

∂Φ u m, J ---------J-(x) , m, l = 1, 2, 3 ∂x l

(C.23)

Eqn. C.24 represents the Jacobian entries. N

N

∑ J=1

CJ ∫ ρ˜ Φ J Φ I dΩ u' m, J + Ω

∑ J=1

N

 d ∂Φ  ˜ρ  ∑ u l ---------J-  Φ I dΩ u' m, J ∫ ∂x l  Ω l = 1

(C.24)

57

Nd

N

+

∑ ∑

∂u m ˜ Φ J ---------Φ ρ ∫ ∂xl I dΩ u'l, J –

N



J = 1l = 1 Ω Nd

N



2 ∂Φ J ∂Φ I

- --------- dΩ ∑ ∑ ∫ --3- µ˜ --------∂x l ∂x m

J = 1l = 1 Ω

∑ J=1 N





u' l, J

N

N

+



N

N d  d ∂Φ J ∂Φ I  - ---------  dΩ u' m, J + ∑ ∑ ∫ µ˜  ∑ --------∂x l ∂x l  Ω l=1 J = 1l = 1

∂ρ J  -g Φ Φ dΩ T ' J – ∫ ------∂T m J I

J=1 Ω N

J=1

Nd

∑ ∑

N





J=1 Γ



Ng

u' l, J

∂ρ J

J = 1k = 1 Ω

∂( f T • n ) m -Φ J Φ I dΓ u' l, J – ∫ ---------------------∂u l

∂( f T • n ) m ---------------------∫ ∂T -ΦJ ΦI dΓ T 'J –

∂Φ J ∂Φ I

- --------- dΩ ∫ µ˜ --------∂x m ∂x l

- g  Φ Φ dΩ ∑ ∑ ∫ -------∂Y k m  J I

J = 1l = 1 Γ N

∂Φ I - dΩ P ' J Φ ∫Ω J -------∂x m

N

Y ' k, J

∂( f T • n ) m

N

-Φ J Φ I dΓ ∑ ∫ ---------------------∂P

J=1 Γ Ng

∂( f T • n ) m

-Φ J Φ I dΓ ∑ ∑ ∫ ---------------------∂Y k

J = 1k = 1 Γ

P 'J

Y ' k, J

Eqn. C.25 is the discrete residual equation for the mth component of the momentum equation at the Ith global node. The overbar symbol has been dropped from the dependent variables that are evaluated from Eqn. C.23 for the sake of clarity. N

N

∂u ∂Φ  d ∂u  ∂Φ  d ∂u  ˜ρ --------m- Φ I dΩ + ρ˜  ∑ u l --------m-  Φ I dΩ – P ---------I dΩ – 2--- µ˜  ∑ -------l  ---------I dΩ ∫ ∂t ∫ 3  ∂x l  ∂x m ∫ ∂x m ∫ ∂x l  Ω Ω Ω Ω l = 1 l=1

(C.25)

Nd

∂u ∂u ∂Φ I + ∫ µ˜ ∑  --------m- + --------l-  --------- dΩ – ∫ ρ˜ g m Φ I dΩ  ∂x l ∂x m  ∂x l Ω l=1



N

d ∂u l ∂ρ˜  ∂ρ˜  + β  ∫ ------ + ∑ ρ˜ ------- + ------- u l u˜m Φ I dΩ – ∫ ( f T • n ) m Φ I dΓ  Ω ∂t l = 1  ∂ x l ∂ x l   Γ

Here, f T • n , which in general may be a function of all of the dependent and independent variables, is the user-supplied function for the normal component of the stress tensor. ( f T • n ) m is the component of that vector in the x m -direction. The density ρ˜ in the time derivative and convection operator is not handled robustly in the Jacobian. Its dependence on T and Y k will be included in the Jacobian in the future. Part of the formalism to do this, i.e., ∂ρ˜ ⁄ ∂T and ∂ρ˜ ⁄ ∂Y k , already exists.

58

If Dirichlet boundary conditions are applied for the mth component of the velocity at the Ith global node, Eqn. C.24 and C.25 are not used. Instead, Dirichlet conditions are applied as described in Section 4.2. These conditions may have Jacobian entries associated with them due to their dependence on other dependent variables. Discrete Terms in the Newton Formulation of the Total Mixture Continuity Equation at the Ith global node The starting point for development of the discrete Jacobian for the total continuity equation is Eqn. C.18, the linearized expression. N

Nd

∑ ∑

Φ J ∂ρ˜ ρ˜ ∂--------- Φ  Φ dΩ u' m, J ∫  ∂xm- + -------∂x m J  I

(C.26)

J = 1m = 1 Ω

N

+ ( ρτ ) e

Nd

∑ ∑

CJ

Ωe

J = 1m = 1 N

+ ( ρτ ) e

Nd

J = 1 l = 1 Ωe m = 1

Nd

J = 1 m = 1 Ωe

l=1

Nd

∑ ∫ ∑

J = 1 Ωe m = 1

∑ ∫

J = 1 Ωe N

– ( ρτ ) e

∂Φ J ∂Φ I --------- dΩ u' m, J ∂x m

∑ ∑ ∫ ρ˜ ∑ ul --------∂x l

N

– ( ρτ ) e

u' m, J

∂u m ∂Φ I ρ˜ Φ J --------- --------- dΩ u' l, J ∂x l ∂x m

Nd

N

+ ( ρτ ) e

Nd

∑ ∑ ∫ ∑ N

+ ( ρτ ) e

∂Φ I

- dΩ ∫ ρ˜ ΦJ -------∂x m

∂Φ ∂Φ ---------J- ---------I dΩ P ' J ∂x m ∂x m

ρJ    ∂------- Φ  ∂T  J 

Ng

∑ ∑ ∫

J = 1 k = 1 Ωe

Nd

∂Φ I  -------g ∑ m ∂xm- dΩ T 'J m=1

∂ρ J    -------- Φ  ∂Y k  J 

Nd

∂Φ I  g m ---------  dΩ Y ' k, J ∂x m  m=1



The subscript e denotes values corresponding to the eth element, e = 1, …, N e . The symbol Ω e denotes the interior of the eth element. An implicit summation over all elements in Eqn. C.26 for the stabilization terms was omitted to enhance readability. The discrete residual expression follows:

59

Nd

∫∑

Ωm = 1

ρ˜ ˜ ∂u m ∂˜ρ   ∂----- + ρ --------- + --------- u m Φ I dΩ  ∂t ∂x m ∂x m  Ne

+

∑ e=1

Nd

( ρτ ) e ∫

Ωe m = 1 Nd

Ne

+

∑ e=1



( ρτ ) e ∫



Ωe m = 1

(C.27)

∂u m ρ˜ --------- + ∂t

Nd

∂u m ∂Φ I --------- dΩ ∂x m

∑ ρ˜ ul -------∂x l

l=1

∂P ∂Φ --------- ---------I dΩ – ∂x m ∂x m

N

Ne

∂Φ I   d ( ρτ ) e ∫ ρ˜  ∑ g m ---------  dΩ ∂x m  Ω e m = 1

∑ e=1

No stress terms other than the pressure gradient are currently included in the stabilization vector. Discrete Terms in the Newton Formulation of the Energy Equation at the Ith global node. The starting point for the development of the discrete Jacobian expression for the thermal energy equation is Eqn. C.20. The specific heat capacity of the mixture, CˆP , the thermal conductivity, λ˜ , the mixture specific enthalpy, hˆk , the production rate of species k due to gas phase reactions, ˙ k , and the volumetric source term, Q˙ , are evaluated at global element nodes, and spatially interω polated with the basis functions as in Eqn. C.22. N

∑ J=1

N

CJ ∫ ρ˜ CˆP Φ J Φ I dΩ T ' J + Ω

N

+

Nd

∑ J=1

N

 d ∂Φ  ˜ρ CˆP  ∑ u l ---------J-  Φ I dΩ T ' J ∫  ∂x l  Ω l=1

∂T -Φ dΩ u' l, J + ∫ ρ˜ CˆP ΦJ -----∂x l I

∑ ∑

J = 1l = 1 Ω N

+

∑ J=1 N

+

Ng

J = 1k = 1

+

˙ k, J  W hˆ ∂ω ------------- Φ Φ dΩ Y ' k, J k , J k  ∂Y k  ∫ J I

Ng



Nd

∂Φ J

- Φ dΩ ∑ ∫ ∑ ∑ ( jk )l Cˆ P, k --------∂x l I

J = 1 Ωk = 1l = 1 N



∑ J=1

60

∑ J=1

d ∂Φ J ∂Φ I  ˜λ  --------- ---------  dΩ T ' J ∑ ∫ ∂x l ∂x l  Ω l = 1

N ˙ k, J )  ∂ ( hˆ k, J ω  g ----------------------------  ∫ Φ J Φ I dΩ T ' J  ∑ Wk ∂T J k = 1 Ω

∑ ∑ N

N

N

∂Q˙ J   --------- Φ Φ dΩ T ' J –  ∂T J  ∫ J I Ω

N

T 'J

Ng

∑ ∑ J = 1k = 1

∂Q˙ J   ------------Φ Φ dΩ Y ' k, J  ∂Y k, J  ∫ J I Ω

(C.28)

Nd

N

+

∑ ∑

T

∂f - Φ Φ dΓ u' l, J + ∫ -------∂u l J I

J = 1l = 1 Γ N

+



T

∂f - Φ Φ dΓ T ' J + ∫ -------∂T J I

J=1 Γ

N

N

T



∂f - Φ Φ dΓ P ' J ∫ -------∂P J I

Ng

∂f - Φ Φ dΓ Y ' k, J ∫ -------∂Y k J I

J=1 Γ

T

∑ ∑

J = 1k = 1 Γ

The discrete residual for the thermal energy equation is given by Eqn. C.29. ∂T -Φ dΩ + ∫ ρ˜ Cˆ p ∫ ρ˜ Cˆ p ----∂t I





Ng

+

Nd

N

 d ∂T ∂Φ I  ∂T ˜ - Φ dΩ + ∫ λ  ∑ ------- ---------  dΩ ∑ ul -----∂x l I ∂x l ∂x l  Ω l = 1 l=1

˜

˜

∫ ∑ W k ( hˆ k ω˙ k ) ΦI dΩ – ∫ ( Q˙ ) ΦI dΩ + ∫ f Γ



Ωk = 1

T

Ng

(C.29)

Φ I dΓ

N

d ∂u l ∂ρ˜  ∂ρ˜  + ∫ ∑ j k • Cˆ p, k ∇T Φ I dΩ + β  ∫ ------ + ∑ ρ˜ ------- + ------- u l Cˆ p T Φ I dΩ  Ω ∂t l = 1  ∂ x l ∂ x l   Ωk = 1

T

In the equations above, f is the user-supplied boundary condition for the normal component of the heat conduction, n • q c . For ideal gases, hˆk is not a function of Y k and therefore does not have to be included in the partial derivative in Eqn. C.28. Quantities in the Jacobian expression, Eqn. C.28, that can be factored outside the integral sign are taken out. In particular, all source terms, whether they occur in the volume of the domain or on the surface of the domain, can be taken out of the integral due to their group finite element representation. Discrete Terms in the Newton Formulation of the Species Transport Equation at the Ith global node The starting point for the development of the discrete Jacobian expression for the species continuity equations is Eqn. C.21. Eqn. C.30 is the Jacobian expression for the kth species. N

CJ

∑ ∫ ρ˜ ΦJ ΦI dΩ

J=1 Ω

N

+

Nd

∑ ∑

N

Y ' k, J +

Nd

J=1 Ω l=1

∂Y k - Φ dΩ u' l, J + ∫ ρ˜ ΦJ -------∂x l I

J = 1l = 1 Ω

∂Φ J

- Φ dΩ ∑ ∫ ρ˜ ∑ ul --------∂x l I N

Nd

∑ ∑

(C.30)

Y ' k, J

∂Φ ∂Φ ˆ k ---------J- ---------I dΩ Y ' ˜D ρ k, J ∫ ∂x l ∂x l

J = 1 l = 1Ω

61

T

N



+

J=1 N





J=1 N

˙ k, J  ∂ω W ------------- Φ Φ dΩ T ' J – k  ∂T J  ∫ J I Ω

Nd

Y

∑ ∑

+

N



+

N

J = 1i = 1

Ng

J=1 Γ

N



∑ ∑

˙ k, J  ∂ω W ------------- Φ Φ dΩ Y ' i, J k  ∂Y i, J  ∫ J I Ω

Y

Y

N

N

∂fk ---------Φ ∫ ∂P J ΦI dΓ P 'J

J=1 Γ

∂fk Φ dΓ T ' J + ∫ ---------Φ ∂T J I

Ng

∑ ∑

∂fk ---------Φ ∫ ∂ul J ΦI dΓ u'l, J +

J = 1l = 1 Γ N

T

Dk  d ∂T ∂Φ I  - ---------  dΩ T ' J ∫ --------2-ΦJ  ∑ -----∂x l ∂x l  Ω T l=1

N D k  d ∂Φ J ∂Φ I  - ---------  dΩ T ' J – ∑ ∫ ---------  ∑ --------∂x l ∂x l  Ω T l = 1 J=1

Y

∂fk ---------Φ ∫ ∂Y i J ΦI dΓ Y 'i, J

J = 1i = 1 Γ

The discrete residual for the kth species continuity equation is described by Eqn. C.31. ∂Y ˜ --------k- Φ I dΩ + ρ˜ ρ ∫ ∫ ∂t Ω



Nd

∂Y k

∑ ul -------∂x l

Φ I dΩ

(C.31)

l=1 N

T

N

D k  d ∂T ∂Φ I   d ∂Y ∂Φ  ˆ k  ∑ --------k- ---------I  dΩ + -------- ---------  dΩ + ∫ ρ˜ D ∫ T -  ∑ -----∂x l ∂x l  ∂x l ∂x l   Ω Ω l=1 l=1 N

d ∂u ∂ρ˜  ∂ρ˜  Y ˜ ˙ k ) Φ I dΩ + ∫ f k Φ I dΓ + β  ∫ ------ + ∑ ρ˜ -------l + ------- u l Y k Φ I dΩ – ∫ (W kω  Ω ∂t l = 1  ∂ x l ∂ x l   Ω Γ

Y

In the equations above, f k is the user-supplied expression for j k • n , the outward facing normal flux of species k. One of the species continuity equations at each global node, call it species m, is replaced by the requirement that the sum of the mass fractions equals one (see Eqn. 2). For this equation, Eqn. C.32 represents the Jacobian, and Eqn. C.33 represents the discrete residual. Ng

– ∑ [ 1 ]Y ' i, I

(C.32)

i=1 Ng

1–

∑ Yi i=1

62

(C.33)

C.5 Definition of Jacobian Block Matrices

Mass matrices [ M ] IJ =

∫ ΦJ ΦI dΩ

(C.34)



[ M ( ρ ) ] IJ =

∫ ρ˜ ΦJ ΦI dΩ

(C.35)



∫ ρ˜ Cˆ p ΦJ ΦI dΩ

[ M ( ρCˆ p ) ] IJ =

(C.36)



Convection matrices Nd

[ C ( U ) ] IJ =

∂Φ J

- Φ dΩ ∑ ∫ ρ˜ ul --------∂x l I

(C.37)

l = 1Ω

[C

(T )

Nd

∂Φ J

- Φ dΩ ∑ ∫ ρ˜ CˆP ul --------∂x l I

( U ) ] IJ =

(C.38)

l = 1Ω (P) [ C m ( U ) ] IJ

(Um)

[ Dl

Nd

( Y ) ] IJ =

(P) [ D l ( U ) ] IJ

(C.39)

∂u m

- Φ dΩ ∫ ρ˜ ΦJ -------∂x l I

( U ) ] IJ =

(T )

(Yk)

∂Φ I

l = 1Ω

[ D l ( T ) ] IJ =

[ Dl

∂u l

- Φ --------- dΩ ∑ ∫ ρ˜ -------∂x m J ∂x l

=



∂T

- Φ dΩ ∫ ρ˜ CˆP ΦJ -----∂x l I

(C.41)

∂Y ˜ Φ J --------k- Φ I dΩ ρ ∫ ∂x l

(C.42)





Nd

=

(C.40)

∂Φ J ∂Φ I

- --------- dΩ ∑ ∫ ρ˜ um --------∂x m ∂x l

(C.43)

m = 1Ω

63

Diffusion matrices (Um)

[ Kl

∂Φ J ∂Φ I ∂Φ J ∂Φ I 2 ∂Φ J ∂Φ I ] IJ = – ∫ --- µ˜ ---------- --------- dΩ + δ ml ∫ µ˜ ∑ ---------- --------- dΩ + ∫ µ˜ ---------- --------- dΩ (C.44) 3 ∂x l ∂x m ∂x k ∂x k ∂x m ∂x l k Ω

[K

(T )





Nd

∂Φ J ∂Φ I

- --------- dΩ ∑ ∫ λ˜ --------∂x l ∂x l

] IJ =

(C.45)

l = 1Ω

[K

(Yk)

Nd

∂Φ J ∂Φ I

- --------- dΩ ∑ ∫ ρ˜ Dˆ k --------∂x l ∂x l

] IJ =

(C.46)

l = 1Ω

[L

(P)

Nd

∂Φ J ∂Φ I

- --------- dΩ ∑ ∫ --------∂x l ∂x l

] IJ =

(C.47)

l = 1Ω

[S

(Y k)

N

( T ) ] IJ =





[E

(T )

(C.48)

N

Ng

] IJ

N

ˆ Tk ˆ Tk  d ∂Φ J ∂Φ I   d ∂T ∂Φ  D D -------  ∑ ---------- ---------  dΩ – ∫ ------2- Φ J  ∑ ------- ---------I  dΩ l = 1 ∂ x l ∂ x l  T l = 1 ∂ x l ∂ x l  Ω T

∂Φ J   d = ∑ ∫  ∑ ( j k ) l Cˆ p, k ----------  Φ I dΩ ∂x l  k = 1 Ω l = 1

(C.49)

Divergence and gradient matrices ∂Φ J

- Φ dΩ ∫ --------∂x l I

[ Q l ] IJ =



[ Q(ρ) l ] IJ = [ R ] IJ =

(C.50)

∂Φ J

- Φ dΩ ∫ ρ˜ --------∂x l I

(C.51)



∂ρ˜

- Φ Φ dΩ ∫ -----∂x l J I

(C.52)



Surface integral matrices (U )

[ G U l m ] IJ =

64

∂ ( f T • n )m - Φ J Φ I dS ∫ ----------------------∂u l

Γ

(C.53)

(Um)

[ GP

(Um)

[ GT

] IJ =

∂ ( f T • n )m - Φ J Φ I dS ∫ ----------------------∂P

(C.54)

∂ ( f T • n )m ----------------------∫ ∂T - ΦJ ΦI dS

(C.55)

∂ ( f T • n )m ----------------------∫ ∂Y k - ΦJ ΦI dS

(C.56)

Γ

] IJ =

Γ

(U )

[ G Y k m ] IJ = (T )

[ G U l ] IJ = (T )

[ G P ] IJ = (T ) [ G T ] IJ

(T )

Γ

T

∂f Φ dS ∫ ---------Φ ∂u l J I

(C.57)

Γ

T

∂f Φ dS ∫ ---------Φ ∂P J I

(C.58)

Γ

T

∂f = ∫ ---------Φ J Φ I dS ∂T

(C.59)

Γ

[ G Y l ] IJ =

T

∂f Φ dS ∫ ---------Φ ∂Y l J I

(C.60)

Γ

Y

(Y ) [ G U l k ] IJ

∂f k = ∫ ----------- Φ J Φ I dS ∂u l

(Y ) [ G P k ] IJ

∂f k = ∫ ----------- Φ J Φ I dS ∂P

(Y ) [ G T k ] IJ

∂f k = ∫ ----------- Φ J Φ I dS ∂T

(Y ) [ G Y l k ] IJ

∂f k = ∫ ----------- Φ J Φ I dS ∂Y l

(C.61)

Γ

Y

(C.62)

Γ

Y

(C.63)

Γ

Y

(C.64)

Γ

65

Volume source term matrices Ng

(T ) [ ℜ T ] IJ



=

k=1

˙ k, J ∂ω ∂hˆ k, J  ˙ k, J ------------ Φ Φ dΩ W k hˆ k, J -------------- + ω  ∂T J  ∫ I J ∂T J

(C.65)



˙ k, J ∂ω (Y ) [ ℜ T k ] IJ = W k  --------------  ∫ Φ J Φ I dΩ  ∂T J 

(C.66)

˙ k, J ∂ω (Y ) [ ℜ Y l k ] IJ = W k  --------------  ∫ Φ J Φ I dΩ  ∂Y l, J 

(C.67)

∂Q˙ J (T ) [ ℑ T ] IJ =  ---------  ∫ Φ J Φ I dΩ  ∂T J 

(C.68)

∂Q˙ J (T ) [ ℑ Y l ] IJ =  ------------  ∫ Φ J Φ I dΩ  ∂Y l , J 

(C.69)









Property variation matrices ∂ρ

-  Φ Φ dΩ ∫ ----∂T  J I

[ B T ] IJ =



[ B Y k ] IJ =

(P) [ B T ] IJ

(P) [ B Y k ] IJ

66

(C.70)

∂ρ

-  Φ Φ dΩ ∫ -------∂Y k  J I



Nd

=

(C.71)

∂ρ

∂Φ I

-  g Φ --------- dΩ ∑ ∫ ----∂T  l J ∂x l

(C.72)

l = 1Ω Nd

=



∂Φ I ∂ρ   -------dΩ g Φ ∫ ∂Y k  l J --------∂x l

l = 1Ω

(C.73)

References

1

Shadid, J. N., Moffat, H. K., Hutchinson, S. A., Hennigan, G. L., Devine, K. D., Salinger, A. G., “MPSalsa: A Finite Element Computer Program for Reacting Flow Problems. Part 2: User’s Manual”, Sandia National Laboratory Report, SAND95–XXXX, Albuquerque, NM (1995).

2

Kee, R. J., Miller, J. A., “A Structured Approach to the Computational Modeling of Chemical Kinetics and Molecular Transport in Flowing Systems,” Sandia National Laboratories Report, SAND86–8841, Albuquerque, NM (1986).

3

Kee, R. J., Rupley, F. M., and Miller, J. A., “Chemkin-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics,” Sandia National Laboratories Report, SAND89–8009, Albuquerque, NM (1989).

4

Coltrin, M. E., Kee, R. J., Rupley, F. M., “SURFACE CHEMKIN (v. 4.0) A Fortran Package for Analyzing Heterogeneous Chemical Kinetics at a Solid-Surface —Gas-Phase Interface,” Sandia National Laboratories Report, SAND90–80033, Albuquerque, NM (1990).

5

T. D. Blacker et. al., “CUBIT Mesh Generation Environment Volume 1: Users Manual”, Technical Report, Sandia National Laboratories, SAND94-1100, May 1994

6

Schoof, L. A., Yarberry, V. R., “EXODUS II: A Finite Element Data Model,” Sandia National Laboratories Technical Report, SAND92–2137, Albuquerque, NM (1994).

7

Rew, R., Davis, G., Emmerson, S., “NetCDF User’s Guide: An Interface for Data Access, Version 2.3,”, UCAR, April, 1993.

8

Shadid, J. N., Hutchinson, S. A., Moffat, H. K., Hennigan, G. L., Hendrickson, B. A., Leland, R. W., “A 65+ Gflop/s unstructured finite element simulation of chemically reacting flows on the Intel Paragon”, Proceedings of Supercomputing ‘94, Washington, DC, pp. 673–679, Nov. 14–18, 1994.

9

Hennigan, G. L., Shadid, J. N., “A Parallel Extension of the EXODUS Unstructured FEM Format,” Sandia National Laboratories Technical Report, SAND95–XXXX, Albuquerque, NM (1995).

10 Shadid, J. N., Tuminaro R. S., “A Parallel Preconditioned Krylov Solver Library: Algorithm Description. Version 1.0,” Sandia National Laboratories Report, SAND–######, Albuquerque, NM (1995). 11 Hutchinson, S. A., Shadid, J. N., Tuminaro, R. S., “Aztec User’s Guide: Version 1.0,” Sandia National Laboratories Technical Report, SAND95–XXXX, Albuquerque, NM. (1995).

67

12 S. Carney, M. Heroux and G. Li, “A proposal for a sparse BLAS toolkit”, SPARKER Working Note #2, Cray Research, Inc., Eagen, MN, (1993) 13 Coltrin, M. E., Kee, R. J., Evans, G. H., Meeks, E., Rupley, F. M., Grcar, J. F., “SPIN: A Fortran Program for Modeling One-Dimensional Rotating-Disk/Stagnation-Flow Chemical Vapor Deposition Reactors,” Sandia National Laboratories Report, SAND91–8003, Albuquerque, NM (1991). 14 Lutz, A. E., Kee, R. J., Miller, J. A., “SENKIN: A Fortran Program for Predicting Homogeneous Gas Phase Chemical Kinetics with Sensitivity Analysis,” Sandia National Laboratories Technical Report, SAND87–8248, Albuquerque, NM (1987). 15 Ethier, C.R. and Steinman, D.A., “Exact Fully 3D Navier-Stokes Solutions for Benchmarking,” Intl. Jrnl. Num. Meth. Fluids, 19, 369-375 (1994). 16 Bird, R. B., Stewart, W. E., and Lightfoot, E. N., Transport Phenomena, John Wiley (1960). 17 Kuo, K. K., Principles of Combustion, John Wiley, NY (1986). 18 Hirschfelder, J. O., Curtis, C. F., Bird, R. B., Molecular Theory of Gases and Liquids, John Wiley and Sons, Inc., New York (1954). 19 Chapman, S., Cowling, T. G., The Mathematical Theory of Non-Uniform Gases, 3rd Ed., Cambridge University Press, Cambridge (1970). 20 Dixon-Lewis, G., “Flame Structure and Flame Reaction Kinetics II. Transport Phenomena in Multicomponent Systems,” Proc. Roy. Soc. A., 307, 111–135 (1968). 21 Gartling, D. K., “NACHOS II: A Finite Element Computer Program for Incompressible Flow Problems. Part 1 – Theoretical Background,” Sandia National Laboratories Report, SAND86– 1816, Albuquerque, NM (1986). 22 Paolucci, S., 1982, “On the Filtering of Sound from the Navier-Stokes Equations,” Sandia National Laboratories Technical Report, SAND82–8257, Albuquerque, NM (1982). 23 Martinez, M. J., “Analysis of Anelastic Flow and It’s Numerical Treatment via Finite Elements, Sandia National Laboratories Technical Report,” SAND84–0320, Albuquerque, NM (1994). 24 Kee, R. J., Dixon-Lewis, G., Warnatz, J., Coltrin, M. E., Miller, J. A., “A Fortran Computer Code Package for the Evaluation of Gas-Phase Multicomponent Transport Properties,” Sandia National Laboratories Report, SAND86–8246, Albuquerque, NM (1986). 25 Williams, F. A., Combustion Theory, Benjamin/Cummings, Menlo Park, CA (1985).

68

26 Coffee, T. P., Heimerl, J. M., “Transport Algorithms for Premixed Laminar Steady-State Flames,” Combustion and Flame, 43, 273–289 (1981). 27 Heimerl, J. M., Coffee, T. P., “Transport Algorithms for Methane Flames,” Combustion Science and Technology, 34, 31–43 (1983). 28 Oran, E. S., Boris, J. P., Numerical Simulation of Reactive Flow, Elsevier, New York (1987). 29 Mitchell, R. E., Kee, R. J., “A General-Purpose Computer Code for Predicting Chemical Kinetic Behavior Behind Incident and Reflected Shocks,” Sandia National Laboratories Report, SAND82–8205, Albuquerque, NM (1982). 30 Kee, R. J., Grcar, J. F., Smooke, M. D., Miller, J. A., “A Fortran Program for Modeling Steady Laminar One-Dimensional Premixed Flames,” Sandia National Laboratories Report, SAND85–8240, Albuquerque, NM (1985). 31 Kee, R. J., Rupley, F. M., Miller, J. A., “The CHEMKIN Thermodynamics Data Base,” Sandia National Laboratories Report, SAND87–8215, Albuquerque, NM (1987). 32 Gordon, S., McBride, B. J., “Computer Program for Calculation of Complex Chemical Equilibrium Compositions, Rocket Performance, Incident and Reflected Shocks and Chapman-Jouguet Detonations,” NASA Technical Report, NASA SP–273 (1971). 33 Moffat, H. K., Glarborg, P., Kee, R.J., Grcar, J. F., Miller, J. A., “SURFACE PSR: A Fortran Program for Modeling Well-Stirred Reactors for Gas and Surface Reactions,” Sandia National Laboratories Report, SAND91–8001, Albuquerque, NM (1991). 34 Hughes, J. R., Franca, L. P., Balestra, M., “A New Finite Element Formulation for Computational Fluid Dynamics: V. Circumventing the Babuska-Brezzi Condition: A Stable Petrov-Galerkin Formulation of the Stokes Problem Accommodating Equal-order Interpolations”, Comp. Meth. App. Mech. and Eng., 59, 85–99 (1986). 35 Tezduyar, T. E., “Stabilized Finite Element Formulations for Incompressible Flow Computations”, Advances in App. Mech., 28, 1–44, (1992). 36 Gresho, P. M., Lee, R. L., Sani, R. L., “On the time dependent solution of the incompressible Navier-Stokes equations in two and three dimensions,” Recent Advances in Numerical Methods in Fluids, Vol. 1, Pineridge Press, Swansea, U. K., 27–81 (1980). 37 Eisenstat, S. C. and Walker, H. F. “Choosing the forcing terms in an inexact Newton method.” Technical Report 6/94/75, Utah State University Math. Stat. Dept., June 1994. To appear in SIAM J. Sci. Comput. 38 Eisenstat, S. C. and Walker, H. F., Globally convergent inexact Newton methods.” SIAM J. Optimization, 4:393-422, 1994

69

39 Walker, H. F., “Summary of Activities and Results at Sandia National Laboratories, June September, 1995. 40 Schunk, P. R., Shadid, J. N., “Iterative solvers in implicit finite element codes,” Sandia National Laboratories Technical Report, SAND92–1158, Albuquerque, NM (1992). 41 Lee, L, Gresho, P., Chan, S., Sani, R., and Cullen, M., “Conservation Laws for Primative Variable Formulations for the Incompressible Flow Equations using the Galerkin Finite Element Method,” Lawrence Livermore Laboratory Report, preprint number UCRL-82868, March, 1981. 42 Hendrickson, B. and Leland, R., “The Chaco User’s Guide, Version 2.0,” Sandia National Laboratories Report, SAND94-2692, Albuquerque, NM (1995). 43 Boyd, J. P., “Chebyshev and Fourier Spectral Methods,” Lecture notes in Engineering, Springer-Verlag, Berlin, Heidelberg, 1989 44 A. C. Hindmarsh, “LSODE and LSODEI: Two new Initial Value Ordinary Differential Equation Solvers”, ACM Signum Newsletter, 15, No. 4, pp 10-11, 1980

70

Nomenclature

b

th

ak

Activity of the k

A Ai

Surface area. Pre-exponential factor in computation of rate constant for reaction i .

bulk-phase species.

c k(n) Concentration of the k

th

surface species in the n

C b(n) Average molar concentration of the n

th

th

surface phase.

bulk phase (mol cm -3).

Cˆ P

Specific heat of the mixture at constant pressure.

Cˆ p, k

Specific heat at constant pressure for species k.

dk

Diffusional driving force for species k.

D kj

Multicomponent diffusion coefficient.

D kj

Binary diffusion coefficient between species k and j.

T Dk

D km

Mixture thermal diffusion coefficient for species k. Mixture-averaged diffusion coefficient.

ˆk D

Effective Fickian diffusion coefficient for use in the Jacobian (cm2 s-1).

Ei

Activation energy for reaction i .

f n • u Dirichlet boundary condition on the normal component of the velocity. f t1 • u Dirichlet boundary condition on one of the tangential components of the velocity, in the direction t1. Vector value of the surface integral boundary condition applied on the normal component of the stress tensor.

fτ Y

fk

Value of the surface integral boundary condition for the normal component of the

diffusion flux of the kth gas-phase species. External force of gravity. g G(n) Molar growth rate per unit of surface area for bulk phase n (mol cm -2 s-1). o

∆H f , j ( T o ) Heat of formation of the jth species at the reference temperature T o . h he

Mixture enthalpy per unit mass. *

Effective element length of element Ω e .

hˆ k

Specific enthalpy of species k (per unit mass).

i I

Vector [ 1, 0, 0 ] . Identity matrix or second order tensor.

T

71

T

j

Vector [ 0, 1, 0 ] .

jk

Diffusive flux for species k (gm cm-2 s-1).

k

Vector [ 0, 0, 1 ] .

T

f

Forward rate constant for the i

ki

r

ki

Reverse rate constant for the i

c

th

th

reaction. reaction.

Equilibrium constant in concentration units for reaction i.

Ki

f

th

l

th

K b (n) First bulk species in the n K b(n) Last bulk species in the n

bulk phase. bulk phase.

f

th

l

th

K s (n) First surface species in the n K s(n) Last surface species in the n th

surface phase. surface phase.

Ln

Film thickness for the n

bulk phase.

N Nd

Number of global nodes. Number of dimensions in the problem.

Ne

Number of elements in Ω .

Ng

Number of gas-phase chemical species; also the number of gas-phase species equa-

NR

tions. Number of elementary reversible or irreversible reactions.

bulk

N phase Number of bulk phases. surf

N phase Number of surface phases. N unk

Total number of solution unknowns.

P Po

Hydrodynamic pressure. Thermodynamic pressure.

q qc

Total heat flux vector. Heat conduction vector.

qr

Radiative heat flux vector.

qi

Rate-of-progress variable for the i

s˙k

*

e

Modified element Reynolds number for element Ω e . Surface production rate of gas- or surface-phase species k due to surface reactions (mol cm-2 s-1).

72

gas-phase reaction (mol cm -3 s-1).

Universal gas constant.

R Re

th

T

Shear stress tensor. Surface traction vector. Temperature (Kelvin). Time.

τ T t Q˙ u Vk Wk

Volumetric source term for the energy equation. Mass averaged velocity (cm s-1). Diffusion velocity of species k. Molecular weight of species k.

W x Xk

Mean molecular weight of mixture. A point in space; x = ( x, y ) in 2-D; x = ( x, y, z ) in 3-D. Mole fraction of species k.

b

X k (n) Bulk mole fraction for bulk species k in the n

th

bulk phase.

[ X k ] Concentration of species k (moles cm-3). Yk Mass fraction of species k. Z k(n) Surface site fraction for surface species k for the n

th

surface phase.

GREEK β

Parameter in residual of continuity equation.

β βi

Coefficient of volumetric expansion. Temperature exponent in computation of rate constants for reaction i .

ε εa

Specific internal energy of the mixture (erg gm-1). User-specified absolute accuracy.

εr

User-specified relative accuracy.

Γ Γn

Boundary of computational domain Ω . Surface site density for surface phase n.

χk

Chemical symbol for the k

ϒ

Shear stress tensor.

ρ µ λ

Mixture density (gm cm-3). Mixture dynamic viscosity. Mixture thermal conductivity.

σk

Number of surface sites covered by the k

ν ki

ν'' ki + ν' ki .

ν′ ki

Stoichiometric coefficient of the k

th

species.

th

th

species.

th

species for the forward direction of the i gas-

73

phase reaction. th

ν″ ki

Stoichiometric coefficient of the k

φ ΦJ

phase reaction. Viscous dissipation term in energy equation. Global finite element basis function at node J .

˙i ω

Volumetric molar rate of production of species i (mol cm-3 s-1).



Computational domain.

74

th

species for the reverse direction of the i gas-

Distribution EXTERNAL DISTRIBUTION: Steve Ashby Lawrence Livermore Nat. Lab. M/S L-316 PO Box 808 Livermore, CA 94551-0808

US Department of Energy OSC, ER-30, GTN Washington, DC 20585 T. Chan UCLA 405 Hilgard Ave. Los Angeles, CA 90024-7009

Rob Bisseling Department of Mathematics Budapestlaan 6, De Uithof, Utrecht PO Box 80.010, 3508 TA Utrecht The Netherlands

Warren Chernock Scientific Advisor DP-1 US Department of Energy Forestal Bldg. 4A-045 Washington, DC 20585

Petter Bjorstad University of Bergen Institutt for Informatikk Thomohlengst 55 N-5008 Bergen, Norway

Doug Cline The University of Texas System Center for High Performance Computing %Balcones Research Center 10100 Burnett Road, CMS 1.154 Austin, Texas 78758

Randall Bramley Dept of CSci Indiana University Bloomington IN 47405 Rich A. Cairncross Mechanical Engineering Department University of Delaware 313 Spencer Laboratory Newark, DE 19716-3140 G. F. Carey ASE/EM Dept., WRW 305 University of Texas Austin, TX , 78712 Steven P. Castillo Klipsch School of Electrical & Computer Eng. New Mexico State University Box 30001 Las Cruces, NM 88003-0001 J. M. Cavallini

Tom Coleman Dept. of Computer Science Upson Hall Cornell University Ithaca, NY 14853 Prof. D. S. Dandy Colorado State Univ. Dept. Agriculture and Chem. Eng. Fort Collins, CO 80523 Prof. J. Derby Dept. of Chemical Eng. and Materials Science University of Minnesota 421 Washington Ave. S.E. Minneapolis, MN 55455 J. J. Dongarra Computer Science Dept. 104 Ayres Hall University of Tennessee Knoxville, TN 37996-1301

75

I. S. Duff CSS Division Harwell Laboratory Oxfordshire, OX11 ORA United Kingdom

111 College Place Syracuse, NY 13244 R. F. Freund NRaD- Code 423 San Diego, CA 99152-5000

Erik Egan Motorola Semiconductor Products Sector 2200 West Broadway Rd., MS350 Mesa, AZ 85202

D. B. Gannon Computer Science Dept. Indiana University Bloomington, IN 47401

Alan Edelman Dept. of Mathematics MIT Cambridge, MA 02139 %[email protected]

Horst Gietl nCUBE Deutschland Hanauer Str. 85 8000 Munchen 50 Germany

Steve Elbert US Department of Energy OSC, ER-30, GTN Washington, DC 20585

Paul Giguere Group TSA-8 MS K575 Los Alamos National Laboratory Los Alamos, NM 87545

H. Elman Computer Science Dept. University of Maryland College Park, MD 20842 R. E. Ewing Mathematics Dept. University of Wyoming PO Box 3036 University Station Laramie, WY 82071 Charbel Farhat Dept. Aerospace Engineering UC Boulder Boulder, CO 80309--0429 J. E. Flaherty Computer Science Dept. Rensselaer Polytechnic Inst. Troy, NY 12180 G. C. Fox Northeast Parallel Archit. Cntr.

76

John Gilbert Xerox PARC 3333 Coyote Hill Road Palo Alto, CA 94304 R. J. Goldstein Mecnaical Enginerring Department University of Minnesota 111 Church St. Mineeapolis, MN 55455 G. H. Golub Computer Science Dept. Stanford University Stanford, CA 94305 Anne Greenbaum New York University Courant Institute 251 Mercer Street New York, NY 10012-1185

Satya Gupta Intel SSD Bldg. CO6-09, Zone 8 14924 NW Greenbrier Parkway Beaverton, OR , 97006 J. Gustafson Computer Science Dept. 236 Wilhelm Hall Iowa State University Ames, IA 50011 Michael Heath Univ. of Ill., Nat. CSA 4157 Bechman Institute 405 North Matthews Ave. Urbana, IL 61801-2300 Mike Heroux Cray Research Park 655F Lone Oak Drive Eagan, MN 55121 Dan Hitchcock US Department of Energy SCS, ER-30 GTN Washington, DC 20585 Fred Howes US Department of Energy OSC, ER-30, GTN Washington, DC 20585 Prof Marylin C. Huff Department of Chemical Engineering University of Delaware Newark, DE 19716 Prof. K. J. Jensen Massachusetts Institute of Technology Dept. Chem. Eng. MIT 66-566 Cambridge, Mass. 02139-4307 Christopher R. Johnson Department of Computer Science 3484 MEB

University of Utah Salt Lake City, UT 84112 David Keyes NASA Langley Research Center ICASE M/S 132C Hampton VA 23681-0001 David Kincaid Center for Numerical Analysis RLM 13.150 University of Texas Austin, TX , 78713--8510 T. A. Kitchens US Department of Energy OSC, ER-30, GTN Washington, DC 20585 Vipin Kumar Computer Science Department Institute of Technology 200 Union Street S.E. Minneapolis, MN 55455 Joanna Lees Intel Corp. Scalable Systems Division CO1-15 15201 NW Greenbrier Parkway Beaverton, OR 97006 John Lewis Boeing Corp. M/S 7L-21 P.O. box 24346 Seattle, WA 98124-0346 T. A. Manteuffel Department of Mathematics University of Co. at Denver Denver, CO 80202 S. F. McCormick Computer Mathematics Group

77

University of CO at Denver 1200 Larimer St. Denver, CO 80204 Robert McLay University of Texas at Austin Dept. ASE-EM Austin, TX 78712 %[email protected] P. C. Messina 158-79 Mathematics & Comp Sci. Dept. Caltech Pasadena, CA 91125 C. Moler The Mathworks 24 Prime Park Way Natick, MA 01760 Gary Montry Southwest Software 11812 Persimmon, NE Albuquerque, NM 87111 D. B. Nelson US Department of Energy OSC, ER-30, GTN Washington, DC 20585 Kwong T. Ng Klipsch School of Electrical & Computer Eng. New Mexico State University Box 30001 Las Cruces, NM 88003-0001 S. V. Patankar Mecnaical Enginerring Department University of Minnesota 111 Church St. Mineeapolis, MN 55455 Linda Petzold L-316 Lawrence Livermore Natl . Lab.

78

Livermore, CA 94550 Barry Peyton Mathematical Sciences Section Oak Ridge National Laboratory PO. Box 2008, Bldg. 6012 Oak Ridge, TN , 37831-6367 Paul Plassman Math and Computer Science Division Argonne National Lab Argonne, IL 60439 Claude Pommerell AT&T Bell Labs 600 Mountain Ave, Room 2C-548A Murray Hill, NJ , 07974--0636 Alex Pothen Department of Computer Science Old Dominion University Norfolk, VA 23529-0162 J. Rattner Intel Scientific Computers 15201 NW Greenbriar Pkwy. Beaverton, OR 97006 Patrick Riley Intel-SSD 600 S. Cherry St., Suite 700 Denver, CO 80222 Ed Rothberg Silicon Graphics, Inc. MS 7L-580 2011 N. Shoreline Blvd. Mountain View, CA 94043 Y. Saad University of Minnesota 4-192 EE/CSci Bldg. 200 Union St. Minneapolis, MN 55455-0159 Joel Saltz

Computer Science Department A.V. Williams Building University of Maryland College Park, MD 20742

Richard Sincovec Mathematical Sciences Section Oak Ridge Nat. Lab. P.O. Box 2008, Bldg. 6012 Oak Ridge, TN 37831-6367

A. H. Sameh CSRD, University of Illinois 305 Talbot Laboratory 104 S. Wright St. Urbana, IL 61801

Vineet Singh HP Labs, Bldg. 1U, MS 14 1501 Page Mill Road Palo Alto, CA 94304

P. E. Saylor Dept. of Comp. Science 222 Digital Computation Lab University of Illinois Urbana, IL 61801

Anthony Skjellum Mississippi State University Computer Science PO Drawer CS Mississippi State, MS 39762

Rob Schreiber RIACS NASA Ames Research Center Mail Stop T045-1 Moffett Field, CA 94035-1000

L. Smarr Director, Supercomputer Apps. 152 Supercomputer Applications Bldg. 605 E. Springfield Champaign, IL 61801

M. H. Schultz Department of Computer Science Yale University PO Box 2158 New Haven, CT 06520

Harold Trease Los Alamos National Lab PO Box 1666, MS F663 Los Alamos, NM 87545

Mark Seager LLNL, L-80 PO box 803 Livermore, CA 94550 Horst Simon Silicon Graphics Mail Stop 7L-580 2011 N. Shoreline Blvd. Mountain View, CA 94043 T. W. Simon Mecnaical Enginerring Department University of Minnesota 111 Church St. Mineeapolis, MN 55455

C. VanLoan Department of Computer Science Cornell University, Rm. 5146 Ithaca, NY 14853 John VanRosendale ICASE, NASA Langley Research Center MS 132C Hampton, VA 23665 Steve Vavasis Department of Computer Science / ACRI 722 Engineering and Theory Center Cornell University Ithaca, NY 14853 R. G. Voigt MS 132-C

79

NASA Langley Resch Cntr, ICASE Hampton, VA 36665 Phuong Vu Cray Research, Inc. 19607 Franz Road Houston, TX 77084 Steven J. Wallach Convex Computer Corp. 3000 Waterview Parkway PO Box 83385l Richardson, TX 75083-3851 G. W. Weigand U.S. DOE 1000 Independence Ave., SW Room 4A-043 (DP1.1) Washington, DC 20585 Olof B. Widlund Dept. Computer Science Courant Inst., NYU 251 Mercer St. New York, NY , 10012

INTERNAL DISTRIBUTION: 1 1 1 1 20 1 1 1 10 50 40 10 1 1 1 1 1

MS 0151 MS 0321 MS1427 MS0601 MS 0601 MS0601 MS 0827 MS 1111 MS 1111 MS 1111 MS 1111 MS 1111 MS 1111 MS 1111 MS 1111 MS 1111 MS 1110

80

Gerold Yonas, 9000 William Camp, 9200 P. Mattern, 1100 P. Esherick, 1126 Harry K. Moffat, 1126 M. E. Coltrin, 1126 J. S. Rottler, 5600 Sudip Dosanjh, 9221 Scott Hutchinson, 9221 John N. Shadid, 9221 Andrew G. Salinger, 9221 Gary L. Hennigan, 9221 Mark P. Sears, 9221 Daniel Barnette, 9221 Steven J. Plimpton, 9221 David R. Gardner, 9221 Richard C. Allen, 9222

1 1 1 1 1 1 1 1 10 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

MS 1110 MS 1110 MS 1110 MS 1110 MS 1109 MS 1109 MS 1109 MS 1109 MS 1109 MS 1109 MS 1109 MS 0441 MS 0441 MS 0441 MS 0819 MS 0819 MS 0819 MS 00439 MS 0841 MS 0833 MS 0841 MS 0843 MS 0834 MS 0834 MS 0834 MS 0826 MS 0826 MS 0825 MS 0825 MS 0437 MS 0827 MS 0827 MS 0827 MS 0827 MS 0827 MS 0827 MS 0827 MS 0827 MS 0827 MS 0827 MS 0834 MS 0835 MS 0835 MS 0826 MS 0750 MS 0750

Bruce A. Hendrickson, 9222 David E. Womble, 9222 Ray S. Tuminaro, 9222 Lydie Prevost, 9222 Art Hale, 9224 Ted Barragy, 9224 Bob Benner, 9224 Robert W. Leland, 9224 Karen Devine, 9224 Courtenay Vaughn, 9224 James Tomkins, 9224 S. W. Attawy, 9225 L. A. Schoof, 9225 T. J. Tauges, 9225 J. Michael McGlaun, 9231 James S. Perry, 9231 Allem C. Robinson, 9231 David R. Martinez, 9234 P. L. Hommert, 9103 Johnny H. Biffle, 9103 E. D. Gorham, 9104 A. C. Ratzel, 9112 M. R. Baer, 9112 A. S. Geller, 9112 R. R. Torczynski, 9112 W. L. Hermina, 1553 T. J. Bartel, 9153 C. C. Wong, 9154 Basil Hassan, 9155 G. D. Sjaardema Dave K. Gartling, 9111 Randy Schunk, 9111 Phil Sackinger, 9111 Mario Martinez, 9111 Mike Glass, 9111 Bob McGrath, 9111 Poly Hopkins, 9111 Jim Schutt, 9111 Melinda Sirmar, 9111 Steve Kempka, 9111 Robert B. Campbell, 9112 Roy E. Hogan Jr., 9111 Mark A. Christon, 9111 Robert J. Cochran, 9114 Greg A. Newman, 6116 David L. Alumbaugh, 6116

1 1 1 1 1 1 1 1 1

MS 9214 MS 9042 MS 9051 MS 9051 MS 9051 MS 9051 MS 9043 MS 9043 MS 9042

1 5 1 2

MS 9018 MS 0899 MS 0619 MS 0100

Juan Meza, 8117 Joseph F. Grcar, 8745 W. T. Ashurst, 8351 Alan Kerstein, 8351 Jackie Chen, 8351 H. Najm, 8351 R. J. Kee, 8745 S. K. Griffiths, 8745 Greg Evans, 8745 Central Technical Files, 8523-2 Technical Library, 4414 Print Media, 12615 Document Processing, 7613-2 For DOE/OSTI

81

View publication stats

Lihat lebih banyak...

Comentários

Copyright © 2017 DADOSPDF Inc.