Comments
Description
Transcript
= OF
N A S A CONTRACTOR REPORT ; o M *o v I = U e r/) 4 z A COUPLED COMPUTER CODE FOR THETRANSIENT THERMAL RESPONSE A N D ABLATION OF NON-CHARRING HEAT SHIELDS A N D NOSETIPS by Curl B. Moyer, Lurry W. Anderson, und Thomus J. Duhm Prepared by AEROTHERMCORPORATION Mountain View, Calif. for Langley Research Center N A T I O N AAL E R O N A U T I C AS N SD P A C AE D M I N I S T R A T I O N W A S H I N G T O N , D. C. OCTOBER 1970 5 A S A CR-1630 I.,' A COUPLED COMPUTER CODE FOR THE TRANSIENT THERMAL RESPONSE AND ABLATION O F NON-CHARRING HEAT SHIELDS AND NOSE TIPS By C a r l B. Moyer, Larry W. Anderson, and Thomas J . Dahm n c-. c - I Issued by Originator. as ReporbNo-. 69 -60 Prepared under Contract No. NAS 1-8501 by AEROTHERM CORPORATION Mountain View, Calif. for Langley Research Center NATIONAL AERONAUTICS AND SPACE ADMINISTRATION ~~ For sale by the Clearinghouse for Federal Scientific and Technical Information Springfield, Virginia 22151 CFSTI price $3.00 - ABSTRACT The report presents the analysis background and numerical details computer code for the transient thermal response and ablation of noncharring heat shields and nose-tips. The controlling sub-code is the in-depth transient heat conduction routine; it treats two-dimensional axisymmetric bodies using the explicit finite-difference method. It is coupled to a completely general thermochemistry sub-code for the surface state and ablation computation. This routine allows for full equilibrium considering a large number of species, and can also account for any number of nonequilibrium effects. The code also features automatic coupled calculations of the pressure around the body, using a variation of modified Newtonian relations coupled to Prandtl-Meyer expansions. The subsonic region pressure description includes empirical corrections necessary on very blunt bodies. In addition, the code automatically computes the boundary layer edge static enthalpy, the recovery enthalpy, and the convective transfer coefficient, employing appropriate computed property values any for edge gas. of a iii TABLE OF CONTENTS Title Section Page iii vii ix ABSTRACT LIST OF FIGURES LIST OF SYMBOLS GENERAL 1 INTRODUCTION 1.1 PROBLEM STATEMENT 1.1.1Background 1.1.2 Present Computer Code 1.1.3 Purpose of This Report 2REPORT COMPUTER PROGRAM AND OVERVIEW OF 4 3 6 6 4 IN-DEPTH TRANSIENT THERMAL CALCULATION 3.1 GENERAL REMARKS 3.2 NODAL LAYOUT AND GEOMETRY 3.2.1 General Pattern 3.2.2 Geometry 3.2.2.1 Location of Nodal Centers 3.2.2.2 Path Lengths for Conduction 3.2.2.3 Side Areas 3.2.2.4 Volume 3.2.2.5 Geometric Effects of Surface Recession 3.2.2.6 Surface Shape 3.3 INTERNAL CONDUCTION PARAMETERS 3.4 IN-DEPTH CONDUCTION SOLUTION 3.5 TEMPERATURES O F SURFACE NODES AND SURFACE POINTS IV-V COUPLING TO SURFACE THERMOCHEMISTRY SOLUTION GENERAL ASPECTS COMPUTATIONAL APPROACHES 4.2.1 Pre-calculated Table 4.2.2 Direct-Coupled Table 4.2.2.1 General Description 4.2.2.2 Some Dangers in the Direct-Coupled Table Approach 4.2.2.3 The Other Independent Variables in the Direct Coupled Table 4.2.3 Summary 4.3 THE CONVECTIVE SURFACE ENERGY BALANCE EQUATION 4.4 CONVECTIVE SURFACE ENERGY BALANCE ITERATION PROCEDURE! 4.5 AVOIDING THE SURFACE STATE SOLUTION 4.1 4.2 5 ACE CALCULATION STATE SURFACE - P-SE IV 5.1 INTRODUCTION OF 5 . 2 USES PHYSICAL PROBLEMS TREATED BY THE ACE PACKAGE 5.3 5.3.1 Closed System in Equilibrium 5.3.2 Open Systems in Equilibrium- Simplified Case Open Systems in Equilibrium- Unequal Diffusion 5.3.3 Coefficients Open Systems in Equilibrium- With Condensed 5.3.4 Removal Material Surface Phase Non-Equilibrium Effects in Open System Calculation 5.3.5 V 6 8 8 9 10 11 11 12 15 16 17 21 21 23 23 26 26 28 31 32 32 34 35 36 36 36 37 37 42 45 49 51 Title Section Page 5.3.5.1 5.4 6 Removal of Selected Species from the System 5.3.5.2 Isolated Species 5.3.5.3 Isolated Subgroups of Elements 5.3.5.4 Low Temperature Surface Equilibrium Suppression 5.3.5.5 Full Treatment of Kinetics CODING DETAILS Basic Solution Technique Restriction on Corrections Ease Species SOME 5.4.1 5.4.2 5.4.3 52 52 53 53 53 55 55 58 58 COUPLING THE TRANSFER COEFFICIENT CALCULATION TO THE SURFACE THERMOCHEMISTRY SOLUTION (111- IV COUPLING) AND TO THE INDEPTH TRANSIENT SOLUTION (111 - V COUPLING) 60 COUPLING THE FLOW-FIELD (OR BOUNDARY LAYER EDGE STATE) SOLUTION (PHASE 11)TO THE SURFACE STATE SOLUTION (I1 - IV COUPLING) AND TO THE TRANSFER COEFFICIENT SOLUTION (I1 - I11 COUPLING) 61 8 THE BOUNDARY LAYER EDGE STATE SOLUTION 8.1 THERMOCHEMICAL COMPUTATIONS 8.2 STATIC PRESSURE DISTRIBUTION 8.2.1 Subsonic Region 8.2.2 Downstream of the Sonic Point 8.2.3 Examples of the Accuracyof the Method 8.3 EVALUATION OF OTHER PROPERTIES 63 63 64 64 68 70 75 9 THE TRANSFER COEFFICIENT EVALUATION 76 OVERVIEW OF COMPUTER PROGRAM EVENTS 82 APPENDIX 87 7 10 A 92 REFERENCES vi LIST OF FIGURES Title Figure No. Page 4 1 Sketch of Problem 2 Sketch of Typical 3 Sketch of One-Dimensional 4 Sketch of Nodal 5 Illustration of Path Lengths 6 Nomenclature of Nodal Sides 7 Sketch of Surface Nodal 8 Sketch of Surface Points 9 Sketch of Surface Geometrical Relationships (exaggerated) 14 10 Sketch 15 11 Sketch of Implicit and Explicit Temperature Links in Finite Difference Solution, for Two Typical Situations 18 12 Representation 21 13 Representation of Course of Independent Surface History, One Body Point of Considered Nodal Center Thermal of 6 Layout Thermal Location for Network Ablating 8 Material 8 9 10 Box Undergoing 12 Recession 13 Resistance Surface RC Nomenclature Energy Terms During Ablation Variables in 24 Sketch Indicating Solution Space and Solution "Sheets" Stored at any Instant During Solution 28 15 Ablation Rate B' Versus Temperature for Carbon in Air 29 16 Ablating 42 17 Diagram 18 Pressure Distribution on 60° a Sphere 19 Pressure Distribution on 20 Pressure Distributionon an 21 Pressure Distribution 22 Heat 23 Overall 24 Surface Energy Balance (Subroutine SURFB) 14 Surface of Mass Nosetip Transfer Balance Showing a on MACABRE Code Flow Cone with Sonic Corner 71 72 Sphere 73 Ellipse Ellipse 74 Prediction 81 an Coefficient 65 Nomenclature 83 Chart Operations vii Flow Chart 84 , , . .,, ., ...... ..- .... .. - - - . . . ..... -. LIST OF SYMBOLS A area of side of nodal box AT a ratioof collision integrals based on a Lennard-Jones intermolecular potential (see Ref. 3 ) ft2 "_ nodalboxsideandcenterlinelengths,Figure 7, Equation (10) in constants effectively defined by Equation (22) "_ a2 b2 constants in Equation ( 2 4 ) Btu/f t2sec OR and Btu/ft2sec Bm pre-exponential factor (see Equation ( 7 4 ) ) f bl f B' for kinetic reaction m irrQ/peueCM) thermochemical defined as (BA ablation parameter -" defined as mc/peueCM,totalablationparameter "_ definedas mr II/peueCM "- "_ same as B' see a,b,c nodal capacity, Equation (14) Stanton number for heat transfer (corrected f o r "blowing", if necessary) C HO 'ki cM and OR Stanton number for heat transfer not corrected for blowing number of k atoms in molecule i Stanton number for mass transfer, corrected , if necessary for "blowing" ix Btu/OF "- SYMBOLS (continued) specific cP heat capacity at constant pressure Btu/lb OR Zi averaged heat capacity (see Equation (A-15) Btu/lb OR C specific Btu/lb OF C see a,b,c D Fick's D constant definedby Equation (62) ft 2/sec binary ft 2/sec SP heat in law diffusion Em activation F radiation Fi,F. 7 - diffusion coefficient energy view ft'/set coefficient for kinetic reaction m factor empirical factors appearing in Equation (62) F ratio of summations defined by Equation (67) f denotes general functional relationship H total enthalpy (sensible+ chemical Hr recovery h static h heat transfer difference hC hw I Btu/lb mole + u2/2) enthalpy Btu/lb coefficient enthalpy of temperature ablating enthalpy gases total number of in the system Btu/lb Btu/lb enthalpy of "_ on temperature based wall adjacent candidate chemical species indices X material at wall to gas the wall phase species Btu/f t Btu./lb Btu/lb "- OR 2sec SYMBOLS (continued) diffusional mass flux of species i lb total diffusional mass flux of element k regardless of molecular configuration K total number of chemical Ki mass fraction of species i elements in the i/ft2sec lb/ft2sec system "- I Kk total mass fraction of molecular configuration K Pi equilibrium constant (see Equation (39)) various k thermal conductivity Btu/ft sec OR kFm forward reaction rate constant tion m (see Equation(72)) L total number of candidate cies in the system L path lengths in a nodal box, Figure 5 Le Lewis L* length scale used in heat transfer analysis, defined in Equation (126) ft M Mach "- m system mi k regardless for kinetic condensed phase of;b k/lb reacvarious spe- ft number number molecular molecular weight corner and lb/lb mole weight cXiVb of species i m node m mass flow rate per unit area thermochemical effects only mC lb/lb mole "- center row number from the surfacelb/ft2sec total mass flow rate of "char" or main ablating lb/ft2sec material per unit surface area, all effects (thermochemical plus condensed phase mechanical removal) xi I element SYMBOLS (continued) flow rate of condensed phase II removed from surface m rL lb/ft2sec mechanically a for N representing species the molecular n node and center P total pi partial pressure of species i lb/f t Pr Prandtl ”- corner heat formula column -” number ”- lb/f t pressure number Btu/ft2sec flux rate of energy input to solid surface by diffusional processes in the boundary layer rate of energy at the surface qcond qrad in out ccnduction into rate of energy radiated away from the sec 0 F/Btu Btu/lb mole constant universal gas forward Btu/ft2sec surface Btu/ft2sec (17), (18) Rm net Rmax maximum of either the as defined in Figure 17 RN material Btu/ft2sec rate of energy input to the surface by radiation from the boundary layer or from outside the boundary layer, same as qrad thermal R resistance, see Equations R solid Btu/ft2sec rate of reaction m (see OR Equation (72) 1 nose nose radius % or R* ft ft radius R* contact thermal resistance, Equation(17); sonic point radius defined in Figure 17 ft2sec OR/Btu, ft r radius ft xii . . .. . .. SYMBOLS (continued) recovery factor surface location, seeAS ft Schmidt number, see Equation(A-17) entropy surface Btu/lb OR running length ft temperature OR failure or mechanical a surface species removal temperature for OR wall (surface) temperature, general term for Ts ,n 0 TW U velocity ft/sec U* velocity used in heat transfer (124)- (126) defined in Equations V nodal V gas X analysis, ft/sec ft3 volume velocity R normal amount of tion ( 4 3 ) condensed streamwise location to surface phase II defined in Equa- or coordinate ft/sec moles I I / moles gas ft mole fraction of species j m o l j/mol general various variable independent defined in variables Section 4.3 Btu/lb Y ratio defined by Equation(13) "_ Y distance ft normal to surface xiii SYMBOLS (continued) axial coordinate unequal diffusion Equation (61) ft quantity for species i, diffusion driving force see Equation (60), and (61) unequal diffusion quantity for element gardless of molecular confiquration k re- unequal diffusion driving potential quantity for elementk regardless of molecular configuration (see Equation ( 5 9 ) GREEK a heated surface absorbtivity (taken equal Eto ~ ) “k j mass fraction of element k in species j lb k/lb j 6 angle between nsrmal radians Y a mass fraction- Z fraction nent, see Equation (60) Y isentropic A denotes AS change A0 time step in finite difference solution E emissivity E error E an empirical diffusion factor exponent (see Equation (63)) rl input multiplicative safety factor step calculation, Equation(19) body centerline and surface weighting -” expo- exponent change in or surface or point location A0 during sec emittance departure from xiv zero various correlation in time SYMBOLS (continued) time sec angle between to a shock velocity vector and a line normal radians constant in Equation (128) dimensionless (90) and (96) dynamic length ratio defined by Equations viscosity lb/ft sec stoichiometric coefficient for species j of kinetic reaction m (see Equation (72)) "_ quantities defined by Equations (A-4), (A-5) , and (A-14), respectively and lb/lb mol; ib mol/lb density 1b1f.1-3 . lb/ft2sec . m = mC mass flow at boundary heat transfer convective mass transfer convectivefilm coefficient Stefan-Boltzmann layer edge film coefficient constant lb/f t2sec lb/ft2sec lb/ft2sec Btu/ft sec energy thicknesses definedby Equation (105) and (106) ft cp parameter defined by Equation(127) "_ m ' rate coefficicht temperature exponent in the for kineticreaction m (seeEquation (74)) SUBSCRIPTS indices for path lengths L (Figure 5 ) , nodal s side areas A (Figure 6) and thermal resistances R (Figure 10) "_ OR& SYMBOLS (continued) denote side lengths, Figure 7 see A , B , C , D denotes back wall, comunicating with reservoir see A , B , C , D denotes nodal column denotes nodal box denotes tofal rial; see m C center line corner amount of "char" or ablating mate- D see A , B , C , D e denotes boundary layer outer edge, or boundary layer edge gas FD denotes fail refers to material failing (being removed) from heated durface in condensed phase, e.g., mechanical removal, melting g denotes H see CH i station i,j denotes any identifiable species: atom, ion, molecule (i reserved for gas phase) int internal contribution Equation ( A - g()A, - 1 2 ) k denotes 9, index of condensed phase removed from surface flat gas disc phase index to thermal conductivity, element xvi species mechanically SYMBOLS (continued) M see CM M match m general index m node mix refers to total mixture property, Equation (A-9) mono refers to pure species, Equations (A-lO), (A-11) mono-mix refers to simplified mixture rule, Equation (A-10) N denotes N see RN n node r see r denotes recovery, radiation rad see ref refers to an Equation ( 6 3 ) res "reservoir" S denotes t total; transition point W denotes wall, i.e., heated surface 0 see point corner and node center row number index center corner and center column r?.umberindex m rE q rad in' qrad out empirical reference communicating heated with value back of wall surface C HO 112 see x1,x2; also denote "earlier" and "later" xvii molecular weight, SYMBOLS (concluded) at times 8 and O s r respectively denotes free stream condition SUPERSCRIPTS (prime) see B ' , BAr B;C ' (prime) denotes at new time 8 ' = 8 ' (prime) at the P reaction products R reaction reactants enthalpy datum TW reference + A8 enthalpy temperature- wall temperature normalized by stagnation point pressure; also see 6,F, denotes total configuration * sonic appearance of element independent of 5 molecular point * * ** x$ denotes the current estimated during an iterative solution x** denotes the new estimated aA iterative solution see 6 xviii value of the value of variable ix the variable xi during SECTION 1 GENERAL 1.1 PROBLEM STATEMENT 1.1.1 Background INTRODUCTION The computationof the thermal and dynamic historya vehicle of entering a planetary atmosphere has a number of important features, which may be crudely summarized as follows: I. Trajectory analysis 11. Inviscid flow field calculation 111. Boundary layer calculation IV. Interactions between vehicle surface and boundary layer V. In-depth response Complete calculation of vehicle response over a period of time would involve a coupled calculation of all five phases, each computation phase providing, in effect, the time dependent boundary conditions for the adjacent phases. A recent survey (Ref.1) has carefully examined the current computational state-of-the-art for each phase and for the coupling between phases. Major deficiencies were identified in Phase 111, boundary layer calculations, where general solution procedures for reacting multicomponent boundary layers did not exist, in Phase IV, where sufficiently general chemical routines did not exist,and in Phase V, where multidimensional in-depth solution programs did not exist. In addition, no 111-IV-V coupled programs existed. Programs with IV-V coupling existed, but were limited to one space dimension for the in-depth response. Efforts to upgrade the computational state-of-the-art will probably not attempt to couple together large "complete" computer routines for all of the five phases; the resulting computer code would be unwieldly, uneconomical, and unreliable. Since, in general, the coupling aspects of the predictive problem are more interesting and important than the degree of computational detail in the individual phases, a more judicious computational approach seeks to couple detailed calculations of only one or perhaps two problem to phases less detailed approximate or correlation calculations of several of the other phases. The resulting computer code would include the important coupling effects without becoming excessively large, and would still offer good detail in at least one phase of the total problem. 1.1.2 Presentcomputer Code l a t t e r approach. The work r e p o r t e d h e r e e x e m p l i f i e s t h i s The code d e v e l o p e du n d e rt h i sp r o g r a mi n c l u d e sl a r g es u b r o u t i n e sf o rt h ed e t a i l e d computationof the t r a n s i e n t i n - d e p t h t e m p e r a t u r e r e s p o n s e o f a n a b l a t i n g h e a ts h i e l d( P h a s e V in the schemeabove)and a l s oo ft h et h e r m o c h e m i c a l the e n v i r o n m e n t g a s a n d t h e h e a t e d s u r f a c e o f t h e h e a t i n t e r a c t i o n sb e t w e e n s h i e l dm a t e r i a l( P h a s e IV). The p r o g r a mf e a t u r e s more s i m p l i f i e da p p r o a c h e s 111), f o r which a con- fortheboundarylayertransportcalculations(Phase ectivetransportcoefficientapproach f l o wf i e l dc a l c u l a t i o n s( P h a s e was a d o p t e d , a n d f o r t h e i n v i s c i d 11), f o r w h i c h e m p i r i c a l c o r r e l a t i o n m o d i f i c a t i o n so fN e w t o n i a np r e s s u r ed i s t r i b u t i o n s around the a b l a t i n g body.Thesecombine tion to define the w i t ht h ei s e n t r o p i ce x p a n s i o n assump- b o u n d a r yl a y e re d g es t a t e . The completedcode a “I1 - 111 - I V - V c o u p l e d c o d e “ , may thus be termed i n which t h e r o u t i n e s f o r p h a s e s I V and V a r e r e l a t i v e l y c o m p l e t e , g e n e r a l and l a r g e , w h i l e t h o s e f o r p h a s e s fied. and v a r i o u s are used t o f i n d t h e p r e s s u r e , I1 and I11 a r e r e l a t i v e l y s m a l l a n d s i m p l i - The n a t u r eo ft h ei n d i v i d u a lp h a s er o u t i n e s may besummarizedas follows : Phase Phase V - In-Depth Response Code C h a r a c t e r i s t i c s 2-D A x i s y m m e t rF i ci n iD t ei f f e r e n c e Transient Heat Conduction Code, E x p l i c i t ,M u l t i - m a t e r i a lC a p a b i l i t y , G e n e r a l Geometry b u t some r e s t r i c t i o n s on g r i d l a y o u t make u s e on verysharpbodiesdifficult, charringmaterials, non- some t y p e so f a n i s o t r o p yc o n s i d e r e d Phase I V - Surface Response General thermochemical state routine allowsforablationof m a t e r i a li n any one any environment, con- sidersfullequilibriumplus any number of heterogeneous and homo- g e n e o u s( b u ts u r f a c ea r e as c a l e d ) kineticallycontrolledreactions; i d e n t i f i e st h e n n o c h e m i c a l l yc o n t r o l lingsurfacematerial,detectsfailing o r m e l t i n g compounds 2 I Phase - Phase I11 BoundaryLayer Code Characteristics Transfer Coefficient routine based onenergyintegraltechnique, laminaronly at present t i m e , boundary layer radiation not evaluated Phase I1 - Flow F i e l d Modified Newtonian pressure for sharp bodies, blending into empirical pressurecorrelationsforblunt b o d i e s combined with Prandtl-Meyer expansionsforsupersonicregions; isentropic expansion used for boundary l a y e re d g ea n dr e c o v e r ys t a t ec a l c u l a t i o n s ,f l o wf i e l dr a d i a t i o n notevaluated Phase I - T r a j e c t o r y Not i n c l u d e di np r o g r a m , time historiesofstagnationpressure andenthalpytobeinput 1 . 1 . 3P u r p o s eo fT h i sR e p o r t T h i sr e p o r tp r e s e n t s a d e s c r i p t i o no ft h ec o u p l e dc o d e . I t describes of t h e codeand theflowofinformationfromonecomputation p h a s et ot h eo t h e r . The r e p o r t a l s o d e s c r i b e s t h e p h y s i c a l modelsused in e a c hc o m p u t a t i o n a lp h a s e ,a n do u t l i n e si n some d e t a i l t h e n u m e r i c a l a n a l y s i s u s e di ne a c hp h a s e . The r e p o r td o e sn o ti n c l u d e a u s e r sg u i d e ;t h i sh a sb e e n p u b l i s h e ds e p a r a t e l y( R e f e r e n c e 2). t h eg e n e r a ln a t u x e 3 SECTION 2 OVERVIEW OC FOMPUTER PROGRAM AND REPORT Figure 1 illustrates the general physical problem considered by the code. Computation Physics Main stream radiation of Calculation flux to wall D e t e r m i n a t i o n of b o u n d a r yl a y e re d g e state Boundary l a y e r e d g e Account f o r b o u n d a r y l a y e r t r a n s p o r t o f mass andenergy t o and from s u r f a c e Determinechemical state at surf a c e and perform energy balance s lution c o u p l i n gt oi n - d e p t h Substrate F i g u r e 1. -In-depththermalresponse calculation SketchofProblemConsidered I nt h ec o u p l e dc o d e ,t h er o u t i n ef o rt h ei n - d e p t ht h e r m a lr e s p o n s e c a l c u l a t i o n s( P h a s e r e f e r r e d t o as MACABRE, BoundaryLayers after"MaterialAblationwithChemicallyActive i n Re-Entry". The couplingbetweenPhases t h r o u g ht h es u r f a c ee n e r g yb a l a n c ec a l c u l a t i o n ,w h i c h s e p a r a t es u b r o u t i n e , It is f u n c t i o n sa st h em a i n ,c o n t r o l l i n gr o u t i n e . V) SURFB. V and I V e f f e c t e d i s performed by a T h i se n e r g yc a l c u l a t i o nn a t u r a l l yr e q u i r e s (among o t h e rt h i n g s )i n f o r m a t i o na b o u tt h ee n t h a l p yo ft h eb o u n d a r yl a y e rg a s e s w a l l ; t h i s i n f o r m a t i o n i s c a l l e d f o r by SURFB when needed a packageof complex t h e r m o c h e m i s t r y s u b r o u t i n e s c o n t r o l l e d by s u b r o u t i n e EQUIL. I na d d i t i o n , SURFB r e q u i r e sv a l u e sf o r a convective adjacent to the and obtained from t r a n s f e rc o e f f i c i e n t ,o b t a i n e da sn e e d e df r o ms u b r o u t i n e pressure,obtained as neededfromsubroutine RUNLP, e n t h a l p y ,o b t a i n e d as neededfromsubroutine PARBL. and RUCH c o n s t i t u t e I V mentsofthecode. 4 - I11 - RUCH, forthelocal and f o r t h e l o c a l r e c o v e r y C a l l st o PARi3L , RUNLP , I1 c o u p l i n g , f i l l i n g o u t t h e c o u p l i n g r e q u i r e - S e c t i o n s 3 through 9 b e l o w d e s c r i b e t h e v a r i o u s c o m p u t a t i o n a l p h a s e s and t h ec o u p l i n gl i n k sb e t w e e n S e c t i o n 10 t h e n p r e s e n t s some o v e r a l l f l o w c h a r t s them. T h i sn e c e s s a r i l yi n v o l v e s much d e t a i l . a recapitulationofthisintroductorysectionand of thecomputercodeoperationswhichshouldillus- t r a t e more c l e a r l y a l l o f t h e i m p o r t a n t f e a t u r e s o f t h e p r o g r a m . 5 I I SECTION 3 IN-DEPTH 3.1 TRANSIENT THERMAL CALCULATION GENERAL REMARKS The in-depth transient temperature calculation serves as the controlling logic unit for the entire coupled computer code. (It also consumes most of the computing time.) This calculation is of the familiar finite difference, lumped-capacity type. The following sections describe the nodal grid employed for this calculation and the basic equations employed. 3.2 NODAL LAYOUT AND GEOMETRY 3.2.1 General Pattern 2, is imagined The geometric shape considered, illustrated in Figure to be divided up into a grid pattern in the customary manner for finite difwithin each grid"box" will be termed fernnce computation schemes. The the a "node"; this is in contrast to terminology which calls eachofcorner grid network a node. The thermal capacity of each node is imagined to be lumped at a single point within the box; this point will be termed the nodal center. It will be convenient to assume for the moment that the nodal center may be located anywhere within the nodal box. Strategies for optimizing the location of the nodal centers will be discussed later. area For convenience, the nodes are imagined to be quadrilaterals, so that the entire nodal network is an assemblage of quadrilaterals. For bookkeeping Insulated Heated Surface Surface Insulated Surface Back Wall Convective Surface h e a t shield Figure 2. Sketch of Typical Nodal Layout 6 I convenience, the network i s imagined t o b e " c o m p l e t e " , even though the shape ofthematerial may d i c t a t e t h a t some of t h e mesh boxes are empty.Both nodes(boxes orcenters)and mesh c o r n e r s are numbered i n a row andcolumn w i l l b e l a b e l e d as an m-n system,where m . d e n o t e s a row and m-n mesh scheme may b e o r i e n t e d a r b i t r a r i l y w i t h * s o t h a t i n g e n e r a l t h e rowr e s p e c tt ot h ep h y s i c a l r - Z coordinatescheme, column d e s c r i p t i o n m i g h t seem t o b e m e r e l y a d e s c r i p t i v e a r t i f i c e ; however, in exploiting the simplicities associated with low c u r v a t u r e b o d i e s a n d h e a t s h i e l d s it hasbeenassumed in constructing the program that the heated surface is a tt h et o p of t h e c o l u m n s ,a si n d i c a t e di nF i g u r e 2. Thus as s u r f a c e recessionoccurs,theboxes a t t h et o po fe a c h column s h r i n k ; t h e o t h e r b o x e s r e m a i nf i x e d .T h i sl i m i t a t i o nt h a tt h eh e a t e ds u r f a c eb el o c a t e da tt h et o p of each column i s merely a c o n v e n i e n c e p r o c e d u r e e x p l o i t e d f o r t h e s p e c i a l case of h e a t s h i e l d s t y p i c a l of f a i r l y b l u n t c o n e s . F o r h i g h c u r v a t u r e b o d i e s s u c ha ss h a r pc o n e st h i sr e s t r i c t i o n becomes i n c o n v e n i e n t ,a l t h o u g hR e f e r e n c e systemwhich d e n o t e s a column.The n g i v e sa ne x a m p l eo ft h es u c c e s s f u lu s eo ft h ec o d e 2 on .a verysharpbody. The s i d e w a l l s o f t h e n o d a l n e t w o r k a r e p r e s u m e d i n s u l a t e d and t h e back w a l l communicatesthrough a s i m p l e h e a t t r a n s f e r c o e f f i c i e n t law ( p l u s r a d i a t i o n ) t o a " r e s e r v o i r " a t Tres. Thus t h eb o u n d a r yc o n d i t i o na tt h i sf a c e doesnotinvolvethermochemistry. t o thelayout of computingthermal c o n d u c t a n c e sb e t w e e nn o d a lc e n t e r st h a tt h e mesh scheme i s n e a r l y o r t h o g o n a l , s o t h a t c o n d u c t a n c e may b e t a k e n a s c o n d u c t i v i t y times s i d e - a r e a d i v i d e d by ** l e n g t h .S e c o n d l y , i t i s presumed t h a to n l yo n em a t e r i a la b l a t e s , and t h a t o n l yo n em a t e r i a l i s exposed t ot h eh e a t e ds u r f a c e .T h i r d l y , i t i s presumed thattheprincipaldirections of t h e r m a l a n i s o t r o p y a r e a l i g n e d w i t h t h e n o d a l mesh. T h i ss i m p l i f i e sc o m p u t a t i o n sa n dr e d u c e si n p u tr e q u i r e m e n t s F . o rt h o s e a p p l i c a t i o n s i n which t h e h e a t e d s u r f a c e i n t e r s e c t s t h e p r i n c i p a l d i r e c t i o n s o fa n i s o t r o p y at "difficult"anglesthisrestrictioncould become a majorinconvenienceand more g e n e r a l schemeswouldhave tobedevisedforthoseproblems. T h r e eo t h e r" c o n v e n i e n c el i m i t a t i o n s "h a v eb e e na p p l i e d of t h e n o d a l g r i d . F i r s t , * it i s assumed f o rt h ep u r p o s e s T h a t i s , t h e n o d a l mesh scheme may b e a b o v e o r b e l o w t h e Z - a x i s l i n e , o r i e n t e d i n anygeneraldirection,and may b e " b e n t " o r s h a p e d i n a n y convenienttotheuser(subjecttothelimitationsdescribedintheparagraphsfollowingbelow). ** A s u i t a b l e more g e n e r a l c o n d u c t a n c e c a l c u l a t i n g s c h e m e c o u l d may b e manner remove t h i s r e s t r i c t i o n w i t h o u t any e s s e n t i a l change t o t h e program. 7 3.2.2 Geometry 3 . 2 . 2 . 1L o c a t i o n of NodalCenters The studysummarized i n Reference 3 showed t h a t t o o b t a i n a node d r o p p i n g schemewhichproducesmooth surface temperature histories, it w a s necessarytoavoidassociatingthermalcapacitywiththesurfacetemperature. Intermsoftheconventionalthermal d i m e n s i o n a l schemewouldlook means t h a t a oneAll t h ec a p a c i t y RC n e t w o r k , t h i s l i k eF i g u r e 3. ~ node 1 w -" " " node 2 - node 3 ~~ node 5 node 4 "'112 -+=>==+ I I - Surf Figure 3 . ofthesurface SketchofOne-DimensionalThermal RC Network node i s l o c a t e d i n t h e c e n t e r o f t h e f i r s t n o d a l z o n e , and zone i s i n t e r p o s e db e t w e e nt h es u r f a c e h a l fo ft h et h e r m a lr e s i s t a n c eo ft h i s t e m p e r a t u r ea n dt h et e m p e r a t u r eo ft h ef i r s tc a p a c i t y lump. Note t h a t f o r m m+l temperature points, in contrast to most Thus f o r an m x n n o d a ln e t w o r ks c h e m e ,t h e r e w i l l be m x n nodal nodes o r b o x e s t h e s y s t e m h a s schemes. * plus n w i t h thermac l apacity) (not associated with thermal capacity) t e m p e r a t u r e s( a s s o c i a t e d s u r f a c et e m p e r a t u r e s . as rN and Z N r t h ec o o r d i n a t e s shown i n . F i g u r e 4 may b e r e p r e - Denoting the n o d acl e n t ecr o o r d i n a t e s f o r node m, n w i t h c o r n e r c o o r d i n a t e s a s sentedin terms of t h e c o r n e r c o o r d i n a t e s by E q u a t i o n s (1) and ( 2 ) : Corner m+l n box m, n center Corner m. n Figure 4 . 1 Corner m + l , n+l Corner m, n + l S k e t c ho fN o d a lC e n t e rL o c a t i o nF o rA b l a t i n g. M a t e r i a l *The c a p a c i t y l o c a t i o n scheme d e s c r i b e d h e r e recommended i n Reference 3 8 . i s animprovementover that r - N, m, n 5, m, n r c,m,n + 'c,m+l,n + 'c,rn+l,n+l + rc,m,n+l 4 - 'c,m,n + 'cIm+l,n + 'c,m+l,n+l + 'c,m,n+l 4 * (Nodes i n t h e a b l a t i n g m a t e r i a l w i l l automaticallybecentered;the MACABRE user has the option to put back-up material nodal centers anywhere i n t h e nodal box. ) 3.2.2.2 P a t hL e n g t h sf o rC o n d u c t i o n w i l l have a g e n e r a l I n t h i s and t h e f o l l o w i n g s e c t i o n s , t h e n o d a l c e n t e r s .oc ro m p u t i nt g hermal conl o c a t i o n rN,m,n' 2 N,m,n f oi lrl u s t r a t i vpeu r p o s e F d u c t a n c e s , it w i l l b e n e c e s s a r y t o h a v e t h e r m a l p a t h l e n g t h s b e t w e e n n o d e s . it i s c o n v e n i e n t f i r s tt oc o n s i d e rp a t hs e g m e n t si n s i d ee a c h box. I n g e n e r a l t h e r e are f o u r p a t hs e g m e n t so fi n t e r e s ta s shown i n F i g u r e 5 . The programcomputes the S i n c em a t e r i a lp r o p e r t i e sa r ea s s o c i a t e dw i t he a c hn o d a lb o x l lengths L mln,B m+l, n GeneralNodal CenterPoint F i g u r e 5. I l l u s t r a t i o no fP a t hL e n g t h s and L i nt h e m d i r e c t i o n as t hde i s t a n c ebs e t w e e tnhne . d acle n t e r mlnlD m + l and m faces of t h e box.(This i s b e c a u s et h e and t h ec e n t e r so ft h e n o d a l c e n t e r i s u s u a l l y on t h e l i n e j o i n i n g t h e c e n t e r s o f t h e s e two f a c e s . ) Thus t h ep a t h s B and D h a v et h el e n g t h s * Except for the first column, i n w h i c h n o d a l p o i n t s are c e n t e r e d o n t h e f i r s t m l l and m + l , l ) u n d e r t h e a s s u m p t i o n t h a t sideofthebox(betweencorners this l i n e o f s i d e s w i l l end a t t h e body s t a g n a t i o n p o i n t . 9 rc,m+l,n 'm, 'c,m+l,n + rc,m,n+l c,m,n 2 N,m.n + r [[ r c,m,n = + =N,m,n = I + rc,m+l,n l2 - 2 N,m, n + zc,m+l,n zc,m,n 2 [[rc,m+,,n+, + Lm,n,c - 2 N,m,n rc,m,n+l 2 +I + z 'c,rn+l,n+l c,m,n+l 2 - r 1'1% i' 1'1% N , m, n N,m. n Side areas Areas of t h e s i d e s of each box are a l s o r e q u i r e d f o r t h e r m a l c o n d u c t a n c e c a l c u l a t i o n s .F o rt h es i d e sl e t t e r e d Figure 6 . 10 - d i r e c t i o np a t hl e n g t h s 'm,n,A 3.2.2.3 l2 - 'c,m+l,n+l 2 = Lrn,n,D n N,m,n 2 +I S i m i l a r l y , for t h e - rc,m+l,n+l + n, B as shown i nF i g u r e NomenclatureofNodalSides 6 , t h ea r e a s are given by e l e m e n t a r y g e o m e t r y ( F i r s t Theorem ofPappus) + r 1" + Only t h e s e two a r e a sn e e d ('c,m,n+l - zc,m,n t o be computed f o r e a c h b o x , s i n c e t h e a r e a s on t h e A and B o f a d j a c e n t othersidesofthenodalboxareidenticalwiththeareas nodes. 3.2.2.4 ,.I+ Volume The volume of t h e e l e m e n t a l box i s by t h e Second Theorem of Pappus, + zc , m + l , n + l + 3.2.2.5 r2c,m+l,n - ~,m,n+l] + Zc,m,n+l Lr c,m,n+l . G e o m e t r i cE f f e c t so fS u r f a c eR e c e s s i o n F o rn o d e sa d j a c e n tt ot h eh e a t e ds u r f a c e ,s i d e B may move due t o s u r f a c er e c e s s i o n( a b l a t i o n ) .T h i sr e c e s s i o nr e d u c e st h et h e r m a lr e s i s t a n c e b e t w e e nt h es u r f a c ep o i n ta n dt h ea d j a c e n tn o d a lp o i n t ,i n c r e a s e st r a n s v e r s e t h e r m a lr e s i s t a n c e s( a sd i s c u s s e db e l o w ) associatedwiththenodalcenterofthenodal , and r e d u c e st h et h e r m a lc a p a c i t y box a d j a c e n t t o t h e s u r f a c e . 11 s o as t o m a i n t a i n t h e The program assumes t h a t t h e s u r f a c e r e c e s s i o n o c c u r s ratios S1,SZ = s u r f a c e points 1-earlier center s u r f a c e node of c e n t e r of node Surfacem,n+l Figure 7. a s shown i nF i g u r e S k e t c ho fS u r f a c eN o d a l 7. a c22 b2 al bl Box UndergoingRecession c1 The b l i n ej o i n st h ec e n t e r so ft h e and s e r v e s t o d e f i n e t h e l o c a t i o n o f t h e s u r f a c e p o i n t m+l,n and m + l , n + la r el o c a t e da c c o r d i n g l y , a sb e f r j x e ,e x c e p tt h a tf o rt h e s en o d e s The n o d a l p o i n t o r c e n t e r S * ' A m planes The moving c o r n e r s andvolumescomputed . and t h e a r e a s m + l and m,n,A Am,n-l,C' i s a l s o presumed t o move so t h a t i s i s the nodalbox always i n t h e a r i t h m e t i c c e n t e r o f * . is Thus t h en o d a lc e n t e r " s q u e e z e d "t o w a r dt h eb a c kw a l lo ft h en o d a lb o xa sr e c e s s i o np r o c e e d s . The p a t hl e n g t h sf o rc o n d u c t i o na r ea d j u s t e da c c o r d i n g l y . Surface Shape 3.2.2.6 by F i g u r e 7 , it is m o s t c o n v e n i e n t t o As n o t e d a b o v e a n d i l l u s t r a t e d considerthesurfacepoints t h en o d a l on t h e l l e a t e d s u r f a c e a s b e i n g a l w a y s l o c a t e d o n column c e n t e r l i n e * , t h a t i s , o nt h el i n ej o i n i n gt h ec e n t e rp o i n t s m + l and m p l a n e s( " p a r a l l e l "t ot h eh e a t e ds u r f a c e ) , where m is t h e row i n d e x o f a nodalbox a t t h e s u r f a c e and m + l i s t h e l o c a l c o r n e r i n d e x o ft h e of t h eh e a t e ds u r f a c e ( s e e F i g . 7) . is .The l o c a t i o no ft h es u r f a c ep o i n t s a l lt h ei n f o r m a t i o nn e e d e da b o u tt h es u r f a c ef o rt h o s ea b l a t i o np r o b l e m sw i t h a s p e c i f i e d ,i n p u ts u r f a c er e c e s s i o nr a t ea s a b o u n d a r yc o n d i t i o n( t h ev a r i o u s ablationproblemboundaryconditionoptionsaredescribedinSection s i n c ef o rt h o s ep r o b l e m s i t i s mostconvenient s i o nh i s t o r ya l o n gt h en o d a lc e n t e rl i n e . t os p e c i f y , The n o d a lg r i d 4 below) as i n p u t , recesas i n p u tt h u ss e r v e s * E x c e p t , as n o t e d a b o v e , f o r t h e f i r s t column o f n o d a l b o x e s , i n which t h e A o fF i g u r e 6 ) o ft h en o d a l n o d a lp o i n t sa r ec e n t e r e d on t h e f i r s t s i d e ( S i d e S i s located a t thecorner m + l , l (Figure 7 ) . boxand thesurfacepoint 12 . . ... . , todefinethevariousnodal column c e n t e r l i n e s , a n d t h e r e c e s s i o n h i s t o r y f o r e a c h column t h e n d e f i n e s t h e h i s t o r y o f t h e s u r f a c e p o i n t s i n a perfectly s t r a i g h t f o r w a r d manner.However, a n o t h e ri m p o r t a n th e a t e ds u r f a c eb o u n d a r y conditionoptioninvolvesnotinputrecessionsbutvariousenergyand chemm a s s loss fromenergybala n c ec o n s i d e r a t i o n s( O p t i o n l, d i s c u s s e d i n S e c t i o n 4.3 below) T h i sa b l a t i o n c a l c u l a t i o nd o e sn o tp r o d u c er e c e s s i o n rates d i r e c t l y , o f c o u r s e ; i n s t e a d it p r o d u c e sr a t e so f mass loss from t h e s u r f a c e . I t w i l l p r o v ec o n v e n i e n t , nevertheless, to adhere to the concept of the surface point which moves a l o n g t h e column c e n t e r l i n e as r e c e s s i o np r o g r e s s e s , T o d e f i n et h es u r f a c ep o i n t motionwithcomputed mass loss r a t e s d e t e r m i n e d f r o m t h i s g e n e r a l e n e r g y balance-determinedoptionrequiresknowledge of t h e a n g l e b e t w e e n t h e l o c a l normal t o t h e s u r f a c e and t h e column c e n t e r l i n e , s i n c e c o m p u t e d mass loss may b e t r a n s l a t e d d i r e c t l y i n t o r e c e s s i o n a l o n g a normal t o t h e s u r f a c e . Rec e s s i o na l o n gt h en o r m a l may b e p r o j e c t e d i n t o t h e n o d a l column c e n t e r l i n e once t h ea n g l eb e t w e e nt h e s et w ol i n e s i s known. The d i r e c t i o n of t h e s u r f a c e normal may c o n v e n i e n t l y b e d e t e r m i n e d f r o m t h e s l o p e o f t h e s u r f a c e , t h a t is t h e s u r f a c e s h a p e , a n d o f c o u r s e t h e d i r e c t i o n of t h e column c e n t e r l i n e i s known. istry information sufficient to calculate surface . I t i s o b v i o u s l yn o ts a f et ou s et h es l o p eo ft h eh e a t e d( t o p )s u r f a c e it is o ft h es u r f a c en o d a lb o xt oo b t a i nt h es u r f a c es l o p e ,s i n c ei ng e n e r a l neitherpossiblenoralwaysdesirabletolayout a nodalgrid "conform" t o t h e " r e a l " s u r f a c e f o r t h e e n t i r e which w i l l p r o b l e mh i s t o r y .T h e r e f o r e t h e MACABRE p r o g r a m i n c l u d e s v a r i o u s s p e c i a l c u r v e - f i t s u b p r o g r a m s w h i c h compute a s u r f a c e s h a p e a t e a c h t i m e stepduringthesolutionafterexamining t h el a y o u to ft h es u r f a c ep o i n t s .R e f e r r i n gt oF i g u r e programprogramcomputes 8 , onecurve f i t sub- the Z F i g u r e 8. s l o p eo ft h es u r f a c e ,d r / d Z , b e t w e e snu r f a cp e o i n tns - 1 S k e t c ho fS u r f a c eP o i n t s a t s u r f a c ep o i n t and n, and n n and n+l. a st h ea v e r a g eo ft h es l o p e s Thus 13 Note t h a t n-1, n , and n + l a r e s u r f a c e p o i n t s , n o t n o d a l c o r n e r p o i n t s . (Problemswithonlyonenodal column h a v eo n l yo n es u r f a c ep o i n t .T h i s r e q u i r e st h ep r o g r a mt oa b a n d o nt h ea v e r a g es l o p es c h e m e and t o u s e t h e n o d a l box heated surf ace slope as the surface slope.) Othercurve fit r o u t i n e s a r e a v a i l a b l e i ns u c c e s s i v ei n t e r v a l s . choiceof whichinvolvefittedquadratics The User’s Manual g i v e sr e c o m m e n d a t i o n sf o rt h e f i t routinesfordifferent body s h a p e s . Specialformulasareusedfortheslopesatthefirst s u r f a c ep o i n t s ,h o w e v e r ,i no r d e rt oh a r m o n i z et h ea r r a y v a l u e sw i t h some b a s i c a s s u m p t i o n s r o u t i n e s .S i n c et h ef i r s t made i n t h e c o n v e c t i v e t r a n s f e r c o e f f i c i e n t column s u r f a c ep o i n t i s presumed t o b e the s t a g - is s e t t o a v e r yl a r g ep o s i t i v ev a l u e . n a t i o np o i n t ,t h es l o p eh e r e assumed t h a t t h e s u r f a c e andsecond of s u r f a c e s l o p e I t is is s p h e r i c a l b e t w e e n t h e f i r s t a n d s e c o n d p o i n t s , i s givenby so that the second point surface slope $1, 1 - y 2 2y = where Y = - 22 r 21 2 With s u r f a c e s l o p e d e t e r m i n e d , t h e s u r f a c e movement AS computed d u r i n g the time s t e p may b e p r o j e c t e d o n t o t h e n o d a l b o x c e n t e r l i n e , Figure 9 . 14 SketchofSurfaceGeometrical R e l a t i o n s h i p s( e x a g g e r a t e d ) and t h e n t h e new n o d a l volumecomputed, as i n d i c a t e d i n t h e e x a g g e r a t e d s k e t c h o f F i g u r e 9. A c t u a l s u r f a c e movement duxing a time s t e p i s l i m i t e d t o a small f r a c t i o n o f the nodal thickness to ensure agood approximationtotheconservationof column i s an e x c e p t i o n , s i n c e t h e s u r f a c e p o i n t Noteagainthatthefirst mass. is l o c a t e d on t h e l e f t c o r n e r , n o t i n t h e c e n t e r o f t h e t o p s u r f a c e . 3.3 INTERNAL CONDUCTION PARAMETERS The t h e r m a l r e s i s t a n c e s b e t w e e n n o d e s a n d t h e t h e r m a l c a p a c i t y o f t h e t i m e i n t e r v a l from t h e m a t e r i a l p r o p e r t i e s o f t h e time. The m a t e r i a l p r o p e r t i e s ( d e n s i t y ( p ) , s p e c i f i c h e a t ( c ) , c o n d u c t i v i t y ( k ) , and e m i s s i v i t y ( E ) a r e i n p u ta st a b l el o o k up f u n c t i o n so ft e m p e r a t u r e .L i n e a ri n t e r p o l a t i o n is emn o d e sa r ec a l c u l a t e de a c h nodecorresponding t o i t s t e m p e r a t u r ea tt h a t ployedformaterialpropertydeterminationattemperaturesintermediateto t h o s et a b u l a t e d .C o n s t a n tt h e r m a lc o n t a c tr e s i s t a n c e s may b es p e c i f i e db e tweenany o ra l ln o d e s . The n o d a lc a p a c i t i e s and r e s i s t a n c e s are c a l c u l a t e d as follows : R m,n,A, e R m,n,B,e - 1 - 1 Lm,n,C Lm,n+l,A km,n+l,B km,n,B + + - 'm+l,n,D kA+l, n, F i g u r e 1 0 shows t hl eo c a t i o norsfe s i s t a n c e s R m,n,A,e + + R* m, n, B R* e m,n,A and R m,n,B,O' The I < Rm, n, B I m, n-1 I m,n m 0 Figure 10. , I Rm,n,A n + l m-l,n SketchofThermalResistanceNomenclature 15 two s i d e s of t h e n o d a l resistanceson the other t h es u r f a c e , however , box are c a l c u l a t e d when the are c a l c u l a t e d .F o rn o d e sa d j a c e n tt o q u a n t i t i e sf o rt h ea d j a c e n tn o d e s # Am ,n+l ,A g e n e r a l l y( r e f e rt oF i g . 6 f o ra r e a nomenclature);forthesenodes Rm,n,A Lm, n , c Lm,n+l, A + Am,n,C kmen, e + %,n+l,~ km,n+l, e Rm*,n,B Am,n, I t s h o u l db en o t e dt h a ta n i s o t r o p i ct h e r m a lc o n d u c t i v i t yv a l u e s IN-DEPTH k (17) and k' (15) and ( 1 6 ) , r e s p e c t i v e l y . a r eu s e di nE q u a t i o n s 3.4 c CONDUCTION SOLUTION The i n - d e p t hc o n d u c t i o ns o l u t i o n is theexplicitfinitedifferencetype o f t e n employed f o rt r a n s i e n th e a tc o n d u c t i o na n a l y s i s . The t e m p e r a t u r eo f node m,n a t t i m e 8' ( T m , n l e l ) i s o b t a i n e d b y a p p l i c a t i o n o f t h e f i n i t e d i f f e r e n c ee n e r g yb a l a n c e and r a t e e q u a t i o n s t o t h e n o d a l Solvingfor 1 Tmlnle,, 1 Rm,n-l,Al 8 Inthe + R o n eo b t a i n s 1 + 1 . m , n , ~R, em , n , ~R, em - l r n r ~ , e program t h i s e q u a t i o n a l l n o d e se x c e p tf o r and f o r backwallnodes. volume. 1 )]e+ T m, n, 8 is usedtoobtain mrnre "new" t e m p e r a t u r e s f o r two nodes i n each column a d j a c e n t t o t h e h e a t e d s u r f a c e The n o d e s n e a r t o t h e h e a t e d s u r f a c e a r e l i n k e d t o thesurfacetemperatureimplicitly andhence a specialprocedure is u s e d f o r t h et e m p e r a t u r eo ft h e s en o d e s ,a sd e s c r i b e di nt h en e x ts e c t i o n . Back w a l l nodes include i n s i d et h eb r a c k e t so fE q u a t i o n zero. 16 (18) a quantity (18) , andhave T m-1 ,n, 8 f o r m a l l ye q u a lt o The e x p l i c i t r e l a t i o n t h e t i m e s t e p s i z e A0 = (18) imposes the f a m i l i a r s t a b i l i t y r e s t r i c t i o n The MACABRE programautomaticallyemploys 8'8 -. on a conservative stability equation for interior nodes: AB = q Normally f o rs t a b i l i t y ,t h ei n p u tp a r a m e t e r n o d e sa r en o tc o n s i d e r e df o r below:back TI i s less t h a nu n i t y S . urface t i m e intervalcalculation,as wallnodesincludethe terms Am , n , D ( h / 2 f w i l l b ee x p l a i n e d ~ U E ~ ~ T : , ~ i n, t~h )e denominator. The a u t o m a t i c s t a b i l i t y c r i t e r i o n c a l c u l a t i o n may b e s u p p r e s s e d f o r a n y node i f t h e u s e r is surethattheallowed time s t e p f o r t h a t node w i l l n e v e r b et h e minimum one f o rt h es y s t e m .T h i ss a v e s some computer t i m e . Alternatively,thestabilitycriterioncalculation ~f it i s n o tu s e d ,t h e n A 8 m u s tb es p e c i f i e d . entirely. 3.5 TEMPERATURES SURFACE OF NODES can b es u p p r e s s e d AND SURFACE POINTS Temperaturesofsurfacepointsaredeterminedeither by assignment (Option 2 ) o r by s u r f a c e e n e r g y b a l a n c e s d e s c r i b e d i n t h e n e x t s e c t i o n ( O p t i o n s 1 and 3 ) . t h en t h The s u r f a c ee n e r g yb a l a n c ed e t e r m i n e st h e column TA new s u r f a c et e m p e r a t u r eo f w i t ha ni m p l i c i ti t e r a t i o nt e c h n i q u e .S t a b i l i t yc o n s i d e r - ations dictate thst the first andsecond node t e m p e r a t u r e s a l s o b e t r e a t e d i m p l i c i t l y , and t h a t a n yt r a n s v e r s eh e a tc o n d u c t i o nl i n k( a c r o s sc o l u m n s ) f o rs u r f a c et e m p e r a t u r e sm u s tb ei m p l i c i t .T h i sl a t t e rr e q u i r e m e n t n o tl i n k e dt r a n s v e r s e l y .F i g u r e c o n d u c t i o np a t h sf o r is a the s u r f a c et e m p e r a t u r ep o i n t s 11 shows t h e i m p l i c i t and e x p l i c i t h e a t complex one t o m e e t ascolumnsrecede:hence are two t y p i c a l s i t u a t i o n s . 17 Boundary Layer Edge SurfaceTemperatures F i r s t Nodes Second Nodes Implicit Links n-1 n n+1 _"""""- BoundaryLayerEdge Explicit Links Surface n-1 n n+1 F i g u r e 11. Ingeneral node S k e t c ho fI m p l i c i t and E x p l i c i tT e m p e r a t u r eL i n k s Two T y p i c a l S i t u a t i o n s inFiniteDifferenceSolution,for terms, t h e e q u a t i o n f o r t h e is new t e m p e r a t u r e o f a surface ( n o t a s u r f a c ep o i n t ) - Tm,n, 8 1 Tm,n+l,e Tm,n,e Rm,n,A, e + Tm,n-l,e Tm,nre R m, n-1, A, 8 (20) This equation links the new t e m p e r a t u r eo ft h en e x tn o d e 18 new f i r s t node temperature down T m - l , n , e , . Tm , n , 6 ' ana to the A n o t h e re q u a t i o na l s ol i n k s T and Tm-l,n,e,; t h i s i s t he n e r g y b a l a n co en t h see c o n d node down. m,n,8' Thishasthe same form as E q u a t i o n (18) w i t h T replaced m,n, e and Tm-l ,n ,e r i n t h e conductance terms : by. T m,n, 8' ana T m - l , n , e + -+ Tm-2,n,~- T m - l , n , ~ Tm,n,e' Tm-l,n,e' Rm-2 ,n , B ,0 Rm-l,n,B,B 1 'm-1 ,nt e E q u a t i o n s ( 2 0 ) and ( 2 1 ) between them h a v et h r e e new unknown temperat i m e s t e p : Ts , n ,e ' , T m , n , e, , and Tm-l e , (where m h e r e d e n o t e st h e row i n d e x o f t h e s u r f a c e m d a l b o x i n Column n ) . The two e q u a t i o n s may becombined t o e l i m i n a t e T m - l , n , e , , leaving a simplelinearlinkbetween thesurfacepointtemperature Ts and t h e s u r f a c e nodetemperature T m,n, 0 ' ' ,nr8' Formally w e have a r e l a t i o n b e t w e e n t h e two unknowns as turesateach (as d e s c r i b e d i n t h e n e x t The s u r f a c ee n e r g yb a l a n c eh a st h eg e n e r a lf o r m section) convectionandchemicalenergy + terms (T s,n,e') radiation to w a l l - radiationout (T s,n,e') -- . l c t i o n a lr e l a t i o n s h i p . w h e r e . t h ep a r e n t h e s e sd e n o t e on a u n i ta r e ab a s i s .R e l a t i o n (22) The e q u a t i o n i s w r i t t e n may b es u b s t i t u t e di n t oE q u a t i o n (23) togivethenon-linearsurfaceenergybalanceEquationfor convectionandchemicalenergy + radiation t o wall - T s,n,e': terms ( T s,n,e') radiationout (T s,n,e' ) 19 When T T h ab s een found from t h ie s q u a t i o nt,h e srnr8' may be found from Equation (22). m,n,e' t h e method f o r f i n d i n g T T is d e s c r i b e di nS e c t i o n 4 below. ( I nO p t i o n 2 , s , n , 0' i s known immediately by user assignment and Equations ( 2 0 ) and ( 2 1 ) s,n,e' then determine 20 new s u r f a c e node temperature F osru r f a ceen e r g b y a l a n coe p t i o n s , Tm , n , e r and T m - l , n , e l a t once.) SECTION 4 IV-V COUPLING TO SURFACE THERMOCHEMISTRY SOLUTION GENERAL ASPECTS 4.1 The d i r e c t c o m p u t a t i o n a l l i n k b e t w e e n t h e s u b s u r f a c e s o l u t i o n a n d t h e is e f f e c t e dt h r o u g ht h es u r f a c ee n e r g yb a l a n c ec a l c u l a - s u r f a c es t a t es o l u t i o n * t i o n .T h i sl i n k a g e solutioninSection was d i s c u s s e df r o mt h ep o i n to fv i e wo ft h ei n - d e p t h 3 a b o v ea n dt h ee n e r g ye q u a t i o nw r i t t e n down i n s c h e m a t i c terms as E q u a t i o n s (23) and (24) detail. The e n e r g yb a l a n c e Figure 1 2 . w h e r et h ei n d i c a t e dc o n t r o l f l u x e sl e a v i n gt h ec o n t r o l . W e now must d i s c u s st h i se q u a t i o ni n i s i l l u s t r a t e di nF i g u r e more 12, R e p r e s e n t a t i o no fS u r f a c eE n e r g y During Ablation Terms volume i s f i x e dt ot h er e c e d i n gs u r f a c e .E n e r g y volume i n c l u d e c o n d u c t i o n i n t o t h e m a t e r i a l , r a d i a - t i o n away from t h e s u r f a c e , e n e r g y i n any f l o w o fc o n d e n s e dp h a s em a t e r i a ls u c h a sm e c h a n i c a lr e m o v a l ,a n dg r o s sb l o w i n ga tt h es u r f a c e .E n e r g yi n p u t st ot h e c o n t r o l volume i n c l u d e r a d i a t i o n i n f r o m t h e b o u n d a r y l a y e r a n d e n t h a l p y f l u x due t o t h e c h a r f l o w r a t e . The f i n a l i n p u t i n t h e s k e t c h i s denoted -9,. I t i n c l u d e sa l ld i f f u s i v ee n e r g yf l u x e sf r o mt h eg a s - p h a s eb o u n d a r yl a y e r . were coupled t o an e x a c t boundary-layer is, of c o u r s e , a complex f u n c t i o no ft h eb o u n d a r y - l a y e rs t r u c t u r e .I f ,o nt h e o t h e rh a n d ,t h ei n - d e p t hr e s p o n s e i s b e i n gc o u p l e dt o a s i m p l i f i e db o u n d a r y l a y e r scheme suchas a c o n v e c t i v e f i l m c o e f f i c i e n tm o d e l , as w i l l bedonehere, t h e n t h e t e r m -qaassumes t h ef o r mo f a correlationequation. Ifthein-depthresponsecomputation solution, the term -9, wouldbeobtainedfromthatsolutionprocedureand The computation of t h e s u r f a c e e n e r g y b a l a n c e r e q u i r e s f r o m t h e i n depthsolution a relationbetween the s u r f a c e t e m p e r a t u r e e n e r g yc o n d u c t i o ni n t ot h em a t e r i a l ,q c o n d . and t h e rate o f As n o t e d i n S e c t i o n 3 above, t h i s relationderivesnacurallyfromthefinitedifferenceenergybalanceforthe *For Option 1 c a l c u l a t i o n s . The s i m p l e rs u r f a c eo p t i o n s 4 . 5 and 4.6 below. discussedinSections w i l l be 21 the node j u s t u n d e r t h e s u r f a c e a n d i s a simple linear equation q cond = a2Tw With t h i s i n f o r m a t i o n , t h e s u r f a c e e n e r g y b a l a n c e c o n s i d e s a t i o n s a l l o w d e t e r - rate minationofthethermochemicalerosion I t w i l l be useful to keep in mind t h a t , mc andsurfacetemperature TW ' from t h i s p o i n t o f v i e w , t h e p u r p o s e a t any i n s t a n t i s t o p r o v i d e i n f o r m a t i o n a b o u t may b e w r i t t e n as q c o n d ( T w ) . The s u r f a c ee n e r g yb a l a n c ee q u a t i o n o ft h ei n - d e p t hs o l u t i o n where The r e l a t i o n q cond = f (Tw) i s d e l i v e r e d by t h e i n - d e p t h s o l u t i o n . O t h e rd e p e n d e n c i e so fi n t e r e s ta r e out F o rt h eo t h e rt e r m s , Tw, -qa , in If the h,, -- 'rad (Tw) w r i t e i n general we may qrad, 'rad out 2' mrr a hi = f u n c t i o n so fb o u n d a r y l a y e r - e d g ee n t h a l p y , p r e s s u r e ,l o c a l boundary-layer solution, l a w sf o rc o n s e r v a t i o n ofchemicalelements, chemicalequilibria and/or kinetic relations, upstream e v e n t s , mc boundarylayertransportaspectsoftheproblemare modeled by a f i l mc o e f f i c i e n ts c h e m e ,t h e nE q u a t i o n( 2 5 )c a nb en o r m a l i z e d on t h e mass t r a n s f e r c o e f f i c i e n t i n t h e c u s t o m a r y 22 manner: + b2' and r e l a t i o n s h i p ( 2 9 ) becomes -9, %'p eue cM * in " p u c eeM hw , hk = f u n c t i o n s of bounda r y l a y e r e d g e enp r e s s ut rhea,l p y , l a w s forconser3.ation of chemical e l e m e n t s , chemical e q u i l i b r i u m and/or k i n e t i c relations, ee ' cM a B; The g e n e r a t i o n o f t h e s e r e l a t i o n s h i p s i s t h eg o a lo fP h a s e s u r f a c et h e r m o c h e m i s t r ys o l u t i o n ,t ob ed i s c u s s e db e l o w .F o rt h ep r e s e n t , IV, the it may an i n i t i a l g u e s so ft h ed i m e n s i o n l e s s mass removal r a t e BA i s o b t a i n e di n some manner. With t h i s BL, t h eq u a n t i t i e s -4, /peueCM, h w , x rirQhQ/peueCM, and Tw a r e b eo b s e r v e dt h a tt h ee n e r g yb a l a n c ec o u p l i n gp r o c e e d sa sf o l l o w s : o b t a i n e df r o mt h es u r f a c et h e r m o c h e m i s t r ys o l u t i o n . The q u a n t i t i e sh c and so o b t a i n e d . The s u r f a c ee n e r g y qrad out b a l a n c e i s t h e nc o m p u t e d ,t h e qcond as a f u n c t i o no f Tw havingbeenprovidedby t h ei n - d e p t hs o l u t i o n .I ng e n e r a l ,h o w e v e r ,t h e sum o ft h e terms w i l l n o te q u a l z e r o ,b u t some a r r o r . An i t e r a t i o n p r o c e d u r e is t h e nu s e dt o select success i v e l yb e t t e r estimates of B L w h i c hd r i v et h ee r r o rt oz e r o .E x p e r i e n c e shows t h a t Newton's,.rocedure, i n which t h e d e r i v a t i v e o f t h e e r r o r w i t h r e s p e c t t o E: i s u s e d t o compute t h en e x tg u e s sf o r g i v e s good r e s u l t s . a r et h e nf o r m u l a t e du s i n gt h e Tw € 3 ; 4.2 COMPUTATIONAL APPROACHES 4.2.1 Pre-celculated T a b l e It is evidentthateachiterationinthesearchfor a s u r f a c ee n e r g y b a l a n c e ,i fp e r f o r m e d as d e s c r i b e da b o v e ,w o u l dr e q u i r e a new s u r f a c e c h e m i s t r y solution,generallyinthenearneighborhoodof many s u c h . p r e v i o u s s o l u t i o n s . If this s t a t e s o l u t i o n is t o b e o b t a i n e d f r o m a l a r g es u b p r o g r a m( r a t h e rf r o m , , t h i sc o n s t a n tr e p e t i t i o no ft h e r m o c h e m i s t r yc a l c u - s a y ,s i m p l ec o r r e l a t i o n s ) l a t i o n s would cqnsume an t h a t a tabularapproachin hand f o rp r e a s s i g n e d BL - l a c c e p t a b l e amount ofcomputer time. which t h e s u r f a c e s t a t e s o l u t i o n s a r e T h i ss u g g e s t s donebefore- v a l u e sw o u l do f f e rs i g n i f i c a n tc o m p u t a t i o n a l economies. As anexampleof t h i sa p p r o a c h ,c o n s i d e ro n e body p o i n t . f o rs i m p l i c i t y ,t h a tc h e m i c a lk i n e t i c sa r en o ti m p o r t a n t :t h e nt h ei n d e p e n d e n t v a r i a b l e si nE q u a t i o n( 3 1 )r e d u c et op r e s s u r e enthalpy H variables. r' and B h . Assume, P , b o u n d a r yl a y e rr e c o v e r y F i g u r e1 3r e p r e s e n t st h es p a c eo ft h e s et h r e ei n d e p e n d e n t The c o u r s eo f a t y p i c a lt r a n s i e n ts o l u t i o nm i g h tb ea sr c p r e s e n t e d 23 . .... .. .. i nt h es k e t c h . As time p r o c e e d s ,s o l u t i o n sp r o g r e s s t F i g u r e 13. RepresentationofCourseofIndependent VariablesinSurfaceHistory, OneBody Point t h r o u g ht h es p a c eo ft h et a b l ei n d e p e n d e n tv a r i a b l e s BA, H r , and P . dotsinthepicturerepresentsolutionssatisfyingthesurfaceenergybalance a s time p r o g r e s s e s .S u r r o u n d i n gt h e s ep o i n t sa r e The a number o f o t h e r p o i n t s e x a m i n e dd u r i n gt h ei t e r a t i o n st os a t i s f yt h es u r f a c ee n e r g yb a l a n c e .F o r s o t o s p e a k ,w i t h i n any i t e r a t i o n , t h e s o l u t i o n p r o c e d u r e f i n d s i t s e l f , cubeformedby t h eb r a c k e t i n gt a b u l a rv a l u e so f d e p e n d e nqt u a n t i t i e s BA, Hr, and P . a The h /peueCM have been prer51 R calculatedforthesetabularpoints;relevantvaluesofthesequantities Tw, -9, /peueCM, hw and f o rt h ec u r r e n ti t e r a t i o nv a l u e s of B;, H r , and P canthenbeformed by cube and t h e s u r f a c e e n e r g y b a l a n c e e q u a t i o n c a l c u l a t e d .I ft h ee n e r g yb a l a n c e is notsatisfiedto some p r e s e l e c t e dd e g r e e o fa c c u r a c y , a new v a l u eo f B h c a nb es e l e c t e da n dt h ep r o c e s sr e p e a t e d . i n t e x p o l a t i o ni n s i d et h e The i n t e r p o l a t i o n f e a t u r e of t h e t a b u l a r a p p r o a c h d r a s t i c a l l y r e d u c e s same time a l l o w i n g The p r e c a l c u l a t e dt a b l ea p p r o a c hh a st h ef o l l o w i n g t h e number o f s u r f a c e s t a t e c a l c u l a t i o n s w h i l e a t t h e s u f f i c i e n ta c c u r a c y . * advantages : 1. I np a r a m e t r i cs t u d i e s ,t a b l e so n c eg e n e r a t e da r eu s e a b l e f o r many d i f f e r e n t p r o b l e m s , y i e l d i n g e v e n g r e a t e r * 24 economy. Forexample, a t y p i c a l 2 0 0 0 t i m e s t e pp r o b l e mw i t ht h eu s u a l 5 iterations which p e r time s t e p would r e q u i r e 1 0 , 0 0 0 s u r f a c e s t a t e c a l c u l a t i o n s , 7 0 9 4 s p e e dc l a s s . A single r e q u i r e 1 t o 3 s e c o n d se a c hf o rm a c h i n e si nt h e 300 s t a t e e n t h a l p y table, on t h eo t h e rh a n d ,m i g h ti n v o l v eo n l ya b o u t 30 improvement. s o l u t i o n s (30 B; v a l u e s times 1 0 p r e s s u r e s ) , a f a c t o ro f M u l t i p l ee n t h a l p yt a b l e sr e d u c et h i sa d v a n t a g e ,b u tu s u a l l yo n l y a few e n t h a l p i e sa r er e q u i r e d . 2. Formostproblems, it i s d i f f i c u l t t o s p e c i f y a p r i o r i an of i n d e p e n d e n t t a b u l a r v a l u e s EA, H,-, - and P. An e x a m i n a t i o n o f t h e p r e c a l c u l a t e d s u r f a c e tables b e f o r e e x e c u t i o no ft h ea c t u a la b l a t i o n - p r o b l e m - p l u s - i n - d e p t h - s o l u t i o n c a nr e v e a l if t h e r e a r e a n y " h o l e s " i n t h e t a b l e s i n areas whereenergy terms a r ev a r y i n gr a p i d l y .D e s i r a b l e table p o i n t s c a nb ea d d e db e f o r et h ei n - d e p t hr e s p o n s er u n . "adequate" array 3. The s u r f a c e t a b l e s a r e f r e q u e n t l y o f i n d e p e n d e n t i n t e r e s t inthemselvesforjudgingtheablationeffectsundervarious conditions. 4. Finally,withoutthetabularapproach,theoccasional vergentsurfacechemistrysolutionwouldstoptheentirein- noncon- d e p t hs o l u t i o np r o c e s s . With t h ep r e c a l c u l a t e dt a b l ea p p r o a c h , s u c hs o l u t i o n sa r ea u t o m a t i c a l l y weeded o u t o f t h e t a b l e s w i t h o u t damage t o t h e s u b s e q u e n t i n - d e p t h s o l u t i o n s . On t h e o t h e r s i d e o f t h e l e d g e r , some b i g d i s a d v a n t a g e s i n t h e p r e c a l c u l a t e dt a b l ea p p r o a c ha r ee v i d e n th o w e v e r : 1. F i g u r e 1 3 , which i s a r e a l i s t i c s c h e m a t i c view Of a typical calculationhistory,indicatesthat m o s to ft h el a b o r i o u s l y c a l c u l a t e d a n da s s e m b l e ds u r f a c es t a t ep o i n t si nt h et a b l e a r en e v e ru s e di nt h ec o u r s e 2. of a g i v e n s o l u t i o n . state solution The " m e c h a n i c a l "l i n k a g eb e t w e e nt h es u r f a c e and t h e i n - d e p t h s o l u t i o n , i.e., thetransferof punchedcard surfacestateoutputtoinput of t h ei n - d e p t hp r o g r a m ,l e a d s t o computingdelaysandoccasionally (wrongdecks,missingpartsofdecks togrossinputblunders , etc.) . 3. Too m-xh s t o r a g e i s r e q u i r e d by t h e l a r g e p r e c a l c u l a t e d table. 4. The s y s t e m c a n n o t r e a d i l y b e g e n e r a l i z e d t o a c c o m o d a t e more p a r a m e t e r s .B o t hs t o r a g er e q u i r e m e n t sa n dc o m p u t a t i o n time grow o u t o f r e a s o n a b l e b o u n d s . T h i s l a s t d i f f i c u l t y wouldprove to be a great stumbling block in t h ea b l a t i o nc o d ed e s c r i b e dh e r e ,w h i c hc o n s i d e r sc h e m i c a lk i n e t i c si n additiontotheotherindependentvariables BA, H,, and P d e s c r i b e d a b o v e . " K i n e t i c s "c o n t r i b u t e a f o u r t hp a r a m e t e ri nt h ef o l l o w i n g manner: The a number o f s i m u l t a n e o u s e q u a t i o n s i n c l u d i n g and k i n e t i c a l l y c o n t r o l l e d r e a c t i o n s Each k i n e t i c a l l yc o n t r o l l e dr e a c t i o ne q u a t i o nh a so n e chemicalstateroutineconsiders elementbalances,equilibriumrelations, ( R e f e r e n c e s 4 and 5 ) . forward r a t e c o e f f i c i e n t , which i n t u r n c o m p r i s e s o n e p r e - e x p o n e n t i a l f a c t o r B m ( f o r t h e masuch r e a c t i o n ) and a t e m p e r a t u r ed e p e n d e n tq u a n t i t y . All 25 I e q u a t i o n sa r en o r m a l i z e d on t h e t r a n s f e r c o e f f i c i e n t desiredindependentvariable B' peueCg t h i s y i e l d s intheelementbalances n o r m a l i z e dk i n e t i cr a t ec o e f f i c i e n t s Bm/peueCM. the and a l s o p r o d u c e s Thus a n yf o r w a r dr a t ec o e f - ficientsusedinthestateroutineareinputinthe form Bm/peueCM, and t h e whole s e t o f Bm/peueCM v a l u e s c h a r a c t e r i z e s t h e r o l e o f k i n e t i c s . 1 n a given s e t of B,'s i s onlyoneparameterzaling the k i n e t i c e f f e c t s . E x ~ a l t l y what t h i sp a r a m e t e r i s c a l l e d i s n o ti m p o r t a n t ; l e t i t betermed B1/peueCM where B1 r e p r e s e n t s t h e f i r s t o f t h e s e t o f k i n - problemthe is f i x e d ; t h u s r e a l l y t h e r e e t i c a l l yc o n t r o l l e dr e a c t i o n sc o n s i d e r e d . Thus t h e f o u r p r o b l e m v a r i a b l e s , w h i c h w o u l d b e t h e i n d e p e n d e n t table, are: v a r i a b l e si nt h ep r e p a r e d 1. A b l a t i o nr a t e 2. K i n e t i c sp a r a m e t e r 3. Recovery enthalpy 4. Pressure BA B1/peUeCM F o u ri n d e p e n d e n tv a r i a b l e sa r ep l a i n l yv e r yt r o u b l e s o m e . t h es a v i n g si nc o m p u t a t i o n time a r e s e r i o u s l y First, compromised, s i n c e t h e number o fp r e p r e p a r e ds o l u t i o n sh a sg o n e up by a l a r g e nvrmber. Second,machine are now i n c o n v e n i e n t ;s e v e r a ld e p e n d e n tv a r i a b l e sm u s t b es t o r e da te a c hp o i n ti na na r r a yw h i c hm i g h ti n c l u d e 30 x 1 0 x 1 0 x 1 0 = 3 0 0 0 0 p o i n t s .T h e s ed i f f i c u l t i e si nt h ep r e - c a l c u l a t e dt a b l ea p p r o a c h lead to the invention of a revisedtabularapproachdescribedinthenext sub-section. s t o r a g er e q u i r e m e n t s 4.2.2 4.2.2.1 Direct-Coupled Table G e n e r aD l escsiption I t is obviousthatthedifficultiesinvolvedintheprecalculated tableapproachto a " f o u r - d i m e n s i o n a lt a b l ep r o b l e m "c o u l db e removed w i t h a " c l o s e - c o u p l i n g ' 'p r o c e d u r ew h i c ho n l yc a l l e df o rc a l c u l a t i o no ft a b u l a r * v a l u e sa st h e ya r en e e d e dd u r i n gt h ea c t u a li n - d e p t hs o l u t i o n . A t any i n s t a n t ,t h es u r f a c e - e n e r g y - b a l a n c e - p l u s - i n - d e p t hr e s p o n s es o l u t i o np r o c e d u r e f i r s t examines t h ec o r n e r so ft h e" s q u a r e "( o r" c u b e " ,d e p e n d i n g number o fi n d e p e n d e n tv a r i a b l e s ) it f i n d s i t s e l f i n , t o on t h e see i f t h e n e c e s s a r y d e p e n d e n tq u a n t i t i e sh a v eb e e nc o m p u t e da n ds t o r e d .I fn o t ,t h ei n - d e p t h routinereturnsto v a l u e s be s u p p l i e d . a m a s t e rc o n t r o lr o u t i n e and a s k s t h a t any m i s s i n g The m a s t e rc o n t r o lr o u t i n ed e t e r m i n e sw h i c hc o r n e r s *Another way o u t , n o t c o n s i d e r e d h e r e , is t o r e d u c e t h e number o f p a r a m e t e r s fromfourto twoby i g n o r i n g k i n e t i c s andbyassumingequaldiffusion coF o rt h i s" c l a s s i c a l "s i m p l ec a s e ,a s is w e l l known, e f f i c i e n t s and CM = CH. theedgestatedoesnotenterthedeterminationofthesurfacestate, and B ' and P a r e t h e onecanusethe sametwo parametersurfacetableinwhich independentquantities,forallsurfacelocations and a l l tfmes. 26 I. " . . I. ."..... ......""""._ needfilling, calls t h e s u r f a c e s t a t e program t o f i l l them,andthenreturns tothein-depthroutine,whichcontinueswiththeproblemsolution. Suchanapproachminimizes the number o f s t a t e p o i n t s c a l c u l a t e d , since the four parameters in the list above shrink to only three because a n de d g es t a t e ) can now time = t h e s o l u t i o n p r o g r e s s e s , a l t h o u g h when t h es e c o n d ,t h i r d ,a n df o u r t hv a r i a b l e s be directly associated with (B1/peueCM t h i s i s done a x i a l s t a t i o n s m u s t b e d i s t i n g u i s h e d f r o m e a c h o t h e r , s i n c e each station w i l l have different t i m e h i s t o r i e s . T h e r e i s a n e t d e c r e a s e of one i n t h e number o f i n d e p e n d e n t v a r i a b l e s . I f t h e t a l e s were t o b e p r e c a l c u l a t e d ,t h et h r e ev a r i a b l e sw o u l d ,o fc o u r s e ,r e q u i r es e p a r a t ea r r a y s s p e c i f i e di na d v a n c e ,s i n c et h et r u eh i s t o r i e s c o n s e q u e n c e so ft h es o l u t i o nd u et o of p u C andedge s t a t e a r e eeM body s h a p e e f f e c t s . With d i r e c t c o u p l i n g , E / p u C and t h ee d g e s t a t e a r e e f f e c t i v e l y known f u n c t i o n so f 1 eeM t i m e foreachstationsincethebodyshape is to be calculated a t each instant. Thus t h ei n d e p e n d e n tv a r i a b l e 1. A b l a t i orna t e 2. A x i aslt a t i o n l i s t becomes: E' T i3m . e Thus t h ed i r e c tc o u p l e da p p r o a c hb r i n g s s y s t e ma n dt h e u s back t o a t h r e e - p a r a m e t e r same generaleconomicsasbefore.Anothercomputational economy i s i n t r o d u c e d s i n c e i n g e n e r a l t h i s a r r a y filledinthecourseof w i l l notbecompletely a t r a n s i e n tp r o b l e m . A d d i t i o n a le c o n o m i e so fc o m p u t e rs t o r a g ea r er e a l i z e ds i n c eo n l y two t a b l e v a l u e s o f t i m e areinvolvedinthecalculationat I n s t e a d of h a v i n g t o s t o r e F i g u r e 13, t h e programneeds o fo n eb o xa si n d i c a t e di nF i g u r e a wholeboxofsolutionsofthesort to store at 14. any one i n s t a n t . shown i n any one t i m e o n l y a d j a c e n t " s h e e t s " A s s o o na st h e time p a s s e so u to ft h e i n t e r v a lb e t w e e nt h e two s h e e t s , t h e e a r l i e r ( O n ) o ft h e s e t w o sheets is d i s c a r d eadn d a new s h e e t a t ( t e m p o r a r i l y empty o f any dependent s o l u t i o nv a l u e s )a d d e d . Thus t h e d i r e c t c o u p l e d a p p r o a c h p u t s t h e "more-than-three-parametert a b l e " p r o b l e mw i t h i nr e a c he c o n o m i c a l l ys p e a k i n g ,a n dr e n d e r s it more 27 4 % 8 dimensionincludes effects of t r a n s i e n t , shapedependentedge s t a t e and p,ueCM 2 4 % +, v) rl d Figure 1 4 . S k e t c hI n d i c a t i n gS o l u t i o nS p a c ea n d S o l u t i o n" S h e e t s "S t o r e da t Any I n s t a n t During Solution e c o n o m i c a lt h a nt h ep r e c a l c u l a t e dt h r e e - p a r a m e t e rt a b l ea p p r o a c h .F u r t h e r more, it r e d u c e sc o m p u t e rs t o r a g ec o n s i d e r a b l ya n dr e d u c e sm e c h a n i c a lh a n d l i n g problemsofpunched tables. However, i n t h i s d i r e c t - c o u p l e d schemetwo i m p o r t a n tp r o c e d u r a l w i l l b ed e s c r i b e di nt h ef o l l o w i n g m a t t e r sh a v et ob ec o n s i d e r e d .T h e s e sub-section. Some Dangers i nt h eD i r e c t - C o u p l e dT a b l eA p p r o a c h 4.2.2.2 is The g e n e r a l a t t r a c t i v e n e s s o f t h e d i r e c t - c o u p l e d t a b l e a p p r o a c h somewhatcompromisedbytwo 1. potentialproblems: The u s e r w i l l probablywant t a b u l a ri n d e p e n d e n t Bh t os p e c i f yi na d v a n c et h ea r r a yo f values.Butwhat if he makes a p o o r c h o i c e , s o t h a tt h ei n - d e p t hs o l u t i o nh a st r o u b l ef i n d i n g a s u r f a c e - e n e r g yb a l a n c ew i t ht h et a b u l a rv a l u e sp r o v i d e d ? Can the p r o b l e m s o l u t i o n b y a u t o - a scheme b e d e v i s e d t o r e s c u e m a t i c a l l yi n s e r t i n q new t a b u l a r E; v a l u e si nu n f o r e s e e nc r i t i c a l ordifficultareasofthetable? 2. What c a nb ed o n ea b o u tt h eo c c a s i o n a ln o n c o n v e r g e n ts u r f a c e s t a t es o l u t i o no b t a i n e d from t h e s u r f a c e s t a t e shows t h a t i t i s n o t p o s s i b l e t o g u a r a n t e e fromcomplex a r r a yu s u a l l y s u r f a c ec h e m i s t r yr o u t i n e s :i n t h es u r f a c es t a t e 28 any s u r f a c e t a b l e some of t h e p o i n t s p r o v e t o b e u n o b t a i n a b l e w i t h o u t some hum.an i n t e r v e n t i o n i n t h e schememust program?Experience a c o n v e r g e ds o l u t i o n program. form o f b e t t e r f i r s t q u e s s e s f o r I nt h ec o u p l e dp r o g r a m ,t h ei n - d e p t h beabletoadjusttothisfact a n dc o n t i n u ec a l c u l a t i n g I of a c o n v e r g e d e n t r y i n t h e s u r f a c e even i n t h e o c c a s i o n a l a b s e n c e s t a t e table. an automaticanswer The c o d e d e s c r i b e d h e r e d o e s n o t c o n t a i n to the f i r s t problem. A somewhat s i m i l a rc o u p l e dc o d e( R e f e r e n c e 6) u s e df o r charringmaterialsofvery complex compositionhas shown t h a t s u c h a f e a t u r e i s less necessarythanmightbesupposed a priori.In w o r t h w h i l es c h e m ew o u l db ev e r yd i f f i c u l tt od e v i s e . any e v e n t , a g e n e r a l l y The codedoescontain, however, a v e r y v a l u a b l e f e a t u r e w h i c h s u c c e s s f u l l y meets most of t h e d i f f i c u l t i e s which a r i s eu n d e rp o i n t number 1. above.Carbon a b l a t i o ni na i r p r o v i d e s a worthwhileexample.Figure o fc a r b o na b l a t i o nr a t e si na i r . 1 5 shows t h e dependence on t e m p e r a t u r e Over a widerangeoftemperature B& i s f o r a l lp r a c t i c a lp u r p o s e si n d e p e n d e n to ft e m p e r a t u r e ;i nt h i sr a n g et h eo n l y u s e f u li n d e p e n d e n tv a r i a b l ei nt h es u r f a c et h e r m o c h e m i c a ls o l u t i o np r o c e s s i s t e m p e r a t u r e ,s i n c e i t wouldbenearlyimpossible t o s e l e c t a sequenceof Bd v a l u e s d i s t r i b u t e d a l o n g t h e p l a t e a u . Therefore , the description presented above, in which BG was usedas o n eo ft h ei n d e p e n d e n tv a r i a b l e s ,s h o u l db ee x t e n d e dt oi n c l u d et e m p e r a t u r e a s an i n d e p e n d e n tv a r i a b l ei np l a c eo f BG, a t l e a s t i n some c a s e s .F i g u r e a l s o shows,however, t h a tt e m p e r a t u r e 15 i s n o t a good i n d e p e n d e n t v a r i a b l e i n some i n s t a n c e s . I I I + K i n e t i cPsl1a t e a u Regime Regime F i g u r e 15. A t h i g ht e m p e r a t u r e s T I I I Sublima.tion Regime A b l a t i o nR a t e B& VersusTemperature f o r Carbon i n A i r B & becomes stronglydependenton T , s o t h a t BA i s t h e preferredindependentvariable. T h i sg e n e r a lp r o b l e m i s r a t h e rt y p i c a l , i t t u r n so u t .C o n s e q u e n t l y thesingleindependentvariablerepresented by B b i n F i g u r e 14 h a s , i n t h e c o d e , a more complexdualnature. The v a r i a b l eb e g i n sw i t h a s t r i n go f a s c e n d i n gs u r f a c et e m p e r a t u r e( n o t B&) values(chosen by t h eu s e r )a n d concludeswith a string of ascending B b values(alsochosen by t h e u s e r ) . 29 MACABRE code calls f o r s u r f a c e s t a t e s o l u t i o n s During the s o l u t i o n p r o c e s s , t h e a t t h e minimum ( f i r s t ) B; Bi s t r i n g a n d n o t e s t h e t e m p e r a t u r e valueineach d i s c o v e r e df o rt h i sp o i n t .T e m p e r a t u r ep o i n t si nt h et e m p e r a t u r es t r i n ga b o v e t h i st e m p e r a t u r ea r ei g n o r e d .I ns o l v i n gf o rt h es u r f a c ee n e r g yb a l a n c ea t a givenpointat a g i v e n t i m e , t h e codecompares f a c et e m p e r a t u r e T o ft h e W t ot h e thepreviously formed s u r f i r s t B l v a l u e sf o re a c h two t e m p e r a t u r e s a t t h e two t i m e s h e e t s c u r r e n t l y s t o r e d . I f TW i s less t h a n e i t h e r o f t h e s e w i l l refer only to two v a l u e s , t h e s u r f a c e e n e r g y b a l a n c e i t e r a t i o n p r o c e s s t h eu s e r - a s s i g n e dt e m p e r a t u r ep a r to ft h et a b l e( w i t h BGls and o t h e rd e p e n d e n t informationateachtabularpointdiscoveredasnecessary by c a l l s t o t h e w i l l refer t ot h e c h e m i s t r yr o u t i n e s ) .O t h e r w i s e ,t h ei t e r a t i o np r o c e s s user-assigned B ' p a r to ft h et a b l e( w i t h T w ' s and o t h e rd e p e n d e n ti n f o r m a t i o n d i s c o v e r e d i n a s i m i l a r manner a s n e e d e d ) . Thisdualnature mechanism s o l v e s many o f t h e d i f f i c u l t i e s o t h e r w i s ea r i s eu n d e rp o i n t 1. above. caution,however.Forexample, The u s e rm u s t w h i c hn i q h t s t i l l e x e r c i s e some i nt h ec a s ei l l u s t r a t e d the by F i g u r e1 5 , code u s e r m i g h t b e g i n h l s s e q u e n c e o f i n d e p e n d e n t v a r i a b l e s w i t h t e m p e r a t u r e e n t r i e ss t a r t i n ga tt h el o w e s tt e m p e r a t u r eo fi n t e r e s t( s a y , e x t e n d i n g up t o ( a s a minimum o b j e c t i v e ) t h e p o i n t w h e r e 53O0R) and Bd b e g i n s t o a b o v et h ep l a t e a uv a l u ef o rt h el o w e s tp r e s s u r et ob ee n c o u n t e r e d . rise I t would "rise b ec o n s e r v a t i v ea n dp r e f e r a b l et oa d dt e m p e r a t u r e se v e nb e y o n dt h e pomt"atthehighestpressureincasethesepointsshouldbeneededduring thesolution. Afterthesetemperaturevalues, rangeof t h ep l a t e a uv a l u e . plateauinthe 1. Bdls may beadded The u s e r m u s t b e c a r e f u l t o p i c k Bd 1s. to span the expected a minimum Bd s l i g h t l y above A Bd v a l u e n e a r t h e s t a r t o f t h e p l a t e a u o r o f f t h e low t e m p e r a t u r e r e g i o n The l o w e s ta s s i g n e d is e x t r e m e l y u n d e s i r a b l e f o r two r e a s o n s : BG d e t e r m i n e st h eb r e a kb e t w e e na s s i g n e d Bg s o l u t i o n s ; a low Bd w i l l t e m p e r a t u r es o l u t i o n sa n da s s i g n e d ineffect"disqualify"alloftheassignedtemper--.%urepoints a c r o s st h ep l a t e a u and l e a v e a l a r g e " h o l e " i n t h e a r r a y o f s u r f a c et h e r m o c h e m i c a ls o l u t i o n s . 2. F o re q u i l i b r i u mc a l c u l a t i o n s ,a s s i g n e d below t h e p l a t e a u h a v e Bd s o l u t i o n s a t Bd v a l u e s a highprobabilityofnon-convergence difficultiesinthechemistryroutines. The s e c o n d p r a c t i c a l d i f f i c u l t y c i t e d a b o v e f o r t h e d i r e c t c o u p l e d t a b l ea p p r o a c hc o n c e r n sp o s s i b l en o n - c o n v e r g e n ts o l u t i o n si nt h ec h e m i s t r y r o u t i n e s .S t r i c t l ys p e a k i n g ,e v e no n es u c hs o l u t i o ns h o u l ds e r v et oh a l t theentiresolutionprocess,butfortunatelythesituationusuallydoesnot demand s u c hd r a s t i ca c t i o n . of iterations to find the 30 The c h e m i s t r yr o u t i n e sa l l o w a c e r t a i n number state solution;failuretoconvergewithinthepre- s e t a l l o w e d number of i t e r a t i o n s i s r a r e a n d u s u a l l y means t h a t t h e s o l u t i o n e i t h e r is oscillating in the very near neighborhood of the true solution or has arrived near to the solution and i s p r o c e e d i n g v e r y s l o w l y t o w a r d it. Therefore,usuallythe MACABRe code may s a f e l y a c c e p t a c e r t a i n number o f ' non-converged answers. Beyond t h i s l i m i t , t h ec o m p u t a t i o n i s s t o p p e d as a are u s u a l l y d u e t o l a r g e safetymeasure,sincerepeatedfailurestoconverge input blunders. 4.2.2.3 The O t h e r I n d e p e n d e n t V a r i a b l e s i n t h e S u b - s e c t i o n 4.2.2.2 o ft h e TJBA abovehasdescribedin Direct Coupled T a b l e some d e t a i l t h e u s e r i n p u t s t r i n go fi n d e p e n d e n tv a r i a b l ev a l u e s .T h e r ea r e ,o fc o u r s e , The f i r s t o f t h e s e two are q u i t e s i m p l e . additional independent variables, but fortunately these i s simply body s t a t i o no rs u r f a c ep o i n t .T h e s e i t s own s t r i n g o f TJBA v a l u e s . Thus wouldappear i n F i g u r e 1 4 a s a i t s own h o r i z o n t a l number o f s t r i n g s o f s o l u t i o n p o i n t s , e a c h w a n d e r i n g i n p l a n e ( t h e t i m e -BA p l a n e ) . a r ek e p ts e p a r a t e ,e a c hp o i n th a v i n g thehistoriesofallthesurfacepoints t i m e dimension, is represented t i m e v a l u e s . The same a r r a y w i l l b e u s e d f o r a l l stations. I t is n o t d i f f i c u l t f o r t h e u s e r t o s e l e c t a s e t ofvaluesappropriate to the transient nature of a givenproblem. The l a s t dimension t o b e c o n s i d e r e d , t h e by a u s e r s e l e c t e d s t r i n g o f any g i v e n time, t h e Thus,inevaluatingthesurfaceenergybalanceat surface energy balance routine refers to s t o r e d t a b u l a r v a l u e s o f T w and enthalpyquantitiesatthe two t a b u l a r v a l u e s o f t i m e s u r r o u n d i n gt h ec u r r e n t time. I ft h er e q u i s i t ei n f o r m a t i o nh a sn o ty e tb e e nc a l c u l a t e d and s t o r e d , theFILLroutinecallsthethermochemicalstateprogramtosupplythisinformation. 4.2.2.1 I f a " p r e v i o u s time" p o i n t( a t andFigure naturallyevaluated 8 i nt h en o m e n c l a t u r eo fS e c t i o n is of t h e l o c a l i n d e p e n d e n t v a r i a b l e s i n t h e 1 4 ) i s beingfilledin,thethermochemicalstate a t storedvalues * c h e m i s t r ys o l u t i o n ,p r e s s u r ea n dr e c o v e r ye n t h a l p y ,a tt h a tp r e v i o u s I f a " f u t u r e time" t a b l e p o i n t is b e i n g f i l l e d i n , t h e F I L L r o u t i n e to obtain the state solution at future local values of pressure e n t h a l p y .I ft h eu s e r i s n o ts u p p l y i n gt h e s ev a l u e sa sf u n c t i o n so f time. is obliged andrecovery time ( a n dt h e r e b ys u p p r e s s i n g some of t h ea u t o m a t i cc o u p l i n gf e a t u r e so ft h ec o d e ; see S e c t i o n 7 b e l o w ) ,t h e nt h e s ef u t u r ev a l u e sm u s tb ee s t i m a t e d .F u t u r e l o c a l p r e s s u r e is e s t i m a t e d as t h e f u t u r e s t a g n a t i o n p r e s s u r e times t h e p a s t ratio of local pressure to stagnation pressure (i.e. , basedonthecorrect . f u t u r es t a g n a t i o np r e s s u r e ,b u tt h ep a s t body s h a p e ) S i m i l a r l y ,f u t u r e l o c a lr e c o v e r ye n t h a l p y is e s t i m a t e da st h ef u t u r es t a g n a t i o ne n t h a l p y times the ratio of the past local recovery enthalpy to the past stagnation enthalpy. * The t r a n s f e r c o e f f i c i e n t p ueCM i s a l s o a n i n d e p e n d e n t v a r i a b l e , a s d i s c u s s e d above; i f k i n e t i c s a r e b e i g g c o n s i d e r e d it i s h a n d l e d d i f f e r e n t l y , a s will be explained below. 31 refers o n l y t o " f u t u r e t i m e " communicationwiththe state r o u t i n e , f o r which an a p p r o x i m a t e f u t u r e v a l u e of H r i s q u i t e a d e q u a t e s i n c e H r h a s o n l y a s l i g h t e f f e c t on t h e s u r f a c e t h e r m o c h e m i c a l state evaluation: inthesurfaceenergy,balanceequationevaluation i t s e l f , the"properly" computed c u r r e n t r e c o v e r y e n t h a l p y i s used.) (Thisapproximation 4.2.3 Summary a surface state Threeprocedureshavebeenconsideredforcoupling s o l u t i o nt ot h et w o - d i m e n s i o n a li n - d e p t hs o l u t i o n :d i r e c tc a l c u l a t i o na t tables, and d i r e c t c o u p l e d t a b l e s f i l l e d i n e v e r yi t e r a t i o n ,p r e c a l c u l a t e d The t h i r d w a s shown t o b e t h e o n l y r e a s o n a b l e a st h es o l u t i o np r o g r e s s e s . p r o c e d u r ef o r a t w o - d i m e n s i o n a li n - d e p t hs o l u t i o n . A s u i t a b l ec o u p l i n g a logicroutinetodeterminewhich table p o i n t s n e e d fillingduringthesurfaceenergysolutionprocess,plus some means o f procedurerequires table i n d e p e n d e n tv a r i a b l e s .I nt h ec o d e ,t h e s p e c i f y i n gt h ev a l u e so ft h e userspecifies f o re a c h a priorithesequenceof body s t a t i o n , p l u s v a l u e sb ec o n s i d e r e da p p r o p r i a t e B' a series o f time p o i n t s a t w h i c hs o l u t i o n s are to be performed. 4.3 THE CONVECTIVE SURFACE ENERGY BALANCE EQUATION w e havebeenabletodiscussthe Up t o t h i s p o i n t , thein-depthsolution and t h e s u r f a c e t h e r m o c h e m i c a l IV-V l i n k b e t w e e n s t a t e s o l u t i o n , which i s e f f e c t e dt h r o u g ht h ec o n v e c t i v es u r f a c ee n e r g yb a l a n c eo p e r a t i o n , w i t h o u t discussingthesurfaceenergyequation , g i v e n by Equations(23) (24) t oc o n s i d e rt h i se q u a t i o ni n , itself except in (25) , ( 3 0 ) , and ( 3 1 ) . more d e t a i l . its generalnature, I t i s now a p p r o p r i a t e As n o t e dp r e v i o u s l y ,t h es u r f a c e l i n k a g e s are b u i l t upon a t r a n s f e r c o e f f i c i e n t a p p r o a c h : t r a n s f e r c o e f f i c i e n t equationsareusedfor mass d i f f u s i o n i n t h e e l e m e n t b a l a n c e e q u a t i o n s and i s employed a transfer coefficient type of surface energy balance equation i n t h e IV-V l i n k a g e . The p r o p e r c h o i c e o f e n e r g y e q u a t i o n c o n s t i t u t e s a subject f a r t o o l a r g e and d i f f i c u l t t o d i s c u s s c o m p l e t e l y i n t h i s r e p o r t . The MACABRE code u s e so n eo ft h r e ed i f f e r e n tf o r m so ft h ec o n v e c t i v ee n e r g ye q u a t i o n , d e p e n d i n go nu s e rc h o i c e .I na s c e n d i n go r d e ro fc o m p l e x i t y ,a n dd e s c e n d i n g o r d e ro f " w e l l f o u n d e d n e s s ", t h e s ee q u a t i o n sa r e (1) F o re q u a l mass d i f f u s i o nc o e f f i c i e n t sf o rb o u n d a r yl a y e r moleculesand Hr - f o r C M = C H ( L e = 1) : ( l + B ' ) h w+ B ' h Fo_ E -__ PeuecH 32 c c Tw4 = + OLWqradin 'eUecH qcond 'eUecH T h i se q u a t i o n i s a " s t a n d a r d " (see, f o r example,Reference 7) one and r e q u i r e s n o d i s c u s s i o n . (2) CM # CH ( L e # 1): F o re q u a ld i f f u s i o nb u t + awqrad i n - 4 FU€TW ~ - 'cue% 'eUecH SO # C forchemicallyactiveboundarylayershasbeen remainsopen(References (3) 7, 8 , 9) CH e f f e c t s , a n dd o u b t l e s s M The accuracy i s frozen. if t h eb o u n d a r yl a y e r (33) 'eUecH T h i se q u a t i o ns u p p o s e d l ya c c o u n t sf o r does = qcond Of Equation ( 3 3 ) much d i s c u s s e d a n d t h e q u e s t i o n . F our n e q u adl i f f u s i o n : 4 FU ETw + - - -OLwqrad -"- 'eUecH "euecH -~ qcond 'eu,% Thisequation,introducedinReference = (34) 10 , a p p e a r s t o b e u s e f u l f o r m o d e l i n gu n e q u a ld . f f u s i o n effects, atleastinfrozenboundarylayers,and currently is b e i n g s t u d i e d i n d e t a i l t o e s t a b l i s h i t s validity(Reference 9). All three o ft h e s ee q u a t i o n s ( 3 2 ) , ( 3 3 ) , and ( 3 4 ) s e r v e t o p r o v i d e i nE q u a t i o n ( 3 0 ) . The conexpressionsforthediffusionenergyflux -q a , whichincludeh. s t i t u e n t s o f -9, ZZ; h y , and 2 Z* h? e iw Wedge gas are g e n e r a t e d by t h e c h e m i s t r y r o u t i n e s a l o n g w i t h hw and Tw. The s u b r o u t i n e FILL,which c a l l st h ec h e m i s t r yr o u t i n e s , thencombines the q u a n t i t i e s i n the a p p r o p r i a t ef a s h i o na sd i c t a t e d by E q u a t i o n ( 3 4 ) o r o n e o f t h e s i m p l e r forms ( 3 3 ) and ( - 2 ) . Forexample,suppose the FILL r o u t i n ed i s c o v e r s a r e q u i r e m e n t f o r a s t a t e s o l u t i o n a t a t a b u l a r B ' and f o r a g i v e n p r e s s u r e C a n dr e c o v e r ye n t h a l p y . FILL c a l l s t h e s t a t e r o u t i n e s , c o n t r o l l e d by E Q U I L , whichperform.the s t a t e s o l u t i o n ,e v a l u a t ek e yp r o p e r t yv a l u e s , and r e t u r n t o FILL v a l u e s Tw , hw , P Z* hiTw , hw, Zz h? FILL t h e n ,r e f e r r i n g e le W . 33 tothetemperature structsthe Tw, looks up t h e a b l a t i n g m a t e r i a l e n t h a l p y h c and con- term of i n t e r e s t . The t e m p e r a t u r e (On t h ea s s i g n e d EQUIL r e t u r n s BA and t h i s i s s t o r e d . ) and s t o r e s i t as o n eo ft h ed e p e n d e n tq u a n t i t i e s Tw i s s t o r e da st h es e c o n dd e p e n d e n tq u a n t i t yo fi n t e r e s t . Tw l o w e r p a r t o f t h e t a b l e , 4.4 CONVECTIVE SURFACE ENERGY BALANCE I T E R A T I O N PROCEDURE A t everysurfacepoint a t which a c o n v e c t i v e e n e r g y b a l a n c e i s t o be o b t a i n e d ,t h ee n e r g yb a l a n c es u b r o u t i n em u s t decidewhether Tw o r B ' i s theproperindependentvariable to iterate s e l e c t a trial v a l u e o f t h i s i n d e p e n d e n t v a r i a b l e refertostoreddependentvaluesofenergyquantity Tw ( o r B; andtemperature Y if in first part of table), inter- p o l a t i n ga sn e c e s s a r y on BA and 8 , and f i l l i n g i n needed FILII t a b l ep o i n t sa sn e c e s s a r yb yc a l l i n q look up aW(TW)and C H , and H r ( e ) EW (Tw) , grad (8), l o c a l F , p l u s CM, in (which i n many cases a r e g e n e r a t e d s u b r o u t i n e sd e s c r i b e di nS e c t i o n s 8 and 9 below) f o r m u l a t et h el e f t - h a n ds i d eo fE q u a t i o n f r o mz e r o( e r r o r (34), n o t e d e p a r t u r e E) if e r r o r i s t o o l a r g e , c o r r e c t , go t o ( c ) variable by s p e c i a l - G e n e r a l l yt h i sp r o c e s sc o n v e r g e sv e r y t r i a l v a l u eo fi n d e p e n d e n t f a s t , u s u a l l yw i t h i nt h r e e o ff o u ri t e r a t i o n s .S t e p( f ) ,t h ec o r r e c t i o ns t e p , a t t h e same t i m e t h e e r r o r i s b a s e d upon t h e i s beingcomputed, Newton-Raphson procedure: the derivative of t h e e r r o r w i t h r e s p e c t t o t h e i n d e p e n d e n t v a r i a b l e a l s of o r m e d ,u s i n gd e r i v a t i v e s 34 is of t h e t a b u l a r q u a n t i t i e s Y and Tw ( o r B b ) i s t h e n computed d i r e c t l y : The c o r r e c t i o n i n t h e i n d e p e n d e n t v a r i a b l e or E . 4.5 AVOIDING THE SURFACE STATE SOLUTION The c o n v e c t i v e b o u n d a r y c o n d i t i o n d e s c r i b e d so f a r i n t h i s s e c t i o n d o e sn o ti n c l u d e a l l p r o b l e m so fi n t e r e s t .C o n s e q u e n t l yt h es u r f a c eb o u n d a r y c o n d i t i o nt r e a t m e n th a sb e e nf o r m u l a t e d as a t r i p l e o p t i o n . The f i r s t o f t h e s e( c a l l e dO p t i o n i nS e c t i o n 4.1-4.4 1) c o n s t i t u t e s the c o n v e c t i v eb o u n d a r yc o n d i t i o nd e s c r i b e d above. A s e c o n do p t i o n( O p t i o n 2) m e r e l yr e f e r st oi n p u t time-dependentquantitiesforuserassignedvaluesofsurfacetemperature a n dr e c e s s i o n rate. A t h i r do p t i o n( O p t i o n 3 ) r e t a i n st h ee n e r g yb a l a n c e aspectsofOption 1, b u t e l i m i n a t e s t h e c o n v e c t i v e h e a T i n g t e r m s a n d d o e s n o ta l l o ws u r f a c er e c e s s i o n .T h i so p t i o n i s u s e f u lf o r "cool-down'' o r "soak - o u t " p o r t i o n s o f p r o b l e m s , a n d c a n a l s o f l u xp r o b l e m si ng e n e r a l( t h r o u g ht h eq r a d be u s e d f o r a s s i g n e d h e a t in v a r i a b l e ) i f s u r f a c e r e c e s s i o n i s n e g l i g i b l e . The h e a t i n g o p t i o n may bechanged a tt h eu s e r ' sc h o i c ea t any t i m e f o r a g i v e n column o r s u r f a c e s t a t i o n , and t h e h e a t i n g o p t i o n a s s i g n ments may v a r y i n anymanneraround t h e body a t a given t i m e . 35 SECTION 5 SURFACE STATE 5.1 CALCULATION - PHASE I V INTRODUCTION Section 4 abovehasdescribedtheconvectivesurfaceenergybalance the u s e made i n t h i s b o u n d a r y c o n d i t i o n calculationofcomputersubroutinesforevaluatingvariousheatedsurface b o u n d a r yc o n d i t i o na n dd e s c r i b e d thermochemical s t a t e v a r i a b l e s , c h i e f l y e n t h a l p y a n d m o l e c u l a r c o m p o s i t i o n . it would be possible to avoid the use of c h e m i s t r y r o u t i n e s f o r t h i s c a l c u l a t i o n a n d i n s t e a d r e l y on e n t h a l p y correlationsandthelike,butforthepresentcoupledcode it w a s s t r o n g l y d e s i r e d t o have a g e n e r a l a b l a t i o n p r o g r a m w h i c h c o u l d t r e a t a n y n o n - c h a r r i n g a b l a t i n gm a t e r i a le x p o s e dt o any a r b i t r a r ye n v i r o n m e n t .F o r t u n a t e l y , suitablethermochemistryroutinesalreadyexistedwhichcouldbecoupledto For a restricted range of problems t h ei n - d e p t hr e s p o n s ec o d e :t h ep r o g r a mr e p o r t e dh e r ed i dn o td e v e l o pt h e s e c h e m i s t r yc o d e s .S i n c et h e s ec o d e sa r ef u l l yd e s c r i b e de l s e w h e r e( R e f e r e n c e s 4 and 5 ) t h i s s e c t i o n physicalproblemstreated w i l l seekonlytogive a g e n e r a lo v e s v i e wo ft h e by t h e c h e m i s t r y r o u t i n e s . as a package are r e f e r r e d t o a s ACE, a f t e r AerothermChemicalEquilibriumcode.This i s an i n a p p r o p r i a t e acronym s i n c e (1) t h er o u t i n e sc o n s i d e rn o n - e q u i l i b r i u mc h e m i s t r y as w e l l a s e q u i l i b r i u m , and ( 2 ) i n t h e c o u p l e dc o d en or o u t i n eh a st h e name ACE; t h ec o n t r o l l i n gr o u t i n e i s c a l l e d EQUIL. F o rh i s t o r i c a lr e a s o n s ,h o w e v e r , w e shallcontinuetocallthethermochemistrypackage ACE. The t h e r m o c h e m i c a lr o u t i n e s USES OF 5.2 ACE I N MACABRE The ACE package i s c a l l e d upon f o r a number o f c o m p u t a t i o n a l t a s k s i nt h ec o u p l e dc o d e . i s i ng e n e r a t i n gs u r f a c e I t s mainuse,ofcourse, thermochemical s t a t e s o l u t i o n s . A s w i l l b ed i s c u s s e db e l o w ,t h e s ea r eo p e n s y s t e mc a l c u l a t i o n si n v o l v i n gi n t e r a c t i o n sb e t w e e ne d g eg a sa n dt h em a i n a b l a t i n g material. terms l i k e hw The s u r f a c ee n e r g ye q u a t i o n s and CKi gas edge h'liy e or CZ; hp (33) and (34) i n c l u d e , h o w e v e r , w h i c ihn d i c a t easn o t h e r . e n e e d e dt h e r m o c h e m i c a lc a l c u l a t i o n :t h em o l e c u l a r make-up o ft h eb o u n d a r y l a y e re d g eg a sm u s tb ei d e n t i f i e da n dc e r t a i np r o p e r t i e se v a l u a t e d( s u c h as Z ? ) 1 . Each c a l l t o ' ACE f o r a s u r f a c e s t a t e s o l u t i o n i n f a c t i n v o l v e s t h e follo8ing computations : 36 a closed system equilibrium calculation on edge gas a l o n e t o d e t e r m i n e Ki and Z? ; t h i sc a l c u l a t i o n 'e done a t a s s i g n e d p r e s g u r e a n d r e c o v e r y e n t h a l p y * . is anopensystemablationcalculation a t assignedpressure, peueCM ( f o r k i n e t i c e f f e c t s ) , and e i t h e r B& o r Tw (depending on l o c a t i o n i n i n d e p e n d e n t v a r i a b l e a r r a y ) : this calculation may i n c l u d e k i n e t i c e f f e c t s a f r o z eenv a l u a t i oonf t e m p e r a t u r e TW. hw and edge gas , Z 2; h p Wedge gas asdescrlbedinSections The q u a n t i t i e s h e from ACE andemployed at h? CZ* 'e , hw, and h? C 2; are o b t a i n e d W 4.3 and 4 . 4 above. (The ACE r o u t i n e may also be used by t h e p r e s s u r e d i s t r i b u t i o n and t r a n s f e r c o e f f i c i e n t r o u t i n e s t o g e n e r a t e c e r t a i n t h e r m o d y n a m i c p r o p e r t i e s and t r a n s p o r tp r o p e r t i e s of t h eb o u n d a r yl a y e rg a s .T h e s ep r o p e r t i e s are s t o r e d for f u t u r e u s e i n a M o l l i e r C h a r t (as w i l l b e d e s c r i b e d i n more detailinSections 8 and 9 below) i n which P andhe a r et h ei n d e p e n d e n t variables;therequired ACE c a l c u l a t i o n s a r e t h e r e f o r e c l o s e d s y s t e m e q u i l i b r i u m c a l c u l a t i o n s of t h e e d g e g a s a t a s s i g n e d p r e s s u r e a n d e n t h a l p y . ) The f o l l o w i n g s e c t i o n s d e s c r i b e t h e manner i n which c a l c u l a t i o n s a r e done i n t h e ACE package. 5.3 P H Y S I C A L PROBLEMS TREATED BY THE ACE PACKAGE The ACE p a c k a g e o f s u b r o u t i n e s t.reats a widerange of thermoc h e m i c a lp r o b l e m s .S i n c et h e s ep r o b l e m sa r ec o m p l e x ,t h ec a p a b i l i t i e so f t h e ACE p a c k a g e c a n p e r h a p s b e s t b e i l l u s t r a t e d by a s e t o f d e s c r i p t i o n s o fi n c r e a s i n gg e n e r a l i t y . The f o l l o w i n gs e c t i o n sd e s c r i b es u c h a graded seriesofproblems. 5.3.1 Closed System Eiqnu i l i b r i u m Closedsystems are d e f i n e d as t h o s e f o r which t h e r e l a t i v e amounts o fc h e m i c a le l e m e n t sa r ep r e s p e c i f i e d . The e d g eg a ss t a t ec a l c u l a t i o n sa r eo ft h i st y p e . Consider K chemicalelements, Good argumentscanbe made f o r u s i n g p r o v e si n c o n v e n i e n ti nt h e code.The any c a s e . Nk, introducedinto a previously the s t a t i c e n t h a l p y i n s t e a d , b u t t h i s d i f f e r e n c ei nr e s u l t s i s s m a l li n 37 e v a c u a t e dc o n t a i n e r .I ng e n e r a l ,t h e s ee l e m e n t s number o fc h e m i c a ls p e c i e s * I f enough t i m e h a s e l a p s e d , Ni (gasphase)and w i l l i n t e r a c t t o form a N R ( c o n d e n s e dp h a s e s ) . s o t h a t thermodynamicandchemicalequilibrium is e s t a b l i s h e d , t h e t h e r m o d y n a m i c s t a t e o f t h e s y s t e m , i n c l u d i n g t h e r e l a t i v e i s completelydetermined i f two i n d e a m o u n t so fc h e m i c a ls p e c i e sp r e s e n t , pendentthermodynamicvariablesare known ** . T h i sc o n d i t i o n mathematicallybyexaminingthegoverningequationsforsuch andshowing that the number of i n d e p e n d e n t e q u a t i o n s may b es t a t e d a system, is equal to the number of unknown q u a n t i t i e s . of t h eg a s e o u sc h e m i c a ls p e c i e s as f o l l o w s : R e l a t i o n se x p r e s s i n gt h ef o r m a t i o n from t h e g a s e o u s c h e m i c a l e l e m e n t s may b e w r i t t e n S i m i l a r l y ,f o r m a t i o no fc o n d e n s e dp h a s es p e c i e sf r o mt h eg a s e o u se l e m e n t s is written: K 'kll (38) Nk K= 1 Intheabove, o fs p e c i e s Cki number of atoms of e l e m e n t k i n a molecule representsthe R (condensed). i ( g a s )o rs p e c i e s If t h e g a s p h a s e s p e c i e s a r e assumed t o i n d i v i d u a l l y b e h a v e a s t h e r m a l l yp e r f e c tg a s e s ,t h e nt h ee q u i l i b r i u mr e l a t i o nc o r r e s p o n d i n gt o r e a c t i o n( 3 7 ) is k=l I n Pi - Cki I n Pk = I n K k= 1 where Pk d e n o t e sp a r t i a lp r e s s u r ea n d f o rt h ef o r m a t i o nr e a c t i o n( 3 7 ) ** 38 K Pi Pi (T) ( T ) i s t h ee q u i l i b r i u mc o n s t a n t of s p e c i e s Ni. Foreachcandidatecondensed " C h e m i c a ls p e c i e s "a su s e dh e r ei n c l u d e sm o l e c u l a r ,a t o m i c ,i o n i c , electronspecies. Duhem's Theorem. (39) and p h a s es p e c i e s K k=l where = i n d i c a t e st h ee x i s t e n c eo ft h ec o n d e n s e dp h a s es p e c i e s N E i ne q u i l i b r i u mw i t hg a sp h a s es p e c i e s ,a n d < NE w i l l not i n d i c a t e st h a tt h ec o n d e n s e dp h a s es p e c i e s bepresentinequilibrium. F o re a c hc h e m i c a le l e m e n ti n t r o d u c e di n t ot h es y s t e m ,t h ec o n s e r v a t i o n of atoms d i c t a t e s that t h e amount ofanyelement k i nt h eg a s andcond e n s e dp h a s e s( r e g a r d l e s so fm o l e c u l a rc o n f i g u r a t i o n )m u s t sum t o t h e t o t a l amount o fe l e m e n t k i nt h es y s t e m .M a t h e m a t i c a l l y ,t h i s may b ew r i t t e n , f o re a c he l e m e n t k , as Mass f r a c t i o n of e l e m e n t k inpu to ht e = sys t e m where nf, x I i=l L s i p i + Ec ~ ~ x Q 11=1 74 i s a c o m p o s i t es y s t e mm o l e c u l a rw e i g h t *d e f i n e d andwhere Xi i s amole xi by f r a c t i o n of condensedphasespecies = m o l e c u l e s of c o n d e n s e d s p e c i e s totalgasphasemolecules i This is themolecularweightappropriate i f c o n d e n s e dp h a s e sa r ep r e s e n t . R d e f i n e da s 9. to theidealgasequation (43) of s t a t e 39 , In addition, for the gas phase species thereexiststherequirementthat system pressure sum t o t h e t o t a l partial pressures must the = P , Mixture thermodynamic properties relatedtothespeciesconcentrations such as specific enthalpy , are b ye q u a t i o n so ft h ef o r m Consider now t h e number o fi n d e p e n d e n te q u a t i o n sf o rt h es y s t e m . is e q u a lt ot h e number o fg a sp h a s ee q u i l i b r i u mr e l a t i o n s( 3 9 ) g a sp h a s es p e c i e s aretrivial I minus t h e number o fe l e m e n t s when i = k ) f o re a c ho ft h e . I na d d i t i o n ,t h e r ee x i s t s K The number of ( b e c a u s eE q u a t i o n s( 3 9 ) (40) a r e l a t i o ns u c ha s c a n d i d a t ec o n d e n s e dp h a s es p e c i e si nt h es y s t e m .N o t e L t h a tt h es y s t e mt e m p e r a t u r e i s c o n t a i n e di m p l i c i t l yi nE q u a t i o n s( 3 9 )a n d t h r o u g ht h et e m p e r a t u r ed e p e n d e n c eo ft h ee q u i l i b r i u mc o n s t a n t s .T h e r e K c o n s e r v a t i o no fe l e m e n t sE q u a t i o n s i n t r o d u c e di n t ot h es y s t e m . t o t h es y s t e mp r e s s u r e (44) (41) , one f o re a c ha t o m i ce l e m e n t The r e q u i r e m e n t t h a t t h e p a r t i a l p r e s s u r e s c o n t r i b u t e so n ea d d i t i o n a le q u a t i o n . (45) sum For any a d d i t i o n a l thermodynamic p r o p e r t i e s o f t h e m i x t u r e ( e n t h a l p y , e n t r o p y , t h e r ee x i s t se q u a t i o n ss u c ha s (40) are . etc.), of t h e I s p e c i e si nt h eg a sp h a s ea r e Considernextthevariablesappropriatetothisformulation problem. The r e l a t i v ec o n c e n t r a t i o n so ft h e g i v e n by t h e P i ' s and t h e amountsof s p e c i e sa r eg i v e n by X L (most o r a l l o f the l a t i o n ,t h ec o m p o s i t es y s t e mm o l e c u l a rw e i g h t a r e oneeach L candidatecondensedphase which may b e z e r o ) . I n t h i s ,% of t h em i x t u r et h e r m o d y n a m i cv a r i a b l e s number o f v a r i a b l e s a n d a v a i l a b l e i n d e p e n d e n t e q u a t i o n s 40 formu- i s also a v a r i a b l e .T h e r e T , P, h , s , e t c . The may besummarizedas VARIABLES NO. OF SUCH VARIABLES EQUATION NUMBER TO. O F SUCH EQUATIONS I (39) I - K L L K (44) 1 1 1 of the type n 1I total variables T h u s ,t h e r ea r e (45) total equations I+L+n+3 two less e q u a t i o n st h a n dynamic s t a t e o f t h e s y s t e m I+L+n+l there a r e v a r i a b l e s , i n d e p e n d e n tv a r i a b l e sa r es p e c i f i e d( e . g . , e l e m e n t a lc o m p o s i t i o n ,t h e nc l o s u r e n and s o i f two P and T ) i n a d d i t i o n t o t h e i s obtainedandthechemicalandthermo- may, i n p r i n c i p a l , b e d e t e r m i n e d . The ACE p r o g r a mp e r f o r m st h i sd e t e r m i n a t i o n .T h a t the e q u i l i b r i u m t h e r m o d y n a m i c a n d c h e m i c a l s t a t e o f i s , t od e t e r m i n e any . c l o s e d s y s t e m , o n e needsonly tofurnishthe ACE p r o g r a m w i t h t h e e l e m e n t a l c o m p o s i t i o n , t h e * candidategaseousandcondensedphasespecies, andtwo thermodynamicproperties. One o f t h e s e p r o p e r t i e s m u s tb ep r e s s u r e ,a n dt h eo t h e r e i t h e rt e m p e r a t u r e ,e n t h a l p y ,o re n t r o p y . program w i l l c a l c u l a t e a n do u t p u tt h e may b e Given t h i si n f o r m a t i o n ,t h e ACE mole f r a c t i o n so fe a c hc a n d i d a t e s p e c i e s , the t e m p c L a t u r e ,e n t h a l p y ,e n t r o p y ,d e n s i t y ,e f f e c t i v em o l e c u l a r w e i g h t ,e q u i l i b r i v m and f r o z e n s p e c i f i c h e a t s , t h e i s e n t r o p i c e x p o n e n t , a n d a few other quantities of potential interest. A= a l r e a d yn o t e d , a n de n t h a l p ys o l u t i o n s . opensystemcalculation MACABRE o n l y c a l l s The e d g e s t a t e upon ACE f o r a s s i g n e d p r e s s u r e i s s a v e dm o m e n t a r i l yw h i l et h e i s performed. *The ACE p r o g r a m m u s t b e s u p p l i e d c e r t a i n b a s i c t h e r m o d y n a m i c d a t a ( d e s c r i b e d i n R e f e r e n c e s 2 and 5) f o r e a c h c a n d i d a t e s p e c i e s t o b e c o n s i d e r e d i n a given i s c o n t a i n e d i n a t h r e ec a r d s e t , one s e t f o r e a c hs p e c i e s . s y s t e m .T h i sd a t a A c e r t a i n amount ofjudgement is r e q u i r e d o n t h e p a r t o f t h e u s e r r e l a t i v e t o w h i c hc a n d i d a t es p e c i e ss h o u l db ei n c l u d e di n a g i v e ns y s t e m .F r e q u e n t l y t h i s judgement i s a v o i d e d b y s i m p l y i n p u t t i n g d a t a f o r a l l s p e c i e s c o n t a i n i n g c o m b i n a t i o n so ft h ei n p u te l e m e n t s . 41 5.3.2 Open Systems i n E q u i l i b r i u m - S i m p l i f i e dC a s e The b a s i c t h e o r y u n d e r l y i n g t h e t r e a t m e n t o f o p e n s y s t e m s may b e s t b e i l l u s t r a t e d by e x a m i n i n g t h e e q u a t i o n s e x p r e s s i n g t h e c o n s e r v a t i o n e l e m e n t s and energy a t t h e a b l a t i n g s u r f a c e . I f t h e b o u n d a r y l a y e r no m a t e r i a l i s removedfrom t e r i z e d by e q u a l d i f f u s i o n c o e f f i c i e n t s , a n d i f the surface in ( i.e. a condensed phase r e m o v a l ) ,t h e nt h e s ee q u a t i o n st a k e , no mechanical erosion or liquid layer on a p a r t i c u l a r l y s i m p l e form f o r e q u i w i l l b ec o n s i d e r e di nt h i ss e c t i o n l i b r i u ms y s t e m s .T h i ss i m p l i f i e ds i t u a t i o n and t h e s e c o n s i d e r a t i o n s of chemical is charae- w i l l beextendedto more g e n e r a l i z e d c a s e s i n Sections 5.3.3 t o 5.3.5. of chemicalelements ( k ) e n t e r i n ga n d l e a v i n g a c o n t r o ls u r f a c ea f f i x e dt ot h ea b l a t i n gs u r f a c e . The s o l i d m a t e r i a l may b ev i s u a l i z e da s moving i n t o t h i s s u r f a c e a t a r a t e 6 . I f it i s assumed C o n s i d e rf i r s tt h ef l u x e s t h a t no m a t e r i a l i s b e i n g removed i n a c o n d e n s e d p h a s e , t h e n t h e s u r f a c e and t h e f l u x e s o f t h e kth chemicalelement Figure 1 6 . may b e i l l u s t r a t e d A b l a t i n gS u r f a c e Terms s u p e r s c r i p t e d by a t i l d e ( - ) Mass Balance representthetotal o fe l e m e n tk i, n d e p e n d e n to fm o l e c u l a rc o n f i g u r a t i o n . as mass f r a c t i o n o r f l u x Thus I 'kW = 'k?iw i=l 42 I I ‘kW = ukijiw i=l where k p e r t a i n st oe l e m e n tk , mass f r a c t i o no e f lement k (47) i s the i p e r t a i n st os p e c i e s i , and o’ki i ns p e c i e s i. F l u x e o s ef l e m e n t k away from thesurfaceconsistofboundarylayerdiffusionandgrossmotionofthefluid a d j a c e n tt ot h es u r f a c ed u et ot h ei n j e c t i o nf l u x AC. From the a b o v es k e t c h , it i s s e e n t h a t c o n s e r v a t i o n o f c h e m i c a l e l e m e n t sr e q u i r e st h a t jkw + ( P V ) K = kw Summing E q u a t i o n ( 4 8 ) o v e ra l le l e m e n t s equation(forthecase kc k y i e l d st h et o t a l mass c o n t i n u i t y when t h e r e i s nocondensedphasematerialremoval) (Pv), An r;l K = (49) mC i m p o r t a n tf u n d a m e n t a lo ft h ep r e s e n tm a t h e m a t i c a lm o d e l i n go ft h ea b l a t i o n p r o c e s s i s t h ee x p r e s s i o no ft h ed i f f u s i v eh e a t a transfercoefficientformulation. In t h i fs o r m u l a t i o nt,h e diffusional and mass f l u x e s i n term -j terms o f is written kW U t i l i z i n gt h i st r a n s f e rc o e f f i c i e n tf o r m u l a t i o n f l u xi nt h ee l e m e n t a lb a l a n c eE q u a t i o n ( 5 0 ) f o rt h ed i f f u s i o n (48) yields 43 .1 ..1 .". ... ,.., -_". I n terms of d i m e n s i o n l e s s a b l a t i o n r a t e (51) for t h e t o t a l B& a n ds o l v i n g mass f r a c t i o n of e l e m e n t k a t t h e w a l l y i e l d s ( f o r e q u a l d i f f u s i o n c o e f f i c i e n t s and no condensedphase material removal) K, Given t h e r e l a t i v e chemicalandthermodynamic kc = a m o u n t so fc h e m i c a le l e m e n t ss p e c i f i e db y( 5 2 ) I the t o theablatingsurface state ofthegasesadjacent may be c a l c u l a t e d f r o m e q u i l i b r i u m r e l a t i o n s . Sincethegases are i n e q u i l i b r i u m w i t h t h e a b l a t i n g s u r f a c e (53) k=l i f 9. r e p r e s e n t st h es u r f a c es p e c i e s ,a n d (54) I n t h ep r e s e n tf o r m u l a t i o n , (54) may be t h o u g h to fa sb e i n gu s e dt oi d e n t i f ys u r f a c es p e c i e s R for f o r a l l o t h e rc a n d i d a t ec o n d e n s e dp h a s es p e c i e s . which(53)applies. The e q u i l i b r i u mr e l a t i o n s f o r g a sp h a s es p e c i e s the form K and t h e p a r t i a l p r e s s u r e s mustobey the relation I & = P i-1 I The e l e m e n t a l mass f r a c t i o n sa d j a c e n tt ot h es u r f a c e l related to the species partial pressures , Pi, K I k" of ( 5 2 )a r e by r e l a t i o n s s u c h a s i s of Thus, i f P i s known and BL i s s p e c i f i e d ( t h i s may b e v a r i e d p a r a m e t r i c a l l y , and as w i l l b e d i s c u s s e d s u b s e q u e n t l y ) i f Tw and Pi number o f unknowns a n a e q u a t i o n s a v a i l a b l e I I. NO. OF SUCH UNKNOWNS UNKNOWNS . . ~~ % .nc Tw - EQUATION :QUATION NO. . I T h u s ,c l o s u r e NO. O SF UCH EQUATIONS AVAILABLE . I (55) I-K K K 1 (53) 1 1 (52) K (56) 1 1 (57) K 1 Unknowns are unknowns, t h e n t h e may b e summarized as -1'0" -11 I+K+2 T 0 Eqquuaa tt ii o on n ss 1 t ' Z K ++ 22 II ++ K - Y i s o b t a i n e da n d ,i np r i n c i p a l ,t h et e m p e r a t u r e and t o t h e s u r f a c e may bedetermined From t h ep r e s s u r e ,t e m p e r a t u r e , andchemical c h e m i c a lc o m p o s i t i o no ft h eg a s e sa d j a c e n t P and EA are s p e c i f i e d . c o m p o s i t i o n ,t h ec a l c u l a t i o no fo t h e rt h e r m o d y a n m i cp r o p e r t i e s( e n t h a l p y , etc.) i s straightfornard. if The ACE p a c k a g e o f p r o g r a m s h a s t h e c a p a b i l i t y t o d e t e r m i n et h e s t a t e o ft h eg a s e sa d j a c e n tt oa na b l a t i n gs u r f a c e i n a f a s h i o ns i m i l a rt ot h a td i s c u s s e dh e r e .I nt h ec o u p l e d MACABRE c o d e ,t h e ACE group i s t y p i c a l l yc a l l e dw i t hp r e s s u r e P known and B ' C a s s i g n e da sa ni n d e p e n d e n tv a r i a b l e by t h eu s e r . The t e m p e r a t u r e T and chemicalandthermodynamic W o t h e r s t a t e d a t a are t h e nd e t e r m i n e d .A l t e r n a t i v e l y , ( i nt h el o w e rt e m p e r a t u r ep a r t so ft h et a b l e s , 4.2,2.2 above)and B' Tw may b ea s s i g n e d as d i s c u s s e d i n S e c t i o n i s determined. C 5.3.3 O p e n S y s t e m s i nE q u i l i b r i u m - U n e q u a lD i f f u s i o nC o e f f i c i e n t s w a s l i m i t e d t o opensystems material i n t h e c o n d e n s e dp h a s e .T h e s es i m p l i f i c a t i o n s were made i n o r d e r t o r e n d e r t h eb a s i ct h e o r y e a s i e r t o e x p l a i n .W h i l et h i ss i m p l em o d e l i s reasonably The d i s c u s s i o n o f t h e p r e v i o u s s e c t i o n withequalspeciesdiffusioncoefficientsandnoremovalofsurface accuratefor many a b l a t i o n s i t u a t i o n s , t h e s e a s s u m p t i o n s a r e i n a p p r o p r i a t e f o ro t h e r s .F o rt h i sr e a s o n ,c a l c u l a t i o n sp e r f o r m e d by t h e ACE package ( a n d h e n c e t h e MACABRE program) are n o t r e s t r i c t e d t o any of t h e s e s i m p l i f i cations. The e x t e n , s i o n t o u n e q u a ld i f f u s i o nc o e f f i c i e n t s w i l l b ed i s c u s s e d first, in this section. 45 The u n e q u a l d i f f u s i o n c o e f f i c i e n t model h a s a l r e a d y b e e n b r i e f l y i n t r o - 4.3 above , i n whichEquations (33) and (34) may becompared. mass f r a c t i o n d u r i n g f o r c e s Ki * e w i t hm o d i f i e dd r i v i n gf o r c e s Zi A s i m i l a rm o d i f i c a t i o nm u s tb e made i n t h e ACE t r e a t m e n t . The d i f f u s i G e flu ofelement k a t the w a l l i s now modeled by duced i n S e c t i o n I n s i m p l i s t i c terms, t h e model r e p l a c e s t h e . * 'kw e e M :z; = p u c -2;) e In (58), i s , i n - e f f e c t , a w e i g h t e da v e r a g eo ft h em o l e Z; of element k . The Z; a r ed e f i n e d andmass fractions by ~ i=l (59) ZTK1-Y z; = I (where y Ref. 10 1 P 2/3, see Ki/Fi i=l where t h e Fi a r e d i f f u s i o n f a c t o r s d e f i n e d by t h e f o l l o w i n g r e l a t i o n f o r thebinarydiffusioncoefficients where 5 i s a c o n s t a n tf o r weakly on t e m p e r a t u r e . a g i v e nt e m p e r a t u r ea n dp r e s s u r ea n dt h e The -0 tj mustobey ( 6 2 ) i no r d e rf o rt h e l a y e rs p e c i e sd i f f u s i o ne q u a t l o n st or e d u c et o b ei n f e r r e d by s i m i l a r i t ya r g u m e n t s .R e f e r e n c e binarydiffusioncoefficientsfor c o r r e l a t e d by ( 6 2 ) . F i is ( 5 8 ) can 11 d e m o n s t r a t e st h a tt h e a v a r i e t yo fc h e m i c a ls y s t e m sa r ea c c u r a t e l y T h i sr e f e r e n c ea l s o r e l a t i o ne q u a t i o nf o rt h e 46 a formfromwhich F depend i boundary shows t h a t a r e a s o n a b l y good cor- when 5 i s t a k e n as t h e s e l f - d i f f u s i o n c o e f f i c i e n t o f 02. A d d i t i o n a l d i s cussionrelativetothisunequaldiffusioncoefficientformulation i s cont a i n e d i n Appendix A. C o n s i d e r a t i o no fu n e q u a ld i f f u s i o nc o e f f i c i e n t sa f f e c t st h es u r f a c e elementalbalancerelationshipswhichareneededto.determinetheequilibrium thermochemical s t a t e a tt h es u r f a c e .S u b s t i t u t i n g( 5 8 )i n t o( 4 8 )y i e l d s and t h e "unknowns" h e r e a r e Kk I and Z; , each of which may beexpressed partiaY pressur& i nt e r m so ft h es p e c i e s 4 and I where 7n 9 defi*ned as is t h e mean m o l e c u l a rw e i g h to ft h eg a sp h a s e and F is a mean Fy i i=l 47 Substitutingtheseexpressionsinto (64) and utilizing the definition of B'C' yields anexpressionforthespeciespartialpressuresatthesurfacein terms of q u a n t i t i e s a t t h e Note t h a t( 6 8 )r e d u c e st o boundary l a y e r e d g e a n d i n t h e m a t e r i a l ( 5 2 ) when t h ed i f f u s i o nc o e f f i c i e n t sa r ee q u a l . When p e r f o r m i n g u n e q u a l d i f f u s i o n c o e f f i c i e n t o p e n s y s t e m c a l c u l a t i o n s , t h e ACE program u t i l i z e s ( 6 8 ) r a t h e r t h a n , e q u a t i o n s .O t h e rt h a nt h i s d i s c u s s e di nS e c t i o n 5.3.2. (52) a st h ee l e m e n t a l t h es o l u t i o np h i l o s o p h y mass b a l a n c e is e s s e n t i a l l ya s The d i f f u s i o n f a c t o r s u t i l i z e d i n t h e ACE pro- gram may b e c a l c u l a t e d i n t h r e e w a y s , a t t h e u s e r ' s o p t i o n a .d i f f u s i o nf a c t o r s e a c hs p e c i e s b. Fi may b ei n p u ti n d i v i d u a l l yf o r i d i f f u s i o nf a c t o r s may b ec a l c u l a t e da c c o r d i n gt o( 6 3 ) withtheuserspecifyingthereferencemolecular weight, c. r e f , and t h ee x p o n e n t E i ft h eu s e rd o e sn o t h i n gs p e c i a l ,t h ep r o g r a m ref = 23.4 and c a l c u l a t e Fi a c c o r d i n gt o( 6 3 )w i t h E will = 0.431. The a c t u a l p r o g r a m i n p u t f o r t h e s e a l t e r n a t i v e s i s discussedinReference 2. It should also be pointed out that the diffusion factors have an effect on theothertransportpropertiescalculated briefly discussed in by t h e ACE program,andthese Appendix A. F o ru n e q u a ld i f f u s i o nc o e f f i c i e n t s ,t h et r a n s f e rc o e f f i c i e n tf o r m u l a - is E q u a t i o n( 3 4 ) .C o n s i s t e n tw i t h( 3 4 ) , t i o nf o rt h es u r f a c ee n e r g ye q u a t i o n f o ru n e q u a ld i f f u s i o nc o e f f i c i e n tp r o b l e m s ,t h es u r f a c et h e r m o c h e m i s t r y q u a n t i t i e s p r e p a r e d by t h e ACE p r o g r a m i n c l u d e t h e q u a n t i t y 2 : 2 i=l hiTW W i na d d i t i o nt ot h eq u a n t i t i e sp r e v i o u s l yd i s c u s s e d . 48 N o t ea g a i n ,t h a tf o r are equaldiffusioncoefficients Similarly, the quantities and frozen hw edge gas are computed a f t e r r e f e r e n c e t o t h e p r e v i o u s a u t o m a t i c a l l y computedand storededgesolution. 5.3.4 all the Open Systems In Equilibrium M a t e r i a l Removal -. With Condensed Phase Surface Considerationstothispointhavebeenbased on t h e a s s u m p t i o n t h a t mass l e a v i n g a n a b l a t i n g s u r f a c e i s i nt h eg a sp h a s e and t h i s flux hasbeendenoted by (pv),.However, f o r many a b l a t i o n s i t u a t i o n s o f i n t e r e s t , some o f t h e m a t e r i a l l e a v i n g t h e s u r f a c e i s i n a condensedphase - e-g., l i q u i d l a y e r removal o r" m e c h a n i c a la b l a t i o n " . The condensedphaseremoval a f f e c t sb o t ht h es u r f a c e m a s s and e n e r g yb a l a n c e s( R e f e r e n c e 1 0 ) . For t h e s u r f a c ee l e m e n t a lb a l a n c e ,t h e r e i s an a d d i t i o n a l f l u x term L 'L'k, L=1 l e a v i n gt h es u r f a c e .T h i sr e s u l t si nt h es u r f a c ee l e m e n t a lb a l a n c e( i n terms o ft h es p e c i e sp a r t i a lp r e s s u r e s - analogous t o( 6 8 )h a v i n gt h ef o r m 49 element k . foreachatomic Summing o v e r a l l e l e m e n t sy i e l d st h eo b v i o u sn e t mass c o n t i n u i t y r e l a t i o n S i m i l a r l y ,t h e r e is an a d d i t i o n a le n e r g yf l u x L term II=1 l e a v i n gt h es u r f a c e . The s u r f a c ee n e r g yb a l a n c ee q u a t i o n becomes ( a n a l o g o u s t o (34) r peuecH ( H r- hw) f r o z e n (2; - to this EA h c W The condensedphasesurfacematerialremoval p o r a t e di nt h e + ) hTW i 2; edgegas model c u r r e n t l y i n c o r - ACE program i s termed a “ f a i lt e m p e r a t u r em o d e l “ .A c c o r d i n g model, a “ f a i l t e m p e r a t u r e ” , T f a i l - & , L. candidatecondensedphasespecies f r o mt h es u r f a c e( i n may b e p r e - a s s i g n e d t o T h i ss p e c i e s a condensedphase) a t a rate w i l l t h e nb e ma if any removed t h es u r f a c et e m p e r a - t u r e is g r e a t e r t h a n , o r e q u a l t o , t h e f a i l t e m p e r a t u r e f o r t h i s s p e c i e s . A d d i t i o n a l comments r e l a t i v e t o t h e c o m p u t a t i o n a l a s p e c t s briefdiscussion of t h e f a i l t e m p e r a t u r e c o n c e p t . The f a i l t e m p e r a t u r e i.e., w i l l follow a is the temperature at which t h e m a t e r i a l “ f a i l s “, w i l l not remain affixed to which a l i q u i d l a y e r i s u n s t a b l e , t h e thetemperatureabovewhichthematerial t h es u r f a c e .F o rm e l t i n gs o l i d sf o r f a i lt e m p e r a t u r e is s i m p l yt h e m e l t t e m p e r a t u r e .F o rt h e n o c h e m i c a l l y erodingsolids,thefailtemperaturemightbe s e t a t some t e m p e r a t u r e c h a r a c t e r i s t i c of t h e m a t e r i a l y i e l d t e m p e r a t u r e a t t h e a p p r o p r i a t e e x t e r n a l l o a d i n g .A l s o ,t h e ACE programhas f a i lt e m p e r a t u r ef o ra l l temperaturemightcorrespond When f a i l i n g o c c u r s , a maximum maximum f a i l to the fail temperature of the substrate. an a d d i t i o n a l e q u i l i b r i u m r e l a t i o n able f o re a c hc o n d e n s e dp h a s es p e c i e s 50 a p r o v i s i o nf o ri n p u t t i n g c o n d e n s e ds p e c i e s .P h y s i c a l l y ,t h i s II f o r which mII > 0 . is avail- A l s ot h es u r f a c e temperature must be consistent with the assignment fail temperatures: ~w ~w ' ' *fail 'fail for the surface species f o r any s p e c i e s f o r w h i c h ma. > 0 Theseadditionalequilibriumrelations,surfacetemperaturerestrictions, andassignedfailtemperaturesaresufficienttosolvefortheadditional mass b a l a n c e r e l a t i o n s . unknowns rhk i n t r o d u c e di n t oe l e m e n t a l Ha. , thermochemicaldata. s o l i d , and i f Tfail < I f Tfail Tmelt, The e n t h a l p y , removed is d e t e r m i n e df r o mt h es p e c i e s o ft h ec o n d e n s e dp h a s eb e i n g Ha. is t a k e na st h ee n t h a l p yo ft h e Tmelt, Ha. i s t a k e na st h ee n t h a l p yo ft h el i q u i d . The conventionforassociatingfailtemperatureswitheachorallcandidate condensedphasespecies w i l l bespecifiedinReference Incommunicatingback to MACABRE, terms ( p v )wHw and 2 * . t h e ACE packageamalgamates the L by s u i t a b l ym o d i f y i n gt h er e t u r n e dv a l u eo f quantity Y asdescribedinSection Hw. Thus FILL c o n s t r u c t s t h e 4.3 without specific reference to failing events. 5.3.5 Non-Equilibrium E f f e c t s i n Open S y s t e mC a l c u l a t i o n To calculatetheequilibriumstate of a c h e m i c a ls y s t e m ,i n f o r m a t i o n r e l a t i v et oa l lc h e m i c a lr e a c t i o n s i s n o tn e e d e d .T h i sf a c tp e r m i t s a significantsimplificationinthe ofthe p r o b l e mf o r m u l a t i o n ,a n dt h o s ef e a t u r e s ACE p a c k a g e d i s c u s s e d t h u s f a r t a k e a d v a n t a g e o f t h e s e s i m p l i f i c a - t i o n s .F o r some s y s t e m so fi n t e r e s t ,h o w e v e r ,t h e etics may n o tb en e g l e c t e d . A g e n e r a ls o l u t i o no f reactionkineticseffectsareimportant l e a s t ' t w or e a s o n s :a )t h e r e e f f e c t s o fr e a c t i o nk i n complexproblems f o r which is p o t e n t i a l l y d i f f i c u l t f o r a t are s i g n i f i c a n tc o m p u t a t i o n a l andbookkeeping problemsassociatedwiththeanalyticaltreatmentofmixedequilibriumand n o n e q u i l i b r i u ms y s t e m s , and b ) f o r theratecontrollingreactionsarenot f o rt h e s er e a c t i o n sa r eu n a v a i l a b l e . many s y s t e m s o f e n g i n e e r i n g i n t e r e s t , w e l l known and/or r a t e c o n s t a n t s The ACE p a c k a g eo fr o u t i n e se f f e c t i v e l y * F o rc a r b o na b l a t i o ni n a i r , t h ef a i lt e m p e r a t u r ec o n c e p t i s o fa l m o s tn o interest.Graphitesloadedwithmetals, however , t e n d t o f o r m t h i n oxidation resistant oxide coatings which g e n e r a l l y a b l a t e by m e l t i n g . The ACE package w i l l p r o p e r l y i d e n t i f y any suchoxidesformed,andthe f a i l t e m p e r a t u r e mechanism a l l o w s them t o b e a b l a t e d r e a l i s t i c a l l y . F a i l . t e m p e r a t u r e s are u s e f u l f o r m e t a l h e a t s h i e l d s , o f c o u r s e . 51 surmounts the f i r s t of t h e s e p r o b l e m s as d i s c u s s e d i n R e f e r e n c e 4 and as w i l l b e b r i e f l y summarizedsubsequently. However, d i f f i c u l t i e s f a l l i n g i n t o t h es e c o n dc a t e g o r y above f r e q u e n t l yp r e e m p t any s u c h " e x a c t " a n a l y t i c a l t r e a t m e n t .I nt h e s ec a s e s ,v a r i o u sa p p r o x i m a t et r e a t m e n t sc a no f t e ny i e l d u s e f u li n f o r m a t i o n .T h e s ea p p r o x i m a t et r e a t m e n t se s s e n t i a l l yc o n s i s to f notallowingcertainspecies,orgroupsofspecies,toreactwith w i s ee q u i l i b r a t e ds y s t e m . A few examplesof thissort an o t h e r ofapproximate t r e a t m e n t w i l l b ed i s c u s s e di nt h ef o l l o w i n gp a r a g r a p h s .T h e s ed i s c u s s i o n s w i l l be followed bya b r i e f summary o f t h e g e n e r a l k i n e t i c s o p t i o n o f t h e ACE program. 5.3.5.1 Removal o fS e l e c t e dS p e c i e sf r o mt h eS y s t e m the c r u d e s t , way t o a p p r o x i m a t ec e r t a i n i s t o simply removefrom t h e e q u i l i b r i u m s y s t e m t h o s e The s i m p l e s t ,a n dp e r h a p s kinetic effects s p e c i e s whose f o r m a t i o n i s , i n r e a l i t y , s u p p r e s s e d by r e a c t i o n k i n e t i c are removed m e r e l y byremoving e f f e c t s .C o m p u t a t i o n a l l y ,t h e s es p e c i e s t h et h e r m o c h e m i c a ld a t ac a r d sf o rt h a tm o l e c u l ef r o mt h es p e c i e st h e r m o chemicaldatadeck.(Reference 2 d e s c r i b e st h i sd e c ki nd e t a i l . ) 5.3.5.2 Isolated Species C e r t a i nc h e m i c a ls p e c i e s c o n c e n t r a t i o n .T h i st r e a t m e n t may b e f r o z e n i n d i v i d u a l l y a t any given i s s o m e t i m e su s e f u lf o rs p e c i e sw h i c ha r e r e l a t i v e l yn o n r e a c t i v eb e c a u s eo fc h e m i c a lk i n e t i ce f f e c t s .F o re x a m p l e , c o n s i d e rt h ef l o wo f a C-0-H gasover c a r b o ns u r f a c e .I s o l a t i o no f would provide amore e q u i l i b r i u mw e r e a r e l a t i v e l y low t e m p e r a t u r e a b l a t i n g H 2 0 a t i t s b o u n d a r yl a y e re d g ec o n c e n t r a t i o n realisticsurfacethermochemistrytablethaniffull assumed ( s i n c e t h e w a t e r - g a s r e a c t i o n r e l a t i v e l ys l o wa tl o w e rt e m p e r a t u r e s ) .I nt h e o fs p e c i e s a t thesurface is ACE p r o g r a m ,t h i si s o l a t i o n i s achieved by t r e a t i n g t h a t p a r t i c u l a r s p e c i e s a s anatomic a f i c t i t i o u s a t o m i c number s o t h a t i t cannot e l e m e n t , anelementwith r e a c tc h e m i c a l l yw i t h any o t h e r members o ft h ec h e m i c a ls y s t e m s . Com- p u t a t i o n a l l y ,t h i sr e q u i r e ss p e c i f i c a t i o no ft h ei n p u to ft h ea p p r o p r i a t e q u a n t i t y of t h i s" e l e m e n t "( a c c o m p a n i e d amounts o f o t h e r a c t u a l e l e m e n t s ) 52 by a p p r o p r i a t e r e d u c t i o n s i n the and t h e a l t e r a t i o n o f t h e s p e c i e s set for that species t h ef i c t i t i o u s" e l e m e n t . " t h e r m o c h e m i s t r yd a t ac a r d one"atom"of , so t h a . t i t c o n t a l n s o n l y 5.3.5.3IsolatedSubgroupsofElements a subgroup of e l e m e n t s w i t h eachother,butthesubgroup is notinchemical e q u i l i b r i u mw i t ht h es y s t e ma s a whole.The i s o l a t i o n o f a subgroupof e l e m e n t sc a nb ea c h i e v e d by c o n s i d e r i n g t h e s e e l e m e n t s t o b e d i f f e r e n t "eleI n some p r a c t i c a l a p p l i c a t i o n s , t h e r e e x i s t s which a r e i n e q u i l i b r i u m m e n t s "w i t hp r o p e r t i e s the same a s , b u t w i t h a t o m i c numbers d i f f e r e n t f r o m , the c o r r e s p o n d i n ga c t u a le l e m e n t s( e . g . , add 1 0 0 t o t h e a t o m i c numbers o f t h ee l e m e n t si n the i s o l a t e ds u b g r o u p ) .F o rt h i st y p eo fc a l c u l a t i o n ,t h e subgroup t o b e i s o l a t e d i s consideredtobeonecomponent,andthe com- p o s i t i o no f this component is s p e c i f i e d a s t h e a c t u a l r e l a t i v e e l e m e n t a l compositionofthesubgroup to be isolated, but in terms o f t h e a l t e r e d e l e m e n t s .I no r d e rt h a tt h e s ea l t e r e de l e m e n t s o t h e r ,t h es p e c i e st h e r m o c h e m i c a ld a t ac a r d which may b ef o r m e df r o mt h e s ea l t e r e de l e m e n t s , atomic number c o n v e n t i o n . 5.3.5.4 s e t m u s tc o n t a i n a s e t o f s p e c i e s i . e . , specieswiththe same Low T e m p e r a t u r eS u r f a c eE q u i l i b r i u mS u p p r e s s i o n A s p r e v i o u s l yd i s c u s s e d ,e q u i l i b r i u m c i o na t may e q u i l i b r a t ew i t he a c h low s u r f a c et e m p e r a t u r e s( e . g . ,b e l o w problems. t u r er a n g e s i s a ni n c r e a s i n g l yp o o r 2000'R) assump- i n some a b l a t i o n However, i n most s u r f a c et h e r m o c h e m i s t r yt a b l e s ,l o w e rt e m p e r a - are u s u a l l y o n l y o f m i n o r i n t e r e s t ( e . g . , d u r i n g a v e r ys h o r t i n i t i a lp e r i o do f a t r a n s i e n ta b l a t i o np r o b l e m ) .N o n t h e l e s s ,t h es u r f a c e tables must f r e q u e n t l y e x t e n d t o t h e s e l o w e r t e m p e r a t u r e s i n o r d e r t o b e c o m p a t i b l ew i t ht h eh e a tc o n d u c t i o n and a b l a t i o ns o l u t i o n .F o rt h e s e c a s e s , t h e ACE p a c k a g ec o n t a i n sa no p t i o n by w h i c hs u r f a c e( h e t e r o g e n e o u s ) e q u i l i b r i u m w i l l besuppressedandsurfacerecessionnulled (fit = 0 ) below some i n p u t l i m i t t e m p e r a t u r e f o r s o l u t i o n s i n t h e a s s i g n e d t e m p e r a t u r e p a r to ft h es u r f a c et a b l e s . The g a s e sa d j a c e n tt ot h es u r f a c ea r e assumed to beinequilibrium,butheterogeneousequilibrium is not required for t h i ss o l u t i o n .F o r many c a s e s ,s u r f a c et h e r m o c h e m i s t r yt a b l e sg e n e r a t e d inthisfashionareclosertorealitythanfullequilibriumsolutions down t o low t e m p e r a t u r e s . The i n p u tc o n v e n t i o nf o ri m p l e m e n t i n gt h i so p t i o n i s d i s c u s s e di nR e f e r e n c e 2 . 5 . 3 . 5 . 5F u l lT r e a t m e n to fK i n e t i c s The p r e v i o u s d i s c u s s i o n h a s d e s c r i b e d v a r i o u s a p p r o x i m a t e t r e a t m e n t s when a g e n e r a l s o l u t i o n , i n c l u d i n g an " e x a c t " t r e a t m e n t o f r e a c t i o n r a t e e f f e c t s , is e i t h e r i m p o s s i b l e o r i m p r a c t i c a l . The ACE p r o g r a md o e s ,h o w e v e r ,p o s s e s st h ec a p a b i l i t yo f a c c u r a t e l y t r e a t i n g a number o f k i n e t i c a l l y c o n t r o l l e d r e a c t i o n s i n o t h e r w i s e e q u i l i b r a t e ds y s t e m s( e i t h e rc l o s e do ro p e n ) . The r e m a i n d e ro ft h i ss e c t i o n ofkineticseffectswhichareoftenuseful 53 w i l l be devoted t o a brief description of t h i s "exact" t r e a t m e n t of k i n e t i c s . is The i n c l u s i o n o f k i n e t i c a l l y c o n t r o l l e d c h e m i c a l r e a c t i o n s ( 3 9 ) from t h e s e t of equa- a c c o m p l i s h e db yr e m o v i n ge q u i l i b r i u mr e l a t i o n s tionsforcertainspeciesparticipatinginreactionsthat are t o b e are r e p l a c e d by k i n e t i c r a t e t h e A r r h e n i u sf o r m )f o re a c hk i n e t i c a l l yc o n t r o l l e dr e a c t i o n . k i n e t i c a l l yc o n t r o l l e d .T h e s ee q u a t i o n s e q u a t i o n s( o f This i s accomplishedby thereactions first i d e n t i f y i n g t h e p r i m a r y r e a c t i v e s p e c i e s i n andsecond,byallow- which are t o b e k i n e t i c a l l y c o n t r o l l e d , rate ing these species to be created or destroyed only via the kinetic e q u a t i o n s .T h i sa p p r o a c hr e q u i r e s a r e l a b e l i n go fs p e c i e st o be c o n s i d e r e d i nt h ek i n e t i c a l l yc o n t r o l l e dr e a c t i o n s .T h e s es p e c i e sa r ec a l l e dp s e u d o elementssincetheybehavelikeelementsexceptthatthey ordestroyed at rahesspecified by t h e r e a c t i o n may be c r e a t e d rate equations. C o m p u t a t i o n a l l y ,t h ei n c l u s i o no fk i n e t i c sr e s u l t si nt h ea d d i t i o n term t o t h e e l e m e n t a l b a l a n c e e q u a t i o n s unknowns t c t h es y s t e m e q u a l i n number t o t h e number o f s p e c i e s whose c o n c e n t r a t i o n s a r e k i n e t i i . e . , thepseudo-elements. The r e l a t i v ec r e a t i o n and c a l l yc o n t r o l l e d , d e s t r u c t i o n r a t e s of a l l pseudo-elements i n a g i v e n r e a c t i o n are r e l a t e d by stoichiometry,however, s o t h e number o f a d d i t i o n a l unknowns remaining is e q u a lt ot h e number o f k i n e t i c a l l y c o n t r o l l e d r e a c t i o n s . The r e a c t i o n r a t e s , fromwhich t h e p s e u d o - e l e m e n t c r e a t i o n o r d e s t r u c t i o n rates d e r i v e , are g i v e n by of a r a t e - o f - c r e a t i o n o r d e s t r u c t i o n f o rt h e s ep s e u d o - e l e m e n t s .T h i sa d d sa d d i t i o n a l f o re a c hk i n e t i c a l l yc o n t r o l l e dr e a c t i o n m o ft h ef o r m where t h e sums a n d p r o d u c t s a r e o v e r t h e s p e c i e s a n d p s e u d o - e l e m e n t s p are t h e s t o i c h i o m e t r i c c o e f f i c i e n t s , c o n s t a n tf o rr e a c t i o n 54 ,K i s t h ee q u i l i b r i u m Pm ( 7 3 ) I and kFm i s t h ef o r w a r d r a t e c o n s t a n t . I n t h e r e a c t a n t s a n dp r o d u c t sr e s p e c t i v e l y I. n N and s u p e r s c r i p t s R and P d e n o t e (72) 1' present formulation, kFm i s r e p r e s e n t e d b y an A r r h e n i u s t y p e f u n c t i o n Cpm i s t h et e m p e r a t u r ee x p o n e n t ,a n d B, i s a c o n s t a n tr e p r e s e n t i n g a m u l t i t u d e of phenomena. Em, C p , and Bm must be input to the program for each kinetically controlled reaction. T h e s ec o n s t a n t ss h o u l db eb a s e do ne x p e r i m e n t a ld a t a . The u n c e r t a i n t y i n , orunavailabilityof,thesedatafor many c h e m i c a l s y s t e m s o f i n t e r e s t i n ablationproblemsfrequentlyrepresent a significant constraint on the a p p l i c a t i o n o f the k i n e t i c s o p t i o n of ACE program t o t h e s e p r o b l e m s . More detaileddiscussionofthecomputationaltreatmentofkinetics is presented i nR e f e r e n c e s 4 and 11. R e f e r e n c e 2 d e s c r i b e st h ei n p u td e t a i l sf o rt h e kineticallycontrolledreactiondata. where Em is t h e a c t i v a t i o n e n e r g y , 5.4 SOME CODING DETAILS are employed i n t h e ACE program t o s o l v e The c o m p u t a t i o n a l p r o c e d u r e s t h ee q u a t i o n s s e t f o r t hi nS e c t i o n5 . 3 ,b r i e f l yd i s c u s s e dh e r e .C o n s i d e r a b l y g r e a t e r d e t a i l i s presentedinReference 4 and i n p a r t i c u l a r , T a b l e I of that reference. 5 . 4 . 1B a s i cS o l u t i o nT e c h n i q u e The b a s i c s o l u t i o n t e c h n i q u e may b e i l l u s t r a t e d e x a m p l e ,a no p e ns y s t e mw i t hu n e q u a ld i f f u s i o nc o e f f i c i e n t s , p h a s em a t e r i a lr e m o v a l ,a n d 5.3.3) . by c o n s i d e r i n g , f o r no condensed no k i n e t i c s ( i . e . , a s d i s c u s s e d i n S e c t i o n F o rt h i ss y s t e m ,t h eb a s i ce q u a t i o n sd e f i n i n gt h e. p r o b l e ma r et h e e l e m e n t a lc o n s e r v a t i o ne q u a t i o n s ( 6 7 ) , t h et o t a lp r e s s u r ee q u a t i o n( 5 6 ) , thereactionequilibriumequations(55), r e l a t i o n( 5 3 ) . and o n e h e t e r o g e n e o u s v a p o r p r e s s u r e shows t h a t t h e r e are as many so c l o s u r e is obtained. The t a b l eo fS e c t i o n5 . 3 . 3 knowns a s unknowns i n t h e s e e q u a t i o n s Summarizing these equations , kPi T P = 0 55 K k= 1 K The number o f unknowns could immediately be reduced by I - K t h r o u g h thedirectsubstitutionof(77), as s o l v e d f o r P i , i n t o t h e o t h e r r e l a t i o n s . I t provesadvantageous,however, t o c o n t i n u e t o t r e a t t h e f u l l s e t of equations,andtosubsequentlyutilizethissubstitutionduringtheitera- t i v e convergenceprocedure. a l g e b r a i ce q u a t i o n s The s o l u t i o no ft h e s es i m u l t a n e o u sn o n l i n e a r i s basedon Newton-Raphson i t e r a t i o n . S i n c e t h i s p r o a more l i n e a r form e t c . ) it i s w e l l t o examine t h e s e t c e d u r e i s a c c e l e r a t e d by c a s t i n g t h e e q u a t i o n s i n t o ( v i at r a n s f o r m a t i o n s ,s u b s t i t u t i o n s , of e q u a t i o n sa b o v e . With t h eb o u n d a r yl a y e re d g e ,c h a ra n dp y r o l y s i sg a s compositiongivenas w e l l a st h e B', ( 7 5 )a n d( 7 6 )a r el i n e a rr e l a t i o n s F i s r e a s o n a b l yc o n s t a n t I. n con9 (77)and(78) are l i n e a r r e l a t i o n s b e t w e e n t h e i n P i , En Pk and thelattervariablebeingapproximatelylinearin 1/T. b e t w e e nt h e trast, Rn Kpi, Pi a n dp r o v i d i n gt h a t The ACE p r o g r a m t a k e s a d v a n t a g e o f t h i s s i t u a t i o n s p e c i e s which a r e s i g n i f i c a n t i n t h e ( 7 6 ) i n terms o f Pi by t r e a t i n g t h o s e mass a n d p r e s s u r eb a l a n c e s( 7 5 )a n d less s i g n i f i c a n t s p e c i e s i n andthe terms of t h e i r i n Pi. The Newton-Raphson p r o c e d u r e , a s a p p l i e d by t h e ACE program,canbe summarized by c o n s i d e r i n g a s e t o f e q u a t i o n s o f t h e g e n e r a l f o r m f . (x1,x2,. . ., x i , . . . ) 3 A t any p o i n t i n t h e s o l u t i o n p r o c e d u r e t h e r e e x i s t s x?, for all the variables r e l a t i o n s and w i l l l e a dt o a s e t of e s t i m a t e s , which w i l l i n g e n e r a l n o t s a t i s f y a non-zerovalue Newton-Raphson methodproceeds evaluatingthechangeineach = 0 a l l ofthe of the f j , namely, t o" d r i v e "t h e s ee r r o r st o w a r dz e r o E j. The by unknown v a r i a b l e , Axi, whichwouldreduce a l lt h ee r r o r st oz e r oi ft h ef u n c t i o n s , f , were l i n e a r . The l i n e a r j approximation i s based on t h e c u r r e n t v a l u e s o f t h e unknown v a r i a b l e s and t h e c o r r e s p o n d i n g a r r a y o f v a l u e s o f t h e p a r t i a l d e r i v a t i v e s a f . / a x i 56 I Thus which is l o c a l l y c o r r e c t a n d i n the l i n e a r approximation. is i n t e g r a t e d t o The s o l u t i o n of (80) is of p a r t i a l d e r i v a t i v e s a p p e a r i n g i n (81) i s s i m p l yt h em a t r i x ( 8 0 ) . In t h e ACE program t h e f o r m u l a t i o n o f t h e p a r t i a ld e r i v a t i v e s u s e s the v a r i a b l e s , Iln Pi Iln T and I n and ( 8 1 ) y i e l d s ,f o re x a m p l e , w h e r et h ea r r a y i n v e r s e of t h e a r r a y i n m which i f t a k e n a s l i n e a r sincethe error. a l l t h e way t o s o l u t i o n y i e l d s desired change i n t h e f u n c t i o n s +n e q u a l ? ye x a c tr e l a t i o no b t a i n e d is s i m p l y t h e n e g a t i v e o f t h e from (82) i s tThechoiceof En P . p e r m i t s a m a t r i x r e d u c t i o n by t h e u s e o f t h e s i m p l e algebraic substituhon, previously mentioned after (78) , p r i o r t o m a t r i x inversion. 57 which i f t a k e n as l i n e a r a l l t h e way t o s o l u t i o n y i e l d s The ACE programuses (85) f o r a l l s p e c i e s s i g n i f i c a n t i n mass balancesand ( 8 3 )f o rt h eo t h e r s . 5 . 4 .R 2 estriction on C o r r e c t i o n s The s e t of c o r r e c t i o n (Axi) * canbethoughtof as a v e c t o r i n t h e s p a c eo ft h ei n d e p e n d e n tv a r i a b l e sw h i c h i s added t o t h e c u r r e n t v e c t o r E x p e r i e n c eh a s shown t h a t approximation x* t o y i e l d a new e s t i m a t e x:* 1 it i s f r e q u e n t l y u n w i s e t o p r o c e e d a l o n g t h i s c o r r e c t i o n v e c t o r t h e f u l l . amount i n d i c a t e d by ( 8 3 )o r (85) R a t h e r , it i s b e t t e r t o p r o c e e d same d i r e c t i o n . way, a l t h o u g hp r e s e r v i n gt h e todepartfromthisvector, . A t o t h e r times, i t i s e x p e d i e n t and seek a n o t h e r b a s e d some v a r i a b l e a n d e l i m i n a t i n g a limited on f r e e z i n g t h e v a l u e o f a correspondingequation. The s c a l i n g o f t h e c o r r e c t i o n v e c t o r i s s u c ha st o l i m i t changes i n the partial pressures of major species to increases of one order of magni- t u d ea n dd e c r e a s e so ft h r e eo r d e r so fm a g n i t u d e ,a n dc h a n g e so ft e m p e r a t u r e t oa p p r o x i m a t e l y 20 p e r c e n t . M o l e c u l a rw e i g h t ,t e m p e r a t u r e correctionsarefrozen and condensedphaseconcentration and a new c o r r e c t i o n v e c t o r g e n e r a t e d i f t h e i n i t i a l set ofcorrectionsindicateexcessivetemperatureormolecularweight e x c u r s i o n s , a c o n t r a d i c t o r yt e m p e r a t u r ec h a n g e ,o rn e g a t i v ec o r r e c t i o n so n newly i n t r o d u c e d c o n d e n s e d s p e c i e s . The f o r m u l a t i o n of t h e s e and o t h e r s c a l i n g a n d f r e e z i n g i s an e s s e n t i a l f e a t u r e o f t h e criteria Because o f t h e s e f e a t u r e s , w e l l f o r m u l a t e d ,p h y s i c a l l yu n i q u e ACE program. convergence i s v i r t u a l l y a s s u r e d f o r problems. 5.4.3 Base S p e c i e s The d i s c u s s i o n o f S e c t i o n 5.3 d e s c r i b e d t h e e q u i l i b r i u m r e a c t i o n e q u a t i o n sa se q u a t i o n sg i v i n gt h ef o r m a t i o n Thus t h e r e a c t a n t s a r e e l e m e n t s of a s p e c i e sf r o ma t o m i ce l e m e n t s . and t h e p r o d u c t s a r e u s u a l l y m o l e c u l e s . T h i s scheme h a s t h e a d v a n t a g e o f f o r m a l s i m p l i c i t y , s i n c e t h e s t o i c h i o m e t r i c coefficientsneededintheequilibriumequationsaregivendirectly 58 by t h e productspecieschemicalformula.This scheme c a nh a v ec o m p u t a t i o n a ld i s - a d v a n t a g e s ,h o w e v e r ,s i n c et h ea t o m i ce l e m e n t s are f r e q u e n t l y n o t p r e s e n t t o any g r e a t e x t e n t i n t h e e q u i l i b r i u m s y s t e m . T h i s i n i t s e l f results i n no d i s a d v a n t a g e . I f one mass b a l a n c e ( e . g . d e f e a tc o n v e r g e n c e . , however , a molecule (e. g. , CO) dominates more t h a n , C and 0) , l o s s of s i g n i f i c a n t f i g u r e s c a n slow o r I t i s more d e s i r a b l e t o w r i t e t h e e q u i l i b r i u m r d a c t i o n s ( 3 8 ) as w e l l as t h e m a s s b a l a n c e s i n (37) and terms o fr e a c t a n ts p e c i e sw h i c h are a c t u a l l yp r e s e n ti na p p r e c i a b l ea m o u n t s .T h e s es p e c i e sa r et e r m e d" b a s e s p e c i e s "( f r o mR e f e r e n c e 1 2 ) s i n c ea l lo t h e rs p e c i e sa r et a k e nt ob e formed from them. The ACE program s e l e c t s t h e b a s e s p e c i e s f r o m t h e c a n d i d a t e s p e c i e s t h e r m o c h e m i c a ld a t aw h i c ha r ei n p u ta ss p e c i f i e di nR e f e r e n c e 2. programautomaticallyselectsasbasespeciesthefirstset s a t i s f y i n gt h er e q u i r e m e n tt h a t thisbasespeciesset may beformedfrom and ( 2 ) t h a t no b a l a n c e dr e a c t i o nc a nb e i n v o l v i n go n l yb a s es p e c i e s . eachelement. (1) a l l o t h e r s p e c i e s One b a s es p e c i e s The of s p e c i e s written may b ec o n s i d e r e dt or e p r e s e n t T h u s , t h eb a s es p e c i e sa r ee s t a b l i s h e d t h e u s e r i n p u t st h ec a n d i d a t es p e c i e st h e r m o c h e m i c a ld a t a . by t h e o r d e r i n which The c a l c u l a t i o n o ft h es t o i c h i o m e t r i cc o e f f i c i e n t s and e q u i l i b r i u m c o n s t a n t s a p p r o p r i a t e t o any s e t o f b a s e s p e c i e s is h a n d l e d a u t o m a t i c a l l y by theprogram. While t h e s e l e c t i o n computationaleconomies whichoccur of b a s e s p e c i e s may b e r e a l i z e d is n o tu s u a l l yc r i t i c a l ,c e r t a i n i f t h e basespeciesarethosespecies i n r e l a t i v e l yh i g hc o n c e n t r a t i o n si nt h es y s t e m . Also t h e t r e a t - ment of a c h e m i c a lk i n e t i c s( S e c t i o n 5 . 3 . 5 ) may havean i n f l u e n c e on t h e s e l e c t i o n of t h eb a s es p e c i e s .T h e r e f o r et h eU s e r s G ' u i d e( R e f e r e n c e 2) recommends c e r t a i n i n p u t p r o c e d u r e s f o r f o r c i n g t h e s e l e c t i o n o f d e s i r a b l e basespecies. 59 SEmION 6 COUPLING THE TRANSFER COEFFICIENT CALCULATION TO THE SURFACE THERMOCHEMISTRY SOLUTION (111 - I V COUPLING) AND TO THE IN-DEPTH TRANSIENT SOLUTION (111 V COUPLING) - 5 S e c t i o n 5 above makes it p l a i n t h a t t h e t r a n s f e r c o e f f i c i e n t a p p r o a c h t o modelingdiffusionalfluxes is builtintothesurfacethermochemistry s o l u t i o np r o c e s s ,s i n c et h ec h e m i s t r yr o u t i n e su s e i n j e c t i o n and t h e B; c o n c e p t d e p e n d s n a t u r a l l y o n idea. An a d d i t i o n a l , B ' as t h em e a s u r eo f mass C the t r a n s f e r c o e f f i c i e n t more e x p l i c i t , l i n k t o t h e t r a n s f e r c o e f f i c i e n t model i s peueCM. The p r e - e x p o n e n t i a lf a c t o r s Bm t h es c a l i n go fk i n e t i ce f f e c t so n are i n p u t d i r e c t l y t o t h e c o u p l e d c o d e , a n d t h e s e are normalizedon peueCM ACE package, as n o t e di nS e c t i o n 5 a b o v e .T h i sr e q u i r e st h a t a i n s i d et h e peueCM v a l u e b e c o m m u n i c a t e d d i r e c t l y t o t h e f o ru s ei nn o r m a l i z i n ga n d B 's. ACE packagewithevery call O t h e rc o u p l i n gl i n k sa p p e a ri nt h es u r f a c e m t e m p e r a t u r ed e p e n d e n c ea n dt h e" b l o w i n gr e d u c t i o n "o ft h et r a n s f e rc o e f f i c i e n t : t h e s e are d i s c u s s e d i n S e c t i o n The I11 - V couplingbetweenthetransfercoefficientand d e p t hs o l u t i o no c c u r s r o f t h es u r f a c ee n e r g ye v e n t s . sentedinSection 60 9 below. the i n course, i n t h e t r a n s f e r c o e f f i c i e n t e q u a t i o n f o r The r e l e v a n te q u a t i o n sh a v ea l r e a d yb e e np r e - 4 . 3 andrequirenofurther comment. SECTION 7 SOLUTION COUPLING THE FLOW-FIELD (OR BOUNDARY LAYER EDGE STATE) STATE SOLUTION (I1 - I V COUPLING) (PHASE 11) TO THE SURFACE COEFFICIENT SOLUTION (I1 - I11 COUPLING) AND TO THE TRANSFER B e f o r ed i s c u s s i n gt h eP h a s e I1 and I11 c a l c u l a t i o n s t h e m s e l v e s , w i l l be of interest to examine t h e c o u p l i n g l i n k s b e t w e e n t h e P h a s e I1 solution,whichreferstothespecificationoftheboundarylayerouter edge s t a t e a r o u n dt h e body ( p r e s s u r e , s t a t i c e n t h a l p y , r e c o v e r y e n t h a l p y , and t r a n s p o r t p r o p e r t i e s ) , and t h e P h a s e I11 andPhase I V solutions. it The b o u n d a r yl a y e re d g es t a t ei n f l u e n c e st h es u r f a c et h e r m o c h e m i s t r y s t a t es o l u t i o n ,a sh a sa l r e a d yb e e nd i s c u s s e di nS e c t i o n thedirect communication of l o c a l p r e s s u r e package. The edge s t a t e i n t u r n as w i l l bediscussedinSection 5 above,through and r e c o v e r y e n t h a l p y t o t h e i s i n f l u e n c e db yt h ec u r r e n t ACE body s h a p e , 8 below. The b o u n d a r y l a y e r e d g e s t a t e i n f l u e n c e s t h e t r a n s f e r c o e f f i c i e n t calculationthroughthedistributionsofpressure body,asmightbeexpected.This a n de n t h a l p ya r o u n dt h e i s d i s c u s s e di nS e c t i o n 9 below. it i s a l s o u s e f u l A l l ofthesecouplinglinksareofinterest,but a tt i m e st ob ea b l et os u p p r e s s some o r a l l o f them.Such c o u p l i n gs u p p r e s sion is particularlydesirableinparametricstudiesattemptingtostudy , one e f f e c t a t a t i m e , r a t h e r t h a n c o u p l e d e f f e c t s it is d e s i r e d t o m a t c he x p e r i m e n t a ld a t ai n s u r f a c et e m p e r a t u r e s ) .T h e r e f o r e ,t h e a l l o wv a r y i n gd e g r e e so fd e c o u p l i n g . foreachsurfacepoint a. and i n s t u d i e s i n which some r e s p e c t( s u c ha sm e a s u r e d MACABRE c o d eh a sb e e nw r i t t e nt o A l l o ft h ef o l l o w i n go p t i o n se x i s t on a body a t anytime: S u r f a c et e m p e r a t u r ea n dr e c e s s i o n r a t e d i s c o v e r e d by e n e r g yb a l a n c e ,o r b .S u r f a c et e m p e r a t u r ea n dr e c e s s i o n f u n c t i o n s of time by u s e r( O p t i o n r a t e a s s i g n e d as 2) Under ( a ) w e have a . 1 .F u l ls u r f a c ee n e r g yb a l a n c ew i t hc o n v e c t i o na n dt h e r m o c h e m i s t r(yO p t i o n a.2. 1) No c o n v e c t i o n ,n or e c e s s i o n ,a s s i g n e dh e a tf l u x( O p t i o n 2) 61 Under (a.1.) w e have L o c a lp r e s s u r ea n d / o rl o c a lr e c o v e r ye n t h a l p ya s s i g n e db yu s e r a.l.1 as a f u n c t i o n of t i m e a . 1 . 2L o c a lp r e s s u r ea n d / o rl o c a lr e c o v e r ye n t h a l p y computedby a p p r o p r i a t es u b r o u t i n e sa u t o m a t i c a l l y F o rb o t h( a . l . l ) a n d 1 a.1. (?) .1 (a.1.2) w e have Localheattransfercoefficientassigned functionof a.1. 1 (7).2 by u s e r as a time Localheattransfercoefficient computedautomatically b ya p p r o p r i a t es u b r o u t i n e s Forboth (a.1. (3)-1) a n d( a . 1 . ($) . 2 ) correctiontothetransfercoefficient. 62 w e haveanydegreeof"blowingreduction" SECTION 8 THE BOUNDARY LAYER EDGE STATE SOLUTION state i s of t h e body s h a p e , shockshape , t h es h o c kl a y e rf l o w ,a n dt h eb o u n d a r yl a y e rf l o w .F o r many e n g i n e e r i n gp u r p o s e s ,h o w e v e r , a detailedanalysisofthissort i s n o tp o s s i b l e , and c e r t a i n s i m p l i f y i n g a s s u m p t i o n s are n e c e s s a r y which b r i n g t h e a n a l y s i s t o a more t r a c t a b l e l e v e l of s o p h i s t i c a t i o n .S i n c et h ea p p l i c a t i o n fortheanalyticaltechniquebeingdiscussedhere i s primarilyre-entry v e h i c l e n o s e t i p s and h e a t s h i e l d s , t h e f i r s t s i m p l i f y i n g a s s u m p t i o n is t h a t theentireregionofinterest i s immersed i n a boundarylayeredgeflowwhich h a sp a s s e dt h r o u g ht h ed e t a c h e d bow shock a t o r n e a r t h e s t a g n a t i o n p o i n t . Thus, it i s n o t n e c e s s a r y t o c a l c u l a t e a shockshapeandshocklayerflow s i n c et h es h o c kl a y e rg a so fi n t e r e s t i s assumed t o p a s s t h r o u g h a normal shock. The n e x ts i m p l i f y i n ga s s u m p t i o n i s t h a tt h eb o u n d a r yl a y e re d g eg a s s t a t e i s determined by a n i s e n t r o p i c e x p a n s i o n f r o m t h e body s t a g n a t i o n condition. With t h ee d g eg a se n t r o p yf i x e da tt h es t a g n a t i o nv a l u e ,t h e completeedgegasthermodynamic s t a t e is determined by t h e s p e c i f i c a t i o n of one more p r o p e r t y ,t y p i c a l l yt h el o c a ls t a t i cp r e s s u r e .T h u s ,t h ed e t e r m i n a t i o no ft h ee d g eg a s s t a t e r e d u c e s t o a problemof firstcalculatingthe s t a g n a t i o ne n t r o p y ,w h i c h i s r e l a t i v e l ys t r a i g h t f o r w a r d ,t h e nd e t e r m i n i n g t h ep r e s s u r ed i s t r i b u t i o na r o u n dt h ea r b i t r a r y body s h a p e , and f i n a l l y evaluatingtheotheredge s t a t e p r o p e r t i e s of i n t e r e s t by c a r r y i n g o u t t h e i s e n t r o p i ce x p a n s i o nt ot h e known l o c a lp r e s s u r e . The t e c h n i q u e su s e dt o carryoutthesesteps are d e s c r i b e d i n t h e s u b s e c t i o n s b e l o w . A c o m p l e t ea n da c c u r a t ea n a l y s i so ft h eb o u n d a r yl a y e re d g e a v e r y complexproblem 8.1 requiringsimultaneousconsideration THERMOCHEMICAL COMPUTATIONS Thermochemicalcomputations are r e q u i r e d i n a number o f i n s t a n c e s i n state calculations. As mentioned i n S e c t i o n 5 , t h e s e c a l c u l a t i o n s a r e c a r r i e d o u t by t h e ACE r o u t i n e . i t hasbeenfoundmostconvenient to store Forboundarylayergasproperties, t h ec o u p l e dc o d e ,i n c l u d i n gb o u n d a r yl a y e re d g e ACE s o l u t i o n s f o r t h e r e q u i s i t e t a b u l a r formandsupply thermodynamicand transportpropertiesin thistableofpropertiestothe programasinput. The o p t i o n a l s o e x i s t s of m e r e l y s p e c i f y i n g t h e d e s i r e d v a l u e s o f t h e i n d e pendent variables in this'Wollier Chad' and c a l l i n g f o r ACE c a l c u l a t i o n s t o d e t e r m i n et h er e q u i r e dp r o p e r t i e s . The p r o p e r t i e s a t t h eg a s s t a t e o fi n t e r e s t d u r i n gt h ec o u r s eo f a solutionarefound by i n t e r p o l a t i o n b e t w e e n t a b u l a r 63 Chart approachhas a number o fa d v a n t a g e s similar t ot h ep r e c a l c u l a t e ds u r f a c et a b l e s :r e u s e a b i l i t y ,e x t e r n a le v a l u a tion of the t a b l e a c c u r a c ya n ds u f f i c i e n c y ,a n de l i m i n a t i o no fn o n - c o n v e r g e d solutions. However, c a r e m u s tb et a k e nt og u a r a n t e et h a ti n t e r p o l a t i o nw i t h i n the chart w i l l y i e l d s u f f i c i e n t l y a c c u r a t e state properties. values. The p r e c a l c u l a t e d M o l l i e r The M o l l i e r C h a r t u s e s p r e s s u r e a n d e n t h a l p y a s i n d e p e n d e n t v a r i a b l e s . Dependentvariablesfound a t e a c h( P - h )c o m b i n a t i o na r et e m p e r a t u r e ,d e n s i t y , P r a n d t l number, v i s c o s i t y ,e n t r o p y , and i s e n t r o p i ce x p o n e n t .L i n e a ri n t e r p o - l a t i o n is usedeverywherewithinthechart,howeverpressuresanddensities i s carried limit t h e number o fp r e s s u r e - e n t h a l p y areconvertedtologarithmicformbeforethislinearinterpolation out.Sincethecomputercodedimensions i t i s recommended t h a t t h e u s e r combinationswhichcanbeinput, check t h e a c c u r a c y of i n t e r p o l a t e d s o l u t i o n s b e f o r e g o i n g a h e a d w i t h p r o b l e m s o l v i n g . Thiscanbeaccomplished by p l o t t i n g t h e p r o p e r t i e s o f i n t e r e s t , d e t e r m i n i n g theregionofmostlikelyerrorwithinthelinearinterpolationapproximation, andcomparingan interpolatedpointwithan"exact"solution from t h e ACE program. i s t y p i c a l l y s t a r t e d by u s i n g t h e The edge s t a t e s o l u t i o n p r o c e d u r e specifiedvalues of s t a g n a t i o n p r e s s u r e and s t a g n a t i o n e n t h a l p y t o d e t e r m i n e t h ee n t r o p yl e v e lf r o mt h eM o l l i e rC h a r t . bodymust 8.2 The s t a t i c p r e s s u r e a r o u n d t h e thenbefound,asdiscussedinthenextsection. STATIC PFiESSUFiE DISTRIBUTION The t e c h n i q u e u s e d t o c a l c u l a t e t h e s t a t i c p r e s s u r e d i s t r i b u t i o n i n c l u d e s a v a r i a t i o n of Newtoniantheory inthesubsonicflownearthestag- nationpoint,matchedto a P r a n d t l Meyer e x p a n s i o n d o w n s t r e a m o f t h e s o n i c point. i nt h e s er e g i o n sa r ee x p l a i n e db e l o w . The methodsused 8.2.1 Subsonic Region The p r e s s u r e d i s t r i b u t i o n i n t h e s u b s o n i c r e g i o n is found by a v a r i a - 1 3 , i n which t h e Newtonian t i o no ft h et e c h n i q u es u g g e s t e di nR e f e r e n c e is modified t o i n c l u d e blunt bodies w i t h s o n i c The b a s i c Newtonian p r e s s u r ee x p r e s s i o n i s pressure distribution relation corners. P = - Pm + (1 - Fa, cos2 B wherebarredquantitiesarenormalizedbythestagnationpressure,and t h ea n g l eb e t w e e nt h e atthepointofinterest. clature. 64 body a x i s ofsymmetryand The sketchbelow B is a normal t o t h e body s u r f a c e w i l l helptoclarifythis nomen- external flow R* = r a d i u s t o s o n i c p o i n t along a s u r f a c e n o r m a l F i g u r e 1 7 . Diagramof N o s e t i p Showing Nomenclature are smoothlycontoured i s subsoniceverywhere ahead of aconvex s h a r p c o r n e r and t h e f l o w e x p a n d s t o t h e s u p e r s o n i c cond i t i o n a t t h ec o r n e r .F o rt h e s es h a p e s ,t h eN e w t o n i a np r e s s u r ed i s t r i b u t i o n is g r o s s l yi n a c c u r a t es i n c e it d e p e n d so n l yo nl o c a ls l o p e . I t is obvious t h a t h ep r e s s u r es h o u l dv a r ys m o o t h l yf r o mt h es t a g n a t i o nv a l u e = 1.0 atthestagnationpointtothesonicpointpressure P = P* a t t h e s h a r p c o r n e r . However, f o r a sphere-coneshapeas shown i nF i g u r e1 7a b o v e ,t h e Newtonian d e s c r i p t i o no fp r e s s u r e would c a l l f o r a c o n s t a n td e t e r m i n e d by t h es u r f a c es l o p ea l o n gt h ee n t i r ec o n es u r f a c e . I t was found i nR e f e r e n c e 13 t h a t a modification of the Newtonian distribution involving an empirical A class of nosetips of particular interest bluntbodieshavinglocalslopessuchthattheflow 65 expressionforpressureon same d i a m e t e r ( s o n i c p o i n t a f l a t d i s c o ft h e d i a m e t e r )w o u l dg r e a t l yi m p r o v ep r e s s u r ep r e d i c t i o n s . The s u g g e s t e dp r e s s u r e is d i s t r i b u t i o ne x p r e s s i o n ,m o d i f i e df o rn o n z e r oa m b i e n tp r e s s u r e , 1 where PFD i s t h e f l a t d i s c p r e s s u r e d i s t r i b u t i o n e x p r e s s i o n , r u n n i n ql e n g t h , n o s er a d i u so r % R*, i s t h en o s er a d i u s , a sd e f i n e di nF i g u r e and Rma, 17. expressionsforthesonicpointpressure s i s t h es u r f a c e i s the maximum o f e i t h e r t h e I t i s n e c e s s a r yt os p e c i f y F* and t h e f l a t d i s c p r e s s u r e FFD I t i s i nt h i sl a s te x p r e s s i o nt h a tt h e i no r d e rt ou s et h i se x p r e s s i o n . p r e s e n ta n a l y s i sd i f f e r sf r o mt h a to fR e f e r e n c e 13. The s o n i c p o i n t p r e s s u r e h a s b e e n c h o s e n t o b e d e s c r i b e d by t h e i d e a l g a se x p r e s s i o n The i s e n t r o p i ce x p o n e n t y i s evaluatedatthestagnationpoint t o r e m a i nc o n s t a n to v e rt h er e g i o no fi n t e r e s t . and i s assumed The assumption of a c o n s t a n t y i n t r o d u c e s a n e g l i g i b l ys m a l le r r o ri n t ot h ea n a l y s i s .A n o t h e rq u a n t i t y of i n t e r e s t i s t h el o c a t i o no ft h es o n i cp o i n t , s*. F o rt h es h a r pc o r n e r e d b o d i e sd i s c u s s e da b o v e , s* i s m e r e l yt h ed i s t a n c et ot h ec o r n e r . However, a s these s h a r p c o r n e r e d b o d i e s a b l a t e o r f o r o t h e r s m o o t h l y - c o n t o u r e d b o d i e s , m u s tb el o c a t e di no r d e rt ou s eE q u a t i o n( 8 7 ) . The Newtonian p r e s s u r e d i s t r i b u t i o n i s employedhere t o g i v e t h e s u r f a c e a n g l e a t t h e s o n i c t h es o n i cp o i n t p o i n t of B * = cos -l F 1-Pm Thus t h e s u r f a c e r u n n i n g l e n g t h 66 s* canbefound by u s i n g t h e (89) known s u r f a c e slopesandrunninglengthsalongthesurface and i n t e r p o l a t i n g b e t w e e n The f l a t d i s c p r e s s u r e d i s t r i b u t i o n u s e d i n R e f e r e n c e x = them. 13 i s s -4 (90) I t w i l l be shown t h a t t h e h e a t t r a n s f e r c o e f f i c i e n t a t t h e s t a g n a t i o n p o i n t depends upon t h e s t a g n a t i o n p o i n t v e l o c i t y g r a d i e n t (see E q u a t i o n( 1 2 0 ) ) . Thisgradientinturndepends upon t h e s e c o n d d e r i v a t i v e o f p r e s s u r e w i t h d i s t a n c ea c c o r d i n gt ot h ef o l l o w i n gr e l a t i o n O'S ( 8 7 ) twice y i e l d s Differentiatinq Equation d2P ds which f o r RN - + -2 2 RNRmax . appioachinginfinityyields n a t i o np o i n ti fE q u a t i o n d2FFD 2 ds a zerovelocitygradientatthestag- ( 9 0 ) i s employed a s i t s t a n d s . T o c o r r e c t t h i s un14 were e v a l u a t e d t o o b t a i n a flat disc stagnation point velocity gradient expression realisticresult,thedataofReference (93) By combiningEquations (911, ( 9 2 ) , and(93) , t h ef o l l o w i n gr e s u l t is obtained 67 This expression provides a more correct stagnation point pressure distribution ( g o ) , has been modified to blend The flat disc pressure distribution, Equation smoothly with the new second derivative requirements, and the definition A has of been x 5 = changed G slightly: p Equations ( 9 5 ) and ( 9 6 ) are the final forms of the f l a t d A > cpressure distribution relations which are used in the present analysis. 8.2.2 Downstream of the Sonic, Point ( 8 7 ) above, the local static Using the expression given in Equation pressure is exactly equal to the Newtonian pressure for all body shapes at the sonic point. The Newtonian description of pressure is then used downstream of the sonic point until the Prandtl-Meyer match point is reached, at which point a Prandtl-Meyer expansion expression is used for the remainder of the running length. The match point is determined by the requirement of continuity of pressure and pressure gradient along the surface: ( 9 7) Prandtl-Meyer It is more pressure : convenient work in terms to of the angle B and dimensionless Prandtl-Meyer The 68 Newtonian pressure expression be candifferentiated to give For a perfectgas,thePrandtl-Meyerpressuregradient - - - is yFM 2 (100) Prandtl-Meyer while the relation between pressure Equatingthe and Mach number i s two g r a d i e n t e x p r e s s i o n s y i e l d s t h e f o l l o w i n g e q u a t i o n f o r t h e match point: T h i s i s an i m p l i c i t e q u a t i o n f o r s o l u t i o n i s foundby Pressuresonthe carryingout B M , s i n c e F and M are f u n c t i o n so f a Newton-Raphson B. The iteration procedure. bodydownstream of the match p o i n t a r e f o u n d by a Prandtl-Meyerexpansionfromthematchpointconditionsto t h el o c a ls u r f a c es l o p e . The g o v e r n i n ge q u a t i o n si nt h i sr e g i o na r e E q u a t i o n (101 ) above r e l a t i n g p r e s s u r e t o Mach number,and which r e l a t e s Mach number t o l o c a l s u r f a c e a n g l e . Once a g a i n ,t h ee q u a t i o n s are i m p l i c i t s i n c e M i s a f u n c t i o n o f B , t h e r e f o r e a Newton-Raphson i t e r a t i o n procedure w a s d e v i s e d t o s o l v e t h e s e e q u a t i o n s . 69 Examples oftheAccuracy 8.2.3 of t h e Method The a c c u r a c y o f t h e p r e s s u r e p r e d i c t i o n t e c h n i q u e d e s c r i b e d a b o v e i s demonstratedwithfoursampleproblems. h a l fa n g l es p h e r ec o n ew i t h a s h a r pc o r n e r The f i r s t o f t h e s e as shown i n F i g u r e i s a 60" 1 8 . With t h i s is fixedatthesharpcorner,thereand Rmax a r e known from t h e geometry. The p r e s s u r e p r e d i c t i o n s h a p e ,t h el o c a t i o no ft h es o n i cp o i n t fore B* , s* , for this shape is compared i n F i g u r e 18 t o t.be d a t a o f R e f e r e n c e 15 , and a l s o t c t h e N e w t o n i a np r e s s u r ed i s t r i b u t i o nf o rt h e same s h a p e . Agreement is s e e n t o b e e x c e l l e n t o v e r m o s t o f t h e s u r f a c e , w h i c h is tobeexpected sincetheempiricalportionsofthetheory were i n v e n t e d t o s o l v e t h i s v e r y problem. The s h o r t c o m i n g so fN e w t o n i a nt h e o r ya r ea l s or e a d i l ya p p a r e n ti n this figure . Figure 1 9 d e m o n s t r a t e st h eu t i l i t yo ft h i sp r e s s u r ep r e d i c t i o nt e c h - of p r e s - n i q u e on a s p h e r i c a ls h a p e .F o rs p h e r e s ,t h eN e w t o n i a nd e s c r i p t i o n s u r ei nt h es u b s o n i cr e g i o n exception. i s g e n e r a l l yv e r yg o o d , and t h .c a s e i s no The d i f f e r e n c e sb e t w e e nN e w t o n i a na n dt h ec u r r e n t l yp r o p o s e d t h e o r y on a s p h e r e i n t h e s u b s o n i c r e g i o n a r e v e r y s l i g h t , h o w e v e r , methodsgive andboth good a g r e e m e n t w i t h t h e e x p e r i m e n t a l d a t a . F i g u r e s 2 0 and 2 1 demonstratetheagreementbetweentheoryand e x p e r i m e n tf o r a b l u n t and a s h a r pe l l i p s e ,r e s p e c t i v e l y .A g a i n ,t h e r e is very l i t t l e differencebetweenNewtonianandthecurrenttheories,with eithergivingsatisfactoryagreementwiththedata. I n summary, t h e p r e s e n t t h e o r y o f f e r s good agreementwithNewtonian theoryforsphericalandellipticalshapes,butalsoprovidesaccurate d e s c r i p t i o n so fp r e s s u r eo v e rf l a t - d i s c - t y p es h a p e sw h e r eN e w t o n i a nt h e o r y f a i l se n t i r e l y .T h i s makes t h et h e o r yp a r t i c u l a r l ys u i t a b l ef o rt h el a r g e cone a n g l e b l u n t s h a p e s b e i n g c o n s i d e r e d f o r p l a n e t a r y e n t r y . 70 1 .c 0.9 0.e 0.7 0 n \ n 0.6 D A T AF R O MR E F ERENCE 15 C . 5I N C H E S 0.5 - 60' 2.283 INCHES 0.4 t 4 0 0.1 0.2 0.3 0 . 5 0.4 0.6 0.7 0.8 0.9 1 .o s/s* Figure 18. Pressure Distribution on a 60° Sphere Cone w i t h Sonic Corner 1 .o I - 0.9 &' _"_ NEWTONIAN THEORY M \ 0.8 I THEORY 0 \ I - - DATA FROM R E F . 1 6 - 6.0 0.7 I 0.6 e 0 0.5 n 0.4 0.3 0.. 2 0.1 0 0.2 0.4 0.6 0.8 1.o 1.2 1.4 S/RN Figure 19. 72 Pressure Distribution on a Sphere 1.6 1 .o 0.8 0.6 0 n \ n 0.4 0.2 0 10 20 30 40 50 60 8-DEGREES Figure 20. Pressure Distribution on an Ellipse 70 80 90 1 .o -. b 0.8 ~ - t M = 4.00 b/a = 0.5 Po = 1 . O ATM H, = 5000 B T U / L B 4 - 0.6 - THEUKY 0 a \ a "" - 0.4 0 NEWTONIAN THEORY DATA F R O M R E F . 1 6 0.2 " 0 10 20 L 30 640 0 -0 " 50 0-DEGREES Figure 21. Pressure Distsibution on an Ellipse 70 I 1 1 80 90 8.3 EVALUATION O F OTHER PROPERTIES Once t h e p r e s s u r e i s assumedequal i s known a t anygiven body p o i n t a n d t h e e n t r o p y tothestagnationvalue,othergaspropertiesofinterest may b e f o u n d b y r e f e r r i n g t o t h e M o l l i e r C h a r t t a b l e s S e c t i o n8 . 1 . as d e s c r i b e d i n Some c a r e i s n e c e s s a r y i n c a r r y i n g o u t t h i s o p e r a t i o n , The M o l l i e r C h a r t l o o k u p s u b r o u t i n e s h a v e b e e n w r i t t e n however. t o usepressureand e n t h a l p ya st h ei n d e p e n d e n tv a r i a b l e si nt h el o o k u pp r o c e d u r es i n c et h e s e a r e t h e known q u a n t i t i e s a t t h e s t a g n a t i o n p o i n t . I n o r d e r t o m i n i m i z e t h e error resulting w h e r ep r e s s u r e from l i n e a r i n t e r p o l a t i o n , p r o p e r t i e s a t o t h e r and e n t r o p y a r e known a r e f o u n d i t e r a t i v e l y body p o i n t s by e s t i m a t i n g e n t h a l p y ,l o o k i n gi nt h e tables a t t h e known P a n de s t i m a t e d h , determining e n t r o p y s , comparingwiththe known s , r e v i s i n gt h e h e s t i m a t e , and so on. T h i sp r o c e d u r ep r e v e n t s were t oi n t e r p o l a t e on t of i n dh . a compounding o f e r r o r s P and h t of i n d s, whichcouldariseif then interpolate P one and s 75 SECTION 9 THE TRANSFER COEFFICIENT &-I energyintegraltechnique EVALUATION is u s e d t o e v a l u a t e t h e n o n a b l a t i n g body shapewhich heat transfer coefficient along the surface of an arbitrary starts a t a s t a g n a t i o np o i n t . The e n e r g yi n t e g r a lt e c h n i q u e is c o m p l e t e l y d i s c u s s e d i n R e f e r e n c e 1 7 , and w i l l b e o n l y b r i e f l y summarize2here. f l a t p l a t e r e l a t i o n s h i p bei s assumed t o b e v a l i d f o r a l l In the energy integral technique, the tween energy thickness and heat transfer rate i s t h e nu t i l i z e dt os o l v et h ee n e r g y t y p e so ff l o w s .T h i sr e l a t i o n s h i p i n t e g r a le q u a t i o nf o ra na r b i t r a r yf l o w . ful to establish method, it i s use- To o u t l i n e t h i s a few b a s i c d e f i n i t i o n s . heat transfer coefficient: a CHI (104 ) p ue (Hr-h W 1 energy thickness : referenceenthalpy: h' = 0.23he + 0.19Hr + 0.58hw(laminar) h' = 0.36he + 0.19Hr + 0.45hw (turbulent) (107) where Hr = h e + ? :u V 2 2 2 = (pr 1 ) 2 = (Pr') ' I 3( t u r b u l e n t ) (laminar) (108) 1 (109) and a l l p r i m e d q u a n t i t i e s a r e e v a l u a t e d a t thereferenceenthalpycondition. 76 It w i l l b e r e c a l l e d t h a t f o r f l a t p l a t e f l o w s w i t h layer,theheattransfercoefficientexpression where a = 0.332 m = 0.5 or no i n i t i a l boundary is I laminar a = 0.0296 turbulent m = 0.2 Combining t h i s w i t h t h e d e f i n i t i o n which r e l a t e s h e a t f l u x t o of C H I r e s u l t s i n x for a flat plate. i s o t h e r m a lw a l l s ,t h ee n e r g yi n t e g r a le q u a t i o n Combining Equations111and F o rt h ef l a tp l a t ew i t hn e a r l y is 1 1 2 a n da s s u m i n gs m a l lv a r i a t i o ni n hw compared w i t h Hr-hwl t h e f l a t p l a t e e n e r g y i n t e g r a l e q u a t i o n c a n b e s o l v e d t o y i e l d r e l a t i o n s h i pb e t w e e ne n e r g yt h i c k n e s sa n d By cornbinihgEquations heatflux a x: 1 1 3 and 111, w e h a v e a n e x p l i c i t r e l a t i o n s h i p b e t w e e n a n de n e r g yt h i c k n e s s : Thisexpression i s now assumed t o b e v a l i d f o r p l a n a r b o d i e s o r b o d i e s ofrevo- lution,witharbitrarypressuregradient. 77 Themore generalboundrylayerenergyequationinintegralform Combining E q u a t i o n s 1 1 5 and114 , t h er e s u l t i n gf o r m ci. is .e i n t e g r a t e d d i r e c t l y to give (x 1 + [rPeue (Hr-hw)] A J x i where X i s a dummy v a r i a b l er e p r e s e n t i n gr u n n i n gl e n g t h .T h i sr e l a t i o nb e t w e e n 9 and x canbecombinedwithEquation 114 above t o g i v e t h e g e n e r a l r e l a t i o n x. b e t w e e nh e a tf l u x( o rh e a tt r a n s f e rc o e f f i c i e n t )a n d s t a r t i n g from a s t a g n a t i o n p o i n t , t h e f i n a l e x p r e s s i o n Forlaminar flows is 0.332p'p'rue(Hr-hw) PeUeCH = (Pr')2/3 '117) { l" S i m i l a r l y ,f o ra l lt u r b u l e n tf l o w , 0.0296 p ' ~ ~ ( " r ) l / ~ ( H ~ -h )lI4 PeUeCH = (Pr')2'3 I W ~ ' ~ ~ ( u ' ) ~ / ~ r ~ 5/4dX / ~ ( H ~ - h ~ ) F i n a l l y ,f o rt r a n s i t i o nf r o ml a m i n a rt ot u r b u l e n tf l o w I p'u'ue(Hr - + 0.0161 JX X t 78 (118) at x = Xt ' For a stagnation point, x * w e use the laminar expression and take the l i m i t as o with r = X = const. HZ-hw = const. p'p' The h e a t t r a n s f e r c o e f f i c i e n t r e d u c e s t o 1/2 Consistentwiththepressuredistribution stagnationpointvelocitygradientcanbe d i s c u s s e d i n S e c t i o n8 . 2 . 1 ,t h e e x p r e s s e d as whichbecomes InthepresentversiolL expressionE , quations of t h e c o d e , o n l y t h e l a m i n a r h e a t t r a n s f e r ( 1 1 7 ) and ( 1 2 0 ) , are i n c l u d e d F . or programming p u r p o s e s , i t i s most c o n v e n i e n t t o w r i t e t h e h e a t t r a n s f e r e x p r e s s i o n i n t h e form where u* * 4 "due R (point ds stagnation 79 or u* = U e away from stagnation point S J (Hr-hw)r2ueds p The running integral term for E* away from the stagnation point is evaluated by treating each linear term in the integrand it as varied if linearly bep', y ' , Hr-hw, r, tween body stations where its value is known. That is, and ue are all approximated as linear functions of s between each pair of of the body stations, thus making direct integration possible. The value integral at each body station is stored and merely added to the increment in R* at each new body station. All primed quantities are evaluated by referring to the Mollier Chart lookup routine with the known pressure and reference enthalpy. AI-I example of the accuracy of the energy integral method is provided in Figure 22, where results of the present theory applied to a spherecone configuration are compared with results from Aerotherm's multicomponent, 18). The two programs nonsimilar, boundary layer analysis program (Reference were supplied with identical pressure and wall temperature distributions Over the length of interest. It is apparent that the energy integral approach gives satisfactory results for the a sphere-cone type of body. nonablating heat transfer coefficientatfor least For use in the MACABRE code, the nonablating heat transfer coefficients computed with the analysis outlined above are corrected to account for the effects of ablation ("blowing").If we denote this nonablating coefficient A = 0) as peueCHOl then both data and simple analyses (see, for (that is, for example, Reference 7) indicate that the dependence of peueCH on is given fairly accurately by where cP= This correction MACABRE code. 80 is built into the 2xi 'euecH0 surface energy balance operations of the I 3.6 I I SPHERE-CONE I BODY 3.2 2.8 2.4 2.0 I u , "1 . 6 01 Q 1.2 0.8 0.4 0 0.4 0.8 1.2 2 . 8 1 . 62 . 4 2.0 S/RN Figure 22. Heat Transfer Coefficient Prediction 81 SECTION 10 OVERVIEW OF COMPUTER PROGRAM EVENTS The preceding sections have described the individual computing phases of the coupled computer code and the linkages between the phases. This secclarify further the tion presents two simple flow charts for the to code operation of the code. Figure2 3 gives an overall chart of code events, and of the surface energy balance sub-package. Figure 24 shows some details 82 ! I I N P U TI ,N I T I A L I Z A T I O N NODAL THERMAL R E S I S T A N C E S AND THERMAL COMPUTE NET HEAT FLUXESTO EACH SUBSURFACE NDDE (BUT RESERVING SURFACE NODES AND NEXT NODE DOWN F O RS P E C I A L I M P L I C I T TREATMENT) ~. CALL SURFACE ENERGY BALANCE PACKAGE TODETERMINE NEW SURFACETEMPERATURES AND R E C E S S I O NR A T E S( S E EF I G U R E2 4 ) COMPUTE NEW TEMPERATURESFORALLINDEPTHNODES,NOTINGSPECIALPROCEDURES FORSURFACENODES AND NEXT NODE DOWN t ~~ OUTPUT, IF (1) T H I S I S F I R S PT A S S (2T ) HIS I S A MANDATORY SURFACE STATE SOLUTIONTIME ( 3 )T H I S (4) I S REGULAR A OUTPUT TIME T H I S I S T HF EI N A T LI M E I F (1) OR ( 2 ) , V O I D OUT ONE OF THE TWO T I M ES H E E T SO F THERMOCHEMICAL DATA, INCREASEINDEXOFPRESENTSETOF 2 MANDATORYSURFACESTATESOLUTIONTIMES BY ONE Figure 2 3 . O v e r a l l MACABRE C o d e Flow C h a r t a3 CALLED . I FROM MACABRE 1 rn + COMPUTE ALLOWABLE TIME INCREMENT AS INFLUENCED BY SPACING OF TIME TABLE ENTRIES, SPACING OF MANDATORY SURFACE STATE SOLUTION TIMES, RECESSION RATES, OUTPUT TIMES * COMPUT CURRENT SURFACE SHAPE, S L O P E S 1 MAIN LOOP OVER SURFACE POINTS (COLUMNS) SET COLUMN INDEX =I 1 1 - m b DETERMINE TIME TABLE ASSIGNED TO THIS COLUMN, LOOKUP ALL TIME DEPENDENT QUANTITIES IN TABLE . i CALL PARBL (I), RETURNS ANY OR ALL OF CURRENT LOCAL PRESSURE, CURRENT LOCAL RECOVERY ENTHALPY, CURRENT LOCAL CONVECTIVE TRANSFER COEFFICIENT, AS REQUIRED ~~ DETERMINE CURRENT LOCAL HEATING OPTION: GO TO RELEVANT "PREPARATIONS" SECTIONS ACCORDING TO OPTION - - OPTION 2 NOT OPTION 3 ~ CONVECTIVE OPTION 1 SURFACE ENERGY BALANCE ~ ~ OPTION ~~~~~ v CONTINUED ON NEXT Figure 2 4 . 84 PAGE Surface Energy Balance Operations Flow Chart (Subroutine SURFB) I” - 1 IF HAVE JUST ARRIVED AT NEXT MANDATORY SOLUTION TIME, LOOK IN RELEVANT TIME TABLE FOR THIS STATION TO FIND FUTURE LOCAL PRESSURE, FUTURE LOCAL RECOVERY ENTHALPY, AND FUTURE LOCAL TRANSFER COEFFICIENT; ALL OF WHICH WILLREBE QUIRED IN FUTURE-TIME CALLS OF THE SURFACE STATE PACKAGE. IF SOME OF THESE QUANTITIESARE NOT ENTERED AT THE FUTURE TIME, USE SPECIAL RULES TO ESTIMATE THEM FROM CURRENT DATA ABLE DATA INCLUDING Tw THIS STATION TO TEMPERATURES LOWIN EST ASSIGNED BA ENTRIES IN TABLES. IF CURRENTTW isGREATER, GO TO BG SECTION OF TABLES AND ITERATE E;.ON IF CURRENTT, IS LESS, GO TO Tw SEC- ASSIGNED BA ITERATION SURFACE ENERGY BALANCE ITERATION SHOWN NOT t CONTINUED Figure 2 4 . ON NEXT PAGE (continued) 85 I 1 IF CURRENT B& IS ZERO HERE, ESTIMATE B& FROM BOTTOMBA’S IN TABLE I A CHECK TABLE ENTRIES SURROUNDING THE CURRENT ESTIMATE OF B& AND THE CURRENT TIME. IF THERE ARE ANY VOIDS AMONG THE DEPENDENT VARIABLE VALUES FOR THESE POINTS, CALL ACE PACKAGE TO FILL IN. FOR SOLUTIONS AT PAST TIME, USE SAVED PAST VALUES OF P and H, IN ACE CALL; FOR SOLUTIONS AT FUTURE TIME, USE PREVIOUSLY UP SET VALUES (OR ESTIMATES) OF P AND H, IN THE ACE CALL. + CONSTRUCT SURFACE ENERGY BALANCE EQUATION, NOTE “ERROR“ (DEPARTURE FROM ZERO). CALCULATE A CORRECTION TO B& BY NEWTONRAPHSON METHOD ~ BALANCED TO ACCEPTABLE TOLERANCE? I t CORRECT B& LIMITING CORRECTIONTO TABLE INTERVAL WIDTHS, REPEAT TO CONVERGENCE - I L ADJUST SIZE e OF SURFACE NODAL BOX ~ ON I I LAST POINT? * INCREMENT COLUMN INDEX: nI I + 1 I - I RETURN TO MACABRE Figure 2 4 . 86 SURFACE I (concluded) APPENDIX A TRANSPORT PROPERTIES FROM THE ACE PACKAGE When t h e ACE package of s u b r o u t i n e s i s used t o g e n e r a t e t r a n s p o r t property data for subsequent use in the transfer coefficient subroutine RUCH, the t r a n s p o r t p r o p e r t i e s a r e c a l c u l a t e d by ACE f r o m e x p r e s s i o n s which d e r i v e f r o m s i m p l e k i n e t i c t h e o r y a n d t h e p a r t i c u l a r m u l t i c o m p o n e n t d i f f u s i o nr e p r e s e n t a t i o np r e v i o u s l yd i s c u s s e di nS e c t i o n 5.3.3. The d e v e l o p m e n to ft h e s ee x p r e s s i o n s is discussedindetailinReference A b r i e f summary o f t h i s d e v e l o p m e n t , a n d t h e r e s u l t i n g e x p r e s s i o n s , a r e 19. presentedinthissection. DiffusionCoefficients - InSection tion for binary diffusion coefficients multicomponentdiffusion 5.3.3, a b i f u r c a t i o n approxima- was mentionedwhichcharacterized phenomena w i t h r e a s o n a b l e a c c u r a c y w i t h o u t u n d u l y c o m p l i c a t i n gt h es y s t e mo fe q u a t i o n st ob es o l v e d .T h i ss i m p l i f i c a t i o n a c h i e v e dt h r o u g h is a correlationforbinarydiffusioncoefficientsofthe form where 5 i s a r e f e r e n c e d i f f u s i o n c o e f f i c i e n t and t h e Fi a r e d i f f u s i o n f a c - t o r s .T h e s eq u a n t i t i e sa r ed i s c u s s e di nd e t a i li nR e f e r e n c e i n c o r p o r a t i o no f ( A - 1 ) i n t h e Stefan-Maxwell r e l a t i o n f o r f l u x e si n d i c a t e st h a tt h ed i f f u s i o nf l u xo fs p e c i e s terms o fo n l yp r o p e r t i e s of s p e c i e s i i 19. The mass d i f f u s i o n may b ew r i t t e ni n and g l o b a ls y s t e mp r o p e r t i e s . j e c t t o a few s i m p l i f y i n ga s s u m p t i o n s( R e f e r e n c e Sub- 19), t h i s e x p r e s s i o n f o r j i may b e w r i t t e n where 87 i s examined i n R e f e r e n c e 1 9 f o r a v a r i e t y o f I t i s shown t h a t t h e B i j c a l c u l a t e d by E q u a t i o n (A-1) re- The a c c u r a c y o f t h i s f o r m u l a t i o n chemicalsystems. present a very substantial improvementoverequaldiffusioncoefficients when compared t o exact v a l u e sc a l c u l a t . calculationofthemixtureviscosity on t h e d i f f u s i o n f a c t o r s g i v e n '1 d i r e c t l y from k i n e t i c t h e o r y . a n dt h e r m a lc o n d u c t i v i t y by (A-1) I a n d t h e s e The is based w i l l bedl-scussed in the following paragraphs. MixtureViscosity - The expression employed by the calculatethemixtureviscosityderives t h e o r y( R e f e r e n c e cussedinReference 20) I ACE program t o fromrigorousfirstorderkinetic subject t o a few s i m p l i f y i n ga s s u m p t i o n s ,a sd i s - 19. 1 I 'mix = i= 1 where ui is t h ev i s c o s i t yo ft h ep u r es p e c i e s i n terms o f t h e s e l f d i f f u s i o n c o e f f i c i e n t s where AZi The vi may be e x p r e s s e d Bii is a r a t i o o f c o l l i s i o n i n t e g r a l s b a s e d on a L e n n a r d - J o n e s i n t e r ( A - 1 ) and (A-7) i n t o (A-6) r e s u l t si nt h e m o l e c u l a rp o t e n t i a l .S u b s t i t u t i n g 88 i. p- following expression for the viscosity of the multicomponent mixture. (A-8) and t h i s i s t h e e x p r e s s i o n u t i l i z e d t o c a l c u l a t e t h e F i x t u r e v i s c o s i t y o u t put by the ACE program. - MixtureThermalConductivity a t o m i cg a sm i x t u r e may b e r e p r e s e n t e d b y ( R e f e r e n c e kmix = k where kmono-mix The t h e r m a l c o n d u c t i v i t y i n mono-mix + is thethermalconductivityin a poly- 20) (A-9) kint a m i x t u r e computed n e g l e c t i n g allinternaldegreesoffreedom andkint is thecontributiontothethermal c o n d u c t i v i t y of t h e m i x t u r e due t o t h e i n t e r n a l d e g r e e s o f f r e e d o m of t h e molecules. A s i m p l i f i e de x p r e s s i o nf o rt h e canbederivedin mono-mixturethermalconductivity a manner s i m i l a r t o t h e p r o c e d u r e p r e v i o u s l y d i s c u s s e d f o r t h em i x t u r ev i s c o s i t y .T h i ss i m p l i f i e de x p r e s s i o n is (fromReference r 19) 1 Xiki mono kmono-mix i=l xi + 1.475 (A-10) where k i is the t h e r m acl o n d u c t i v i t y otfh e pure species a liln t e r n adl e g r e e o s ffr e e d o m otfh e molecule. i n terms o f t h e pi i neglecting The ki may b e x p r e s s e d a sp e r (A-11) The c o n t r i b u t i o n t o t h e t h e r m a l c o n d u c t i v i t y f r o m t h e i n t e r n a l d e g r e e s o f 89 freedom may b e e x p r e s s e d as ( f r o m R e f e r e n c e 1 9 ) (A-12) By combining (A-1) kmix e with (A-9) through ( A - 1 2 ) , t h em i x t u r et h e r m a lc o n d u c t i v i t y as may b e w r i t t e n [ pl 1 5 xi R 4 Fi " x.F l i (6Ali i=l 1.475 +- 5 p1 - (A-13) where p 1 and vq a r eg i v e n L by ( A - 4 ) and ( A - 5 ) r e s p e c t i v e l y ,a n d (A-14) 2 T a. c P Thus, = zicpi (A-15) i= 1 t:he e x p r e s s i o n u t i l i z e d t o c a l c u l a t e t h e m i x t u r e t h e r m a l ( A- 13) i s conductivity output by t h e ACE program. A l s o c a l c u l a t e d and o u t p u t by t h e ACE program are t h e P r a n d t l a n d Schmidtnumberswhich are d e f i n e d h e r e a s Pr = L Cp-frozen k (A-16) (A-17) The t r a n s p o r t p r o p e r t i e s c a l c u l a t e d by t h e ACE program are a l l b a s e d on t h e b i f u r c a t i o n a p p r o x i m a t i o n f o r t h e b i n a r y d i f f u s i o n c o e f f i c i e n t s e x p r e s s e di n (A-1). casediffusion andthermodynamic 90 T h i s is so even f o rc l o s e ds y s t e mc a l c u l a t i o n s( i n which phenomena n e e d n o t b e c o n s i d e r e d t o c a l c u l a t e t h e c h e m i c a l s t a t e o ft h es y s t e m ) and f o r o p e ns y s t e mc a l c u l a t i o n sf o r which e q u a ld i f f u s i o nc o e f f i c i e n t sa r e assumed ( S e c t i o n5 . 3 . 2 ) . From t h e it may b e o b s e r v e d t h a t t h e p r o p e r t i e s c a l c u l a t e d a r e h i g h l yd e p e n d e n to nt h ed i f f u s i o nf a c t o r s , Fi.Three a l t e r n a t e methods f o r p r e s c r i b i n g t h e Fi were d i s c u s s e d i n S e c t i o n 5.3.3. The u s e o f t h e d i f f u s i o n equationspresented, f a c t o rc o r r e l a t i o n ( 6 3 ) w i t hr e s i d e n tv a l u e so f nref and E (whichwere d e r i v e dp r i m a r i l yf r o mc o n s i d e r a t i o no fs p e c i e sd i f f u s i o nc o e f f i c i e n t s )s h o u l d r e s u l t i nr e a s o n a b l ya c c u r a t ev a l u e so fo t h e rt r a n s p o r tp r o p e r t i e s . Altern a t e l y ,t h ec o r r e l a t i o n (63) may be u s e d w i t hv a l u e so f nref and E derived by c o r r e l a t i n g a v a i l a b l e t r a n s p o r t p r o p e r t y d a t a f o r t h e p a r t i c u l a r s y s t e m of i n t e r e s t . I f t r a n s p o r t p r o p e r t i e s o f maximum a c c u r a c ya r et o be c a l culated,thenthediffusionfactorsshouldbeinputindividuallyforeach s p e c i e s .T h e s e data may b eo b t a i n e df r o mt a b u l a t i o n ss u c ha sR e f e r e n c e 21. 91 REFERENCES 1. Bartlett, E. P., Rindal, R . A . and Kendall, R. M.: A Critical Evaulation of Recent Developments and Future Requirements for the Prediction of Ablation of Manned Reentry Vehicles at Superorbital Velocities. Vidya Division, Itek Corporation, Palo Alto, California, Vidya Report No. 1 7 4 , February 4 , 1 9 6 5 . 2. Aerotherm Corporation, Mountain View, California, Input Guide, Computer Program for Material Ablation with Chemically Active Boundary Layers in Reentry (MACABRE), September, 1 9 6 9 . 3. Moyer, C. B.: Axi-Symmetric Transient Heating and Material Ablation Program (ASTHMA), Description and User's Manual, Aerotherm Corporation, Mountain View, California, Report No. 6 8 - 2 7 , January 15, 1 9 6 8 . 4. Kendall, R. M.: An Analysis of the Coupled Chemicall.? Reacting Boundary Layer and Charring Ablator. Part V; A General Approach to the ThermoEquilibrium-Nonequilibrium, Homogeneous or chemical Solution of Mixed Heterogeneous Systems. NASA CR-1064, June 1 9 6 8 (also published as Aerotherm Report No. 66-7, Part V, March14, 1 9 6 7 ) . 5. Powars, C. A., and Kendall, R. M.: User's Manual, Aerotherm Chemical Equilibrium (ACE) Computer Program, Aerotherm Corporation, Mountain View, California, May, 1 9 6 9 . 6. Kendall, R. M., Bartlett, E. P., Rindal, R. A. and Moyer, C. B.: An Analysis of the Coupled Chemically Reacting Boundary Layer and Charring Ablator: Summary Report, NASA CR 1 0 6 0 , June, 1 9 6 8 (also published as Aerotherm Report No. 66-7, Part I, March 1 4 ,1 9 6 7 . 7. Spalding, D. B.: New York, 1 9 6 3 . 8. Moyer, C. B. and Rindal, R. A.: An Analysis of the Coupled Chemically Reacting Boundary Layer and Charring Ablator. Part11; Finite Difof Charring Materials Conference Solution of the In-Depth Response 1968 sidering Surface Chemical and Energy Balances. NASA CR-1061, June (also published as Aerotherm Report No. 66-7, Part 11, March 14, 1 9 6 7 ) . 9. Bartlett, E. P., and Grose, R. D.: An Evaluation of a Transfer Coefficient Approach for Unequal Diffusion Coefficients, Aerotherm Corporation, 6 9 - 5 0 , June 3 0 1 , 969. Mountain View, California, Report Convective Mass Transfer, McGraw-Hill Publishing Co., 10. Kendall, R. M., Rindal, R. A., and Bartlett, E. P.: A Multicomponent Boundary Layer Chemically Coupled to an Ablating Surface. A I M Journal, Vol. 5, No. 6 , June 1 9 6 7 , pp. 1 0 6 3 - 1 0 7 1 . 11. Rindal, R. A., Clark, K. J., Moyer, C. B., and Flood, D. T.: Experimental and Theoretical Analysisof Ablative Material Response in a LiquidPropellant Rocket Engine, Aerotherm Corpcration, Palo Alto, California, NASA CR-72301, September 1, 1 9 6 7 . 12. Browne, H. N., Williams, M. N., and Cruise, D. R.: The Theoretical Computation of Equilibrium Compositions, Thermodynamic Properties, and Performance Characteristics of Propellant Systems.U . S . Naval Ordnance TP 2 4 3 4 , NAVWEPS Report 7 0 4 3 , Test Station, China Lake, California, NOTs (ASTIA AD2 4 6 5 9 1 ) . June 1 9 6 0 . 92 13. Love, E. S., Woods, W. C., Rainey, R. W., and Ashby, G. C., Jr.: Some Topics in Hypersonic Body Shaping, (AIAA 7th Aerospace Sciences Meeting, New York, January 20-22, 1969) A I M Paper No. 69-181. 14. Boison, J. C., and Curtiss, H. A.: An Experimental Investigation of Blunt Body Stagnation Point Velocity Gradient, ARS Journal, 29, Vol.No. 2, February 1959, pp 130-135. 15. South, J. C.: Calculation of Axisymmetric Supersonic Flow Past Blunt Bodies with Sonic Corners, Including a Program Description and Listing. TN NASA D-4563, May 1968. 16. Belotserkovskiy, 0. M., ed.: TT F-453, June 1967. 17. A Study McCuen, P. A., Schaefer, J. W., Lundberg, R. E., and Kendall, R. M.: of Solid-Propellant Rocket Motor Exposed Material Behavior. Vidya Division, Itek Corporation, Final Report No. 149, February 1965. 18. Bartlett, E. P. and Kendall, R.M.: Nonsimilar Solution of the Multicomponent Laminar Boundary Layer by an Integral Matrix Method. Aerotherm Corporation Final ReportNo. 66-7, Part 111, March 1967 (also NASA CR 1062, June 1968). 19. Bartlett, E. P., Kendall, R. M., and Rindal, R. A.: An Analysis of the Coupled Chemically Reacting Boundary Layer and Charring Ablator. Part IV; A Unified Approximation for Mixture Transport Properties for Multicomponent Boundary Layer Applications. NASA CR-1063, June 1968 (also published as Aerotherm Report No. 66-7, Part IV, March 14, 1967). 20. Hirschfelder, J. O., Curtiss, C. P., and Bird, R. B.: Gases and Liquids. John Wiley and Sons, 1954. 21. Svehla, R. A.: Estimated Viscosities and Thermal Conductivities of Gases at High Temperatures. NASA TR R-132, 1962. NASA-Langley, 1970 - 33 CR-1630 Supersonic Gas Flow Around Blunt Bodies. NASA Molecular Theory of 93