Tobacco Documents Online

/ PREDICTION OF SECONDARY VORTEX FLOWFIELDS GENERATED BY AN INTERJ!ACTING MULTIPLE FR'EE-JET CONFIGIJIRATION by A. J. Baker,, J., A. Orzechowski, and' G. E. Stungis presented as Tecfinical Paper AIAA 83-fl289 at ttne AIAA 21st Aerospace Sciences ~!{eet~ing, Reno, NV, Jarnuary 1993, SSubmitted:, May, 1983 Revisedt March, 198'4 Synoptic: August, 1994 Department of Engineering Science and Mechanics University of Tennessee 1Knoxvillle, TN 3799b-203'0'

PRE'DICTION. OF SECONDARY VORTEX FLOWFI'ELDS GENERATED BY AN INTERACTING MULTIPLE FREE-3ET CONFIGURATION A. 3. Baker•, J. A. Orz.ecFiiowski* *, and G. E. Stcmgis+ ABSTRACT A penalty finite element numerical algorithm, f'or soiutiion: of the three- dimensional parabolic Navier-Stokes equations for subsonic turbulent flows, is examined for prediction, of- secondary vortex flowf ields induced by a multiple free-jet configuratiion. The combined action of decay of the initial high speed jets, turbulence level, induced entrainment from the: farfield,, and geometric discreteness for a four-jet configurafiion is predicted to produce a persistent system of eight counter-rotating vortex pairs in the plane normal to the jet axis. The magnitude of the inducedi transverse vortex vellociity components can approximate tu% of the jet initial velociAy. The results of a range.oi numerical predictions are interpreted and compared with available experimental data. L INTRODUiCTION! An important problem class in steady subsonic aerodynamics iss characterized' by the merging of' viscous andl (perhaps) turbulent unidirectiional flows fol'liowing abrupt termination o1 a surface off separation. The classic exampie is the confluence of upper and lower surface boundary layers at an airfoil trailing, edge. A second illustration, is a jet issuing into a quiescent *'IBM Professor of' Engineering Science and Mechanics, University of Tennessee, Knoxville, TN, Associate ffellow AIAA. *'* Principal Progra!mmer, Computational Mechanics Consultants, Inc., Knoxville, TN +Research Consultant, Brown & Willliamson Tobacco Corp., Louisvillle, KY. 1.

chamber or merging with a coflowing stream with distinct initial momenta. Each problem def~inutibn corresponds to abrupt, infusion of an (axial) momentum defect, the relief of which can result in generation of substantial perpendicular velocity components usually termed "entmainrnent " This is state& rnathematically by the continuity equatibn, which requires that an incompressible (or small Mach numberY velbcity (vector field) be always divergence-f ree. Two- dimensional or axisymmetr'ic definit'ibns for laminar flow and certain turbulent flows are amenable to exact andl/or similarity analyses, e.g., the inf'inite slot jet (SchlGchtingl')„ triple deck theory for the trailing edge problem (Melhick2) and/or an, isolated jet downfield of the potential core (Schetz3). These analysis techniques do not usually generalize to three-dimensiional! geometries, or to most two-dimensional turbulent nearfield flows of broader interest in aerodynamics. The characterization of'such flows generates the need to1 develop1 a suitable discrete: approximate (numerical) solution to an appropriately simplified form off the governGng, Navier-Stokes equation system. The abiding character of an aerodynamic flbw in this class is the predominance o2 a preferred (axial) flbw direction. An elementary extension of the convent'iional boundary layer order-af- magniRude analysis confirms that diffusion processes parallell to this preferred - direction are of the order Re Ys smaller than: all convection and transverse plane diffusion processes, where Re is the characteristic Reynolds number. The deletion of these terms yields the so-called thin-layer Navier-Stokes (TLNS) approximation, which corresponds to ai singular perturbation definition since the deletion of' the highest (second) order derivatives in, the axial direction removes ~I the ability (requirement) to specify a, downstream outflbw boundary condition. Provided a suitable approximation procedure exists to enforce the remaining, fully el'liiptic pressure coupling, the thin-layer equation system can be further

siimplifiedl to the "parabolic Navier-Stokes" (P;VS)i equations. A discrete approximate solution technique for the PNS definition utilizes the first order axial convection term to construct an iniitial-valve solution procedure. The PNS equation system remains an elliptic problem, definition in the plane transverse t+oo the axial coordinate direction, and definifiion, of a multiple-sweep viscous- inviscid interaction algorithm " permits an enforcerment of the subsonic flow three-dimensional elliptic pressure couplling, cf., Baker and Orzechowski4. As developed in this paper, a penalty finite element numerical sollution algorithm is well suited to construction of' discrete approximation solutions to a class of' subsonic, three-dimensional multiple free-jet problem definitions. The PNS ordering confirms that the continuity equation governs to first order the development of the transverse plane velocity distribution, while the two transverse plane momentum equations describe first order modifications to the otherwise uniform static pressure distribution. The sole known initial condition is typically the jet velocity (ratio, to the freestream), and'the farfield transverse plane boundary conditions must admit self -generation of a transverse mass flux distribution (corresponding to development of the associated entrainment field),. The key elements of the established penalty finite element algorithm applied to the free-jet problem class include, 1) establishment of a quasi-linear pressure. Pousson, equation, with complementary and particular solution fields that readily admit associated farfiel& boundary condition speci f iicat ions, and! 2) a specific form for the continuity constraint that provides an efficient and accurate procedure for gilobal communication of the local momentum defect relaxatiion.. For an initally turbulent jet, a parabolized' form of the two-equation turbuJent. kinetic energy-isotropic dissipation functiion, (4ty c) different ial equation system is employed, coupled with an appropriately ordered algebraic Reynolds stress 3'

c+onstit~ut~iveequatibn. This penalty finit'eelerment ail'goriitihm, is developed and app'lied to prediction ofal substantial secondary vortex fibwfield assaciartea with a eloseliy coupled multiple firee-jet problem definition. ll. PROBLEM' STATEMENT The three.-dimensional PNS equation system describing thesteady, subsonicc turbulent flow of an isothermal variable-density fluid to the principal scale of ordering4 is a j ~ J- _ 9 X ~ Ltw !Y= a.a,Xl~ 1 ~u !; l~ + d ~ t Lr k) = 3x k+~zR ~ k`°iC = 0 L(k) = a,x! [p u lk] + 2 XR I P uI k+'TP 2u! + pu~- +p£= 0. X a ~ L(c) = aX,! [PulJ + axI Cp utE + QCE €~Jt ax' 1 au 1 2 E 2 +CE Pulu~, k ax~ + e p~ = 0 4 a~ ~l+pu~ - e ax~ ak~ k- ~~ £ uj~ - v~ i ~.»I, / a x li (4)

The variables appearing in equat'ions' 1-S have their usual interpretation: in fluid mechaniks, and the superscript bar' denotes time averaging5. The tensor index summation convention is implied, with xl, aligned with the principal flow direction and' 1< j,c 3 an&2 c(k,,R) < 3. The turbulence kinetic energy k is the trace of' the Reynolds stress tensor,, c is the isotropic dissipation function, and Re is the characteristic Reynolds number. An algebraic Reynolds stress constitutive equation is employed to close equations 1-5. The PNS: approximation yields the significant order off components of u_u.' proportional to Re-Y1, and' the Reynolds stress modeli 1l simplif ied to, this order yields4 T- I> k3 a U l ulu,l, = Clk - C2C4-~ ax E k3 a'u1 2 "iui = C3k - C2Cw~ aX e 2 3 a,; C3k - C2C4 ~ a x3 uluZ = ulu3, _ u7u3 = 2 3u1 ,2, +(a~x1,, k3 au l a u' 1 - C2C4 -~ x2 ax3 (6) The various coefficients in equations 4-5' are modell constants with suggested' "standard"' valiues6 C = 1. C- 1= 1.3' C I = 1!.44 and! C2 = 1.92. The stress k ~ E ~ E E 5

rnodell coefficients Ca, l< a< 4, equation 6 are defined in terms of' two: empirical constants in the form 22(C401 - 1) - 60 C02 - S). C 1 - 33(C 01 - 2C 02) C2 = 33(C 0,1 - 2C02, Y C 44~C ~1 - 22C41 CL2 - t28iC4~2 - 4' = 1165(C~ i - 2C 02)2 C?., + 110 (7) where C'~1 = 2.8 and' Ct2'= 0.45' are suggested "standard" values, Hanjalic and Launder6. The resultant values for the a are { 0.94, 0:067, 0.56, 0.068 Summing the first three terms in equation 6 and dividing by two then yields 1'' _ 1 1t3 a!u li a7u 1 uiuj S i9 _ 2 (Cl + 2C3)k - C2C4 C Z _"~ a:x1 f . . . C (k ~ 2 aul aw l ..-,i ~ l~ 1.033 - O.OOS {£ ax , aX Jr s' ~. (9) The second term in equation 8 is uniformly non-negative and couplles transverse plane mean flow, strain rates into the Reynolds normad, stresses. Numerical experience 4'7 indicates the magnitude of' this term is at most a few percent off k; hence, equation 8 approximates an identity for the parabolic flow problem class. 4(3C 2 1') 11I(C ~ - 2C~ 2J - 22(C~1__1) - 1'2(3C 2 - 1). 6

I{ QI. NUMERICAL SOLUTION ALGORITHM In the primitive derived form, equations 1-3 do not represent a well-posed initial-boundary value statement for the subsonic flow problem class. As a consequence of' the PNS: ordering, simplificatibni both transverse momentum equations are independent of, the corresponding velocity component,, confirming the assertion that the continuity equation governs these effects to first, order.. Subsonic flow PNS algorithms can be distinguished by the manner selected' to establish a well posed numerical statement8+!'2, and each invol!ves augmentation of, equations 1-3 with the higher order (' d _ Re-Y, ) transverse momentum equation sety a~~ La~ k) -' ~ Ip ul,~]I' +' xR Pu~u k" Re ~ z~, '' 0 (9) The use of equation 9 requires inclusion of the appropriate O(d `) terms in the algebraic, Reynolds stress model, equation 6, which are (10) The penalty finite element aligorithm12 i's well suited to the f'ree-jet problem, dass boundary condition definitions. A tmansverse plane Pbisson equation for pressure is establishedl by forming the divergence of equation 3 plus equation 9, yielding I:(p ) = a / L(u k k i. 2- a a a a - _~+ aX~ aX~(p u~ + a xi ~"' iu ~ 2. 1 ~ -o~ - e 2 a Xi 7

For a turbulent flow, the second and third source terms in equation ll are negligibly small compared to the Reynolds stress term, while for a laminar flbw, their inclusion yields predicted (higher order) transverse plane pressure distributions that correctly balance mean flow eonvectioniand diffusion. The solution for the qwasi=linear equation l l is cast in the form p (xi) = pc(xi) + pP(xli,x,) ('12'). The complementary solution -pC satisfiies the homogeneous form of equation 11 with @irichlet farfi+eld boundary conditions determined by the farfield inviscid flbw. For the subsonic free-jet problem, pc reduces to a homogeneous constant. The particular solution p p satisfies the non-homogeneous equation 11 subject to a farfiield homogeneous boundary cond"utiony i6e., pp (x1,xR)= 0. The particular solution pressure field is stored at select axial stafiibns, during a PNS solution, for computation of, a,p%ixl for use in equation 2 during the next, PNS solution sweep. This multi+-pass PNS procedure has been termed "partially paraboliic," and satisfies the cited requirement of a three-dimensional elliptic pressure coupling upon convergence of p(xl) to a stationary f ieldi The second key aspect of' the penalty finite element aligorithm,, applied to the free-jet problem class,, is the functional form selected to enforce the continuity equation. Letting, u~ denote the semi-discrete approximation to u!i, the exact solution to equations 1, 2, and 3+9, a finite element algorithm employs 0 a weighted residuals statement to extremize the associated semi-discrete TJ ~ approximation error12. This algorithm statement for the transverse mormenturn. ~ ~ . 8 ~ ~.

equation must be "penaliaed"' to enforce discrete approximate satisfaction of' the continuity equation. The resulting theoretical statement is, , Lr R) f La (u R)~ +B'~ 2 a { Mk} L(p h) ={ 0} (13). t ! R t where the collurnn matrix {hlktxj)} is the associated finite dimensiional' subspace used'to define the approximation, i.e., u1(X i) z u~(xz,xl) _' U{Llk(x~)}T{ UL(xl)}e Cl4) e In equation 14, the subscript "e"' denotes pertaining, to the finite element domain Re, and the elements of {' UL(xl)} e are the x}-dependent nodal values of uh on the solution domain discretization nh = U'R e x xl. The, functional, form o2 equation 13 is determined as a direct extension of the cllassical analysis for the (linear) Stokes problem, see reference 1'2. Ch.5. i However, direct replacement of' L(p h) with V• p u in equation 13 would be unsatisfactory for the free-jet problem dass, since V hi is only locally supported using the conventional truncated polynomiat basis { Nk},. For a jiet' (or wake) flow, the momentum defect relaxation process is highly intense within the immediate vicinity of the initiation process,, and the intensity decays roughly exponentially into the si'milarity, region. The penalty functional' form must therrforef~acil'utateglobal communication (t'~hroughoutAZ'2, thet'ransverse plane)~ of' local momentum modification processes. This req}uirement is satisfied by definition of the measure of, LGphy by an elliptic boundary value problem on R2. The obvious choice is the harmonic functiion, ~ satisfying the Poisson equation.

t(~' )= a'-2L _ 2 . az. 0 u iY = 0' (15) 2x~ r A first integral of' equation 13' confirms tham ~ exhibits the character of' a scalar velocity potential. Therefore, the boundary condition statement for 4~ on the boundary a!li;' of Tt2 is, ~ (~ ) = a~ + b ~I •~iLI= 0 (1~6). where a= 0 yield5 a non-through flow constraint' while b= 0 defines a porous boundary. Since ~ is arbitrary to within al constant, 0 = 0 may be define& anywhere on a boundary segment where b= 0. Hence,: for equations 15-16 ~h -> 0 in any norrn as p~ approaches a divergence free field. The implementation of ~, h as, the penalty measure is embedded within the iteration procedure associatedl with the PNS: algorithm solution statement. For the remaining initial-valiue dependent variables (u i, kh, Eh), as wel4 as the boundary-value variables (pP; f h)„ and letting qh denote the generalizedi dependent variable set„ the weighed residuals semi-discrete approximation error extremizatiion statement: is12 JR2kt { N (x d} t(q) = { 0 } ('17) Equations 13 and 17, for u h kh and c h each produce a coupled ordinary differential equation, system of the form, IA(~1) dl{ ~i} r { g(qi)) = { 0~ }' li N 10

which is employed to complete the Taylor series on, the interval xji l- xj = e x1, {Fl } E {QI}jil & x1,°LV,. (19) _ { p} dxl j+14 to define the non-linear algebraic statement on the discrete approximation {Ql(x1)} (where I is the discrete index denoting,members of qh) For the Poisson variebles p~ and ~h, equatiorn 17 produces { FI} _{ 0} directly. The soliution of the combined non-linear algebraic equation system is obtained using a modified Newton iteration alg,oritfirn of the form 13QQIJ,I ' I p 1{'~~QI }'~+ 1=-{ Fl'~} P I In equation 20,, p is the iteration index, and the p+ist discrete approximate iterate is obtained in terms of'the solution as {QI},p+l ' {'Ql}J+ 1 + {aQ{}j+l (21) For the PNS' algorithm, implementation1'3, the Newton algorithm Jacobian [' J( {FI })), is replaced by two sparse matrices. The initial-value dependent variables u h lrh, andi E h, i.e.,, {Ql }; i< I< S, are solved as multiple right side substitutions to equation 20 using the u i Jacobian [ J I 1) , where [Jl'ljl = a{' F1}' R7_qff N The f ield variables p~ and (P h, Le., { Qd }, 6 < 1< 7, are thereafter solved as ~ rnultiple right side substitutions using the p Jacobian ( 3661 . The six P 0411 components ofuj~uj, are subsequent'lycalculated using an adgebraicassembly~ 11 4

procedure eqwivalent to solving equation 20 using [ J88] with multiple right side substitutions. The PNS ~ algorithm timing uti'izes Ahe sequence (31. 1] ,[ J66j , ,and [398], , with update of the non-linear Jacobian at' each iteration. Hence, the pth level of the Poisson field 4 P' is calculated using the pth approximation to the velocity field' uP. This predictioct is'added algebraically to the previous p-1s4lutions to construct the functional form of, the penalty term in equation, 13 as I _ 61 2 ~ {Nk) E f: h (23) l i'_ 1 R !f a. R2-£ Hence, each iterate f ~, modifies the previous p-1 contributions to the.penalty term, as the correction requiredl to move the next iterate puprl towards satisfaction of the divergence-free requirement. In the limit as p becomes large, Vh- p+I approaches zero in some measure, hence II ~~+li 0 due to: the definition statement, equations 15-16. In practice, following a few extra iterations to homogenize the inital-condition discrete approximationerror, the penalty f irnite element PNS algorithm converges to VQI Irnax < 3. x 10-4 in - typically 4-5 iterations, for S z L in equation 23. The intrinsic measure for ~, h P is the energy semi-norm,, a~h' aO h 11~pl~~F = ~ R z ~~ ~ (24) and typically 114 p 11E < OC10-S) at convergence. This measure can be interpreted in terms of, kinetic energy, since o~A, P b p(ul- u~) is a direct measure of discrete approximation error. For comparison, the kinetic energy of' the mean: flow is «100), and the kinetic energy o,f the turbulent kinetic energy field is 12

0(10 -2'): Hence, the discrete approximate energy error in the velocity field (measure) deviation from exact divergence-freeness at convergence of the Newton iteration is nominally insignif icant. IW: DISCUSSION AND RESULTS The penalty finite element solution algorithm for the PNS equation system, as operational in the CMC:3DPN'S computer programl3-1S, is well documented for a variety of subsonic two- and three-dimensional aerodynamics problem definitions. Within the inviscid interaction, framework, results are published for laminar and turbulent flows in a three-dimensional juncture region geometry and in, a square cross-section ductk. Data on convergence history and satisfaction of the continuity requirement.are presented, as well as detailed, quantitative comparisons between experimental'. data, and PNS' prediction of distributions of the Reynolds stress tensor, equations 6 and 10, for a flow witti a secondary vortex structure. Similar comparisons are reported7 for a turbulent two- dimensiional airfoil~ wake geometry, including, excellent quantitative agreement with experimental data for mean veloci'ty profiles and Reynolds stress distributions over the entire strong interaction region. The robustness of the penalty algorithm for the subsonic PNS interaction problem, class is thus well substantiated, as well as the pertinence of the algebraic Reynolds stress clbsure model implementation using publishe& correlation coefficients. The specific multiple free-jet problem dass of interest is graphed in Figure 1. Four jets, symmetricalPy disposed about the circumference of a circle of radius r= R, exhaust into the half space xli/R > 0 with an initial veliocity (dustribution) strictly parallel to t'he xl-axis. Dependent upon design parameters, this four jet system can, undergo a rapid loss of' initial momentum, with the 13'

consequence that a substantial transverse plane velocity f ield must become induced. The close proximity of the initial jets yields dominating interference effects, such that' the resuDtant three-dimensional velbcity field exhibits (according to the PNS sollution) a systematic secondary vortex distribution. Figure 1 also graphs in perspective the PNS algorithm prediction of a typical induced vortex flowfieldl evolution, in a symmetric quadrant of the transverse plane on 0:25 < xl/>;t' < 1.5.. The engineering, design requirement for the multi-jet configuration is to induce a rapid mixing, ofl fluid, initially contained within the interior region r < R, into the exterior region r<Pt. Experimental data conf'ums tham t'his multi-jieR configuration operating at design conditions efficiently accomplishes this requirement6 A video recording, of the device in operation, and with the multi'- jet system shut-off,, was shouvn during the oral presentation16 at the AIAA 21st Aerospace Sciences Meeting. Using smoke tracer photography,, these visual comparisons provided strong qualitative evidence of the robustness of the mixing, with, the jets in operation. Conversely, the retained' integrity of the smoke tracer column verifies the persistence of a unidirectional fliow with the jiets inoperable. Figure 2a is a still frame from the video verifying this persistence of the smoke tracer column, with the multi-jet system off. Figure 2b is aa corresponding, stilli framewit'~h the jjetsystem, operating. Eventhough the inlet flow conditions for both tests are strictly unidirectional (to the right, parallel to the experiment' axis), the annihilation of the column of smoke appears almost complete with, the multi-jet system operating. If' the mixing process was diffusion dominated, a well defined spreading strvcture would be evident. That it is not indicates that a convection process must dominate, i.e., the multi-jet 14

system must self-create a substantiall multi-dimensional velocity field from a unidirectional inlet flow conditiion.. A measure of the magnitude of' this generated velocity f'ieldl can be estimated in Figure 3„ which is a close-up view of the strong, interaction region! immediately downstream of the plane xl/rt = 0, see Figure 1; and 2b). The initial smoke column is uniformly distributed' on r < R at, xl/R = 0 for all experiments. With the multi-jet system operating, this column is eroded tola conical focus by xl/R :r 1', whereupon stringers of smoke appear to, be propelled outwards and'downstream aibng planar trajectories. Assuming the flow is steady, these streaklines must identify corresponding streamtubes, hence the corresponding tangent angle is a measure of induced flow angle. Defining p as the angle between the xli axis and! a streakline trajectory, Figure 3 yields -39° < 0< 27° as the range of' flow angularity. WhiFe this does not quantify the nature of the flowfield induced by the multi-jet system, it does give strong evidence that it must be substantially three-dimensional. The requirement of the PN'S analysis is to quantitatively predict the inducedl flowfield and to characterize important design parameters. The characteristic dimension of the multi-jet problem is R = 4mm. The jet flows, and the co-parallel interior and.' exterior (r > R) inlet unidirectional flows, are initiated by, applying a reduced pressure everywhere on the half-space x1JR > -1. The resultant induced initial condition velocity field is, ~ uo(xt +xl/D = 0) = 5o(x~,x11LD = 0)a + Oj + 01~ (25) 15

The design initial' velbcity for eacK jet is ua' = 12 rn/s, while on the interior regibn~ (r < R), uo x 0:25 m/s, andl im the exterior regiom (r > R), u o= 0.1 rnLs. Therefbre, the jet' system velocity ratio is X> S0. The initiali cross-section of each jet is non-circular,, as formed by the region interior to the approximate bisectiion ol a nominal 1 mm diameter hole by the circle r= R, see Figure li. This yields a hydraul'iic diameter dh z I mm4 Since each jet flow channel! thus contains two interior, corners approximating right angles, and each is approximately 20 mm in length and has quite rougta walls, and the jet Reynolds number is Redh z 0:6 x 104, the initial jets probably possess a turbulent structure. Hence, the requirement exists to estimate ko(x,,xl = 0) > 0 and o(xI,xl = 0) > 0, for the PNS analysis,, whilk everywhere exterior to the initial' jets ko Z 0 : eo. This completes identification of' the initial-conditions that must be specified for a PNS analysis. The other requirement is to select an appropriate region of R2 for application of boundary conditions. Computations could be executed on the entire transverse plane, but symmetry in the initial conditions andl geometry permits use of a quarter-plane, see Figure 1. The transverse plane region used for most PNS computations is the domain 0 < x./R < 2.0, as graphed in Figure 4 with boundary O-A.-B-C. A non-uniforrn discretization U'Re was defined to permit "adequate" definition: of the initial jet region using: an M = 19 x 19 mesh,, the nodal lbcat!ions of' which are shown as dots. For later comparison with tlhe smoke streakline datay an, inert' species mass fraction was initialized at' all nodes of the interior region, denoted as "smoke column"' inN. Figure 4. The corresponding PN'S algorithm boundary conditions are summarize in Table 1. 16

Tablle 1 PN5' Algorithm Boumda;ry Conditions for Multi-Jet Symmetric Quarter-Plane Domain, Figure 4 Boumdary Boundary Condifion Scgment Vanishing Neumann Diridhlet Comments O- A u I, a 3, k,, E, pp, m u 2= 0 Non-porous O-& ul,u2,k,e,pp,~ u3=0 Non-porous A- C u! i, u 2, u 31 k, ~ ~= 0' = pp Porous B - C ul,u2,u3,k,4 ~ ~ =0=pP Porous While the basic penalty finite element al'goriRhm is well' verified,, benchmark tests of specific pertinence to the free jet problem are required as relates to initial-condition specification for k andi e and flow symmetries. The simple two-dimensional slot flow geometry,, Figure 5, serves the requirement. Figure 6 summarizes initial cond'utions, and the penalty algorithm prediction for two representative downstream distributions of ur 1, u 2, k and~ E, for an assumed symmetric slot jet flow at velocity ratio a= 50. The initial conditions are uo = 10 mi/s on 0< x2/Hf < 1 and w,o = 0:2 rn/s on 1< x2/'H f< 2.5, where. Hf = 5: cm is the slot half-width. Vanishing normal derivative boundary conditions are appllied'. for all variables at x2/H f = 0 and x2/H,f = 2.5,, except ~= 0 at x2/Hf, = 2.5, and an M = 50 element uniforrm, discretization is employed. No other specificat'ions are required for a liaminar flow siinulation,, but both ko and o distributions are required for a turbullent flow simulation. In the absence off' definitive detailed experimental data, one procedure to self-ini'tialize levels is based on an order of magnitude assumption for a"turbulent eddy viscosity"' vt = C4 k2/'e, see equation 6. For "parabolic" aerodynarnic flows, the extremum level of ko (non-dimensio,nad), may be assumed' of' the ordler 17

O('10 `'), i.e., the kinetic energy in the turbulence field approximates a few percent (at most) of the mean flow kinetic energy. Since C4, is a constant,, defining v~a_' O(10p); where v is the (laminar) kinematic viscosity and (say). 1< p < 3,: defines a corresponding leveJl of eo Secondly, since 1/e appears throughout the turbulence closure equation system, it is preferable numerically to define a smalll but non-zero farfield level. Numerical experimentation~ with the slot- and free-jet pr©blems confirmed the penalty PNS algorithm will rnaintain levels of co/'Re z O(',10-9) and k° c (X10-4), with vanishing normal derivative boundary conditions, yielding vt/v z 1 in the farf'ield as an approximatiion to a nonrturb;ullent region of' flow. Figure 6 highlights the essential aspects of the PNS algorithm prediction for the two-dimensional slot-jet flow for the specific initial' conditions k~ = 0.05, vo/v = 10 and A = 50: The initial distribution for ii ° is assumed a step function on the nodes of' U Ite; hence, non-zero levels of' ko and eo are correspondingly defined at one node only. The PNS algorithm is marched downstream five steps, using A x l,/H f z 10-3, to smooth somewhat these initial conditions, whereupon thee penalty algorithm for ~h is initialized. The energy norm for 4~ h reaches O00-5) in approximately 110 more integration steps; hence, the entire algorithm iss furnctional at A xl/'H1, -` 10 2. Approximately 150 (progressively larger) integration steps are required to reach xl/'H f= 0.5; where the extremum levelss of k and € are computed, Figure 4c)-d), using, the standard closure model constants. Thereafter, the level of' both variabUes monotonically decreases as the jet spreads liaterally, with the PNS algorithm reaching, the final station xl/H f = 1.0 after 50 additional integraRibn, steps: The PNS integration step size would continue to increase in proceeding, further downstream as the solution field becomes progressively smoother. Changing the initial conditions, to say vtlv - 16A 18

li0?, would not alter siignifiicantly the solution appearance, but woutdl simpl'y compress the history over a shorter axial displacement due to the more vigorous initial mixing levels. Since the multi-jet problem of interest is not of circular initial cross- section, see Figure 4, and the resultant flowfield is not axisymmetmic, Figure ll, a second verificatiion, test must evaluate a three-dimensional problem of a nominally circular cross-section jet interpolated onto the nodes of a non-uniform rectangular discretization. Figure 7 graphs the M! = 13 x 25 symmetric half- plane d'ascretization, withispan 6D x 2D; for a circular jiet, of initial diameter D= 0.05 mL At nodes interior to the circle r= Di/2, the discrete initial specification for jet velocity is wo = 30 m/s, while uo = 3 mi/s everywhere on the exteriior, r> Di/2, yielding a jet initial velocity ratio a= 10. No other initial conditions are required for a laminar flow simulation; for a turbulent flow, the initial condition specifications at nodes on r < D/2 was ko = 0.01 and v©/v = 110, while on the exterior, ko = 10-4 andi v~/'v = 1. The three farfield boundary segments: are porous, while the jjet bisector is a symmetry plane, recall Table 1. Figure 8 graphs the computed transverse plane velbcity distributions at xl/D = 0.25i 0.5' andl 1.0, for the laminar flow simulation; Figure 9 presents comparisoni data for the turbulent flow case. The direction of each arrow is parallel to the local streamline, and't'he length is proportional to tfie magnitude J scaled to the lbcal non-dimensional extremurn val'ue u~'_ ~ u~ax I / The t~. initial jet diameter is also noted. Both PNS solutions exhibit an excellent approximation to nominal radial entrainment, except that by x1/D: = 1'.0,, the laminar flow prediction exhibits lbcal perturbations about the jet boundary.. These are induced by the step interpolation of the iniTiall data, coupled' with, the 19

relatively srmall' magnitude computed for u R; see Table 2, which agree well with. the order of Re~ - 10-3. While little visual distinction exists between the laminar and turbulent solutions,, since the vector plots are scaled on the local maximum, the magnitude of the entrainment field for the turbulent simulation is 20-30 times that of the laminar flow, see Table 2, and a distinct maximurn is predicted at xl,/D = 0.5: Table 2'. Local Predicted Entrainment Velocity Maxirna Circular Initial Cross-Sectibn 3et, Uo - 30 m/s t, Axial Station Transverse Velocity Maximum, u R xI/b Laminar Flow Turbulent Flow 0. I 0. 0. 0.2'5' 0X015 0.038 0.5 0.0014 0.058 0.75 0.0015 0.052 1.0 0.0017 0:044 1.25 0.0018 0.034. One additional feature in Figures 8-9 warrants comment. The largest transverse velocity components are predicted directly adjacenro to the jet periphery, with the flow directed radially inwards on r y D/2, radially outwards on r< Dl2,andl vanishiragly small withini the interior regi+on~ (jet potential core). This distribution is the direct response to edge erosiion, of the potential, core, hence au 1/axl < 0, while the adjacent acceleration of' the coflowing, exterior stream yields 3&1i/'axl > 0. In two dimensions, the continuity eqWation respectively yields a1Z 2/ax2 > 0 and au 2/ax2 < 0. Figure 6b) ciearly confirrrnsf this action for the slot-jet probleran„ which is illustrated in Figures 8-9 by the counter-directed vectors about the cirde of diameter D. The majority of P'NS analyses for the multiple jet problem, Figure 1, were conducted using, the symmetric quarter-plane discretizatiion, Figure 4, and 20

boundary conditions defined in Table 1. As noted, each initial jet'probably contains a turbulent structure, hence l~'> 0 and E~ > 01 must be defined. Since the jets are of small initial cross-section,, the potentia[ core defined by the initial conditions erodes very rapidly for any vjS) < 5;, hence the secondary vortex flowfield develops very rapidly, reaching a mature structure within 3- 4 mm dbwnstream from the injection plane, see Figure 1. This extremely rapid evol'utaon creates a numerical stability problem for the PNS penalty algorithm, since the exterior co-flowing stream uo = 0.1 m/s was of such small magnitude. (The matrix [A('-)] in equation 18 contains u i(x, ,xi), and itsessenti+al inverse is required to evaluate { Ft}, equation 19.) This stability problem was circumvented by adding a uniform background velocity u. to~ the initial condition distribution for ~ andl ua. This preserves the initial condition defined velocity strain rate distribution on i1;2, for any u.., at the expense of translating the axial! prediction station by some ratio of um to ul(xl). Nurnericall experimentatiion~ confirmed that um/u~ > 0.1 yielded a stable integration, and uJb?= 0.2 was utilized for allltests. Therefore, the PNS,algorithm initial condirtibnspecification is uj = 1.2(u0 ) _ 14.4 m/s, uo =(Oi25 + 0.2' uo) = 2.65 m/s and u~ = (0.1 + 0.2 uo) = 2.5 m~/s, for bot'~h, turbulent and laminar flow, ppredictions. For the turbulent flow base case, the initial condition assumptions were ko = 0,005 and vo/v = 35 at all nodes inside the jet see Figure 4, and' ka = 0.0001 and vo/v = 3 everywhere else. Figure 10 graptas, the PNS prediction of the resultant nearfield transverse piane velbcity distribution on 0.15 < xli/R < 1.5, i.e., from 0.6 to; 6 mm downstream from the initial condition plane., The q}uarter circJe locates r = R = 4 rnm,, and, the initial jet cross-section is sketchedl in Figure l0a). As in Figures 8'-9, for clarity the vectors of eachl plot are scaledi om the local predicted maximum u~ noted in the 21

legend. At x 1/R - 0.15, Figure 10a), u m ~= 0.085 arnd': the velocity field exhibiRs a nominall radial distribution analagous to that predicted for the discrete circular jet. Some minor grid}induced perturbations are evident,, recall Figure 8b), but these disappear by x1/R = 0.5, Figure 10b). Over this interval, the local' maximum velocity has decreasedl substantiallly to uM' = 0.058, and the radiali pattern now exhibits a bifurcation and' imminent initiation of a vortex pair about the jet at r= R. The maximum transverse vellacity remains nominally constant to xl/'R' = 1.0,, Figure 10c); but the single vortex has noNV become a double pair with foci symmetrically displaced about r = R and ad jacent to the jet. This double vortex pair matures by x1/R = L.5, Figure 10d), andi has increased somewhat in strength to w~n = 0.067. The PNS algorithm predicts this double vortex pair configuration to persist nominally unchanged, to xl/R = 11. Figure 11 is a composite assembly of' the predicted transverse plane distribution u~(xt) on the entire domain span at, x1i[R = 1.5. The four vortex pair systems external to r = R reach to the solution domain boundary aR, yielding the corresponding prediction of influx/efiflux distribution oni M. Thee four vortex pairs generated internal to r = R' contain the largest' velocity magnitude, and penetrate almost to the center. Recall that the device design requirement is to quickly mix the fluid interior to R into the external regpon.. The rrnulti-jeC (decay) system is predicted by the PNS solution to accomplish this by self'-generating, a vortex field distribution of' sufficient magnitude tqo penetrate to the core of the region and then convect this material across the boundary r = R' to be transported away by the external vortex pair systems. This PN5 predict'ion, seems pl'ausible in describing a mixing mechanism that is sufficiently vigorous to corroborate the video data. The requirement is tp, 22

1 verify the impactl of the assurnptions made principally regarding the initial' turbulence level6The eZore, PNS executions wer repeatedi for the initiial- condirtiorn assumptions km = 0.0025, ~o/v = 9, and 19 = 0 hence ~a = 0, i.e., laminar flow. Bottn solutions predicted generation o2l a paired vortex field! in the transverse plane. The caorresponding plots of the transverse velocity vector nearfieldf distributions at xl/R = 1.5' are shown in. Figure 112, and are now scaled to1 the local maximum of the~ base case NIP = 0.067) of Figure lOd). The results obtained for the half-initial turbulence level are visually indistinguisfiable from the base case, although u~ = 0.059 is about, 15% smaller. The l,iiminar flow prediction yields u~= 0.009, an 85% reduction from the base case, and results in generation of a single weak barely-distinguishable vortex pair with centroid at each i side of the jet on r = R. Figures 13-14 summarize the results of' these PNS predictions for jet axial velocity decay and maximum transverse velocity magnitude. For the laminar flow prediction, the jet stiil1 retains its potential core at x1;/ft = 1.5, Figure 13, m since u 1 - ~ ulmax'/% = L Nevertheless, the edge erosion of' this core does induce a entrainment velocity field throughout R2, Figure 12b, with an extremum magnitude u~ = 0.014 prediicted. at xi/lt = 0:4, Figure 14. In distinction, u~ i< 11 by xl/'R ' - 011 for both, initial turbulence level assumptions, Figure 13, and the essential effect of the smaller initial assumpt'ion! is to displace the jet decay curve downstream by the nominal distance Axi,/R ~~ 0.5, Figure 13. Figure 14 confirms that both initial turbulent, assumptions produce an extremum in u~ on xl/It' < 0..25,, and that, the smaller k0 yields a proportiona~lly smaller maximum. Thereafter, however, the curves for u,.rn are nominally parallel and the vortex patterns are qualitatively identical, indicating, that the initial condition assumption (error)' appears of small consequence in the comparison. 23

This may be d'irectly aRtribute6l to use of the k-C difterential equation dosure system, wherein local source and sink mechanisms exist to rapidly, adjust the initial field to the mean flow, strain rate distribution. Table 31 summarizes the comput+ed' maximum turbulent~ eddy viscosity vt = C4 k2/e , non- dimensionalized by V, for the two turbulent simulations, which quant.i'fies the rapidity of the initial condition adjustment process for the mult'a-jet problem def iniRiion. Table 3 Evoltrtion of Maximum Tbrbulent Eddy Viscosity Axial Station Maximum Eddy Viscosity Vt/v x 1/P' ko = 0.00'5' k° = 0.0025 0. 35. 9. 0.25 64. 25. 0.5 52. 31. 0.75 42. 28'. 1.0 38. 25. 1.5' 33. 21. It still remains to assess qualitatively, if possible, that the base PNS turbulent flow prediction correlates with, the experimental data, which is limited to smoke tracer flow visualization. For this purpose, an inert fliuid species masss fraction equation L(Y) was added to the PNS set, of, the fprm L(Y) _~ L~ l Y'' 4~ x.. ~ uQ Y- v~1 II = 0 (26) 24

The boundary condition for Yhi isvanushing normal derivative on O.-A-B-C, Figure. 4, and the initial condition is {Y'(xl = 0)} =(1 } everywhere interior to r= R' excluding the jet. For orientation, Figure 15 graphs a perspective composite view of the initial conditions for u i and Yo, on the solution iniiriation, plane xl = 0, and superimposed on Figure 1 Sa is the locator radius r = R. Figure 16 graphs in composite perspect!iwe the. PNS predicted evolution of the mass fraction distribution on 0< xi/R -< Ll.w (wherein plotting of the base plane zero level has been suppressed for clarity:), The tracer levelI is observed to first be preferentially removed from the innrnediate neighborhood ot each jet, Figures 1i6a)-b), yielding four ridges duminating in a, peak by xl/R = 1.0. Thereafter, the interior region vortex field, Figure 1'1!, contiinuously extracts the tracer from the interior region, (r < R),, and the exterior vortex field propels the extracted material out of the solution domain bourndar allbng radial rays nominally bisecting each jiet.. These radial spokes of tracer material may correlate with the smoke flow filaments observed in the photographic data, Figure 3. These data were employed before to estimate the range -3'9iD < 0 < 27° for the velocity tangentt vector. Figure 17 is a graph of, the extremum computed: velocity vector tangentt angle 01, on r<' R, and 0e on r > R, as computed for the PNS base initial turbulent simulation on 0 < x l/R < 6 according to the equation _ ~ u a'e , . (x _(sg'n uI 'rn~ax~ tarn 1 Rmax 1 (26) j ul ma,x(xl) On 0 < xl/R < 0.25, the extiremum u~' and ue were counterdirected in the immediate vicinity of the jet, recall Figure 10a). Thereafter, the local maxima are directly adjacent to each other at the intersection of the jet edge with r = R', 25'

Figure 1'0d). The data in Figure 1!7 indicate that by xl,/R = 2,, the maximum: velocity vectoc tangent' angle on r< R approximates 250, which is certainly within the range of' the data, at' xl/R = 2, Figure 3. The PNS maximum tangent angle increases to ~l'e . 400 by xl/R = 6, mainly in response to the monotonicc decline in the denominator of' equatiam 26 1 recall Figure 13. These comparisons certainly tend to further substantiate t'he assertion that the PI*1S1 algorithm~ sollution has yielded a qualitativel'y and quantitatively validl estimate of the rrnulti--jet geometry induced flowfield. As stated, the multiple dual vortex pair transverse plane velocity distribution,, as illustrated' in Figure 11, is predictedl to persist nominally unchanged to xl,/R = 11. However, the PNS algorithm begins encountering aa stability problem by xl,/R = 10, since aR this axial l'ocation the magnitude of the transverse velocity is, essentially equal to the axial velocity, in violation of the basic PNS ordering assumption. Insipient instability was detected by t'he penalty algorithm at xl/R = 8, whereupon 11 10 Q11E, no longer exhibited a monotonic decline with the iteratiion index p. However, the level of, 11~ ~E at', convergence remained O(10-5) to xl/'bt = 10; whereupon it began to steadily further increase as xl/R ' increased. Hence, the energy norm of the Poisso variable ~'hi is an excellent' measure ol solution robustness, and can be usedl with confidence to predict when, the theoretical assumptions are becoming violated. This is considered' an important feature of the penalty PNS aligorit'hm, for application to the free-jet problem dass, In this regard as well, the computational cost of the PNS simulatiiion are modest for the multi-jet problem. A representative turbulent execution on 0 < xl,/R < 1.5' required approximately 50 integration steps,, averaged 26

iterations/step, used 250,000 words of', central memory and required 850 cpu seconds in scaliar mode on, a Cyber/203 computer. The basic discretization of R2 e employed 576 triangular elements with 367 nodes and 16 degrees of freedom/node. The execution computer time was oniy doubled in executing the PNS solution over a distance seven times larger, i.e., on,0 < xl/R< 11. SUMMARY AND CONCLUSIONS A penalty finite element numericall solution algorithm, developed for three-dimensional aerodynamic flow problem classes governed by the parabolic Navier-Stoltes equations, has been examined for applicationi to analysis of strong inrteraction regions in jet' type flows. The combination of theoretical~ decisions, i:n the design of this algorithm have been documented for applliicatiion, to analysis of, a specific multiple free-jet configuration. The various decisions required to specify a turbulent flow simuiation, were presented and discussedi and examinations made to: evaluate the impact of initial condition (error) on the floww pred"uction. Comparisons were also made assuming laminar fliow,, wherein the ~ order of the predicted transverse velocity field'rmagnitude agreed well with, the PNS theoretical arguments. Itn comparison, dependent upon initial condition specifications, the turbulent flow simulations produce transverse velocity levels -ys 10-20 times larger than the order Re . The mul!ti-jieR PNS turbulent f'low, ccalculation, predicts self-generation of a higihly detailed transverse plane pair vortex structure. The robustness of' the associated mixing processes appears correlated with the available photographic and video experimental data. Tests were condlucted to verify that the character of ' the vortex structure was not deminated by initial condition assumptions. Certainly, a wide parameter range exists over which the P(^IS algorithm could be 2 7'

employed as a computational laboratory tool. For example, although not' reported herein, coarser grid solutions have been executed on the entiirc transverse plane of span 4D x 4D to compare flowfield evolution for one jiet, non- operable. The coarse-grid four jet solution agrees to within ©'u R= 2lu'9I6 with the finer grid quarter-plane solution, and the non-symmetric three-jet solution indicates the resultant mixing action is diirninished' somewhat but not destroyed. Before proceeding, with . such design, studies however, it would be of utmost importance to obtain high quality detailed experimental' data, on u i and ui ~ ~ distributions for a base configuration. The aNailability ol such data is crucial to el'imination of errors associated with initial condition assumptions. The results of' the benchmark computationai' experiments can serve to guide the laboratory experiments in diagnosing completeness (and accuracy) of the crucial data set. In this manner, computational and experimental mechanics analysis can serve each other in a synergism that will rapidly mature. V'L. REFERENCES. Schlichting, H., Boundary Layer Theory, 7th Ed., McCraw-Hill, New York, 1979. 2. Melnick, R. E., and Chow, R., "Asymptotic Theory of! Two-Dimensional Trailing Edge Flows,"' Technical Report NASA SIP-347, 1975, pp. 177-249. 3. Schetz, J. A., Injection and' Mixing in Turbulent' Flow, Progress in Astro. & Aero., V. 68, AIAA,, 11980. 4. Baker, A. J., and! Orzechowski, J. A., "An Interaction Algorithm For Three- Dimensionall Turbulent Subsonic Aerodynamic Juncture Region Flow," AIAA Jl, V. 21, No. 4, 1983, pp~ 524-533. 5. Cebeci, T., and Smith, A.M.O., Analysis of Turbulent Boundary Layers, ~ Academic Press, New York, 1974. 6. B. E., "A Reynolds Stress Model of Turbulence Hanj,3lGc, K. and Launder ca , , and its Application to Thin Shear Flows," .J. Fluid Mech., V. S2, Pt. 4, ~ Pp1 609-638'„ 1972. ~ ~. 28

7. Baker, A. J., Yu, J. C., Orzechowski, J. A., and Gatski, T. B. "Prediiction And Measurement of Incornpressible Turbulent Aerodynamic Trailing Edge Fliows," AIAA Journal, V. 20, No., 1, 1982, pp. 51-57. Dodge,, P. R. and Lieber, L. S., "A, Numerical Method'. For the Solution of' Navier-Stokes Equation For a Separated Flow," Technical Paper AIAA- 77-170,,1977. 9. Patankar,, S. V., Numerical Heat Transfer and Fluid Flow, M'cGraw- 10. Hi1llHemisphere, NY, 1'92i'0.. Briley, W. R'., and McDonald, H., "Analysis and Computation of Viscous Subsonic Primary and Secondary Flows," Technical Paper AIAA-79-1453, 1'979. 1'1. Mikhaily A. G. and Ghia K. N: "Analysis and Asymptotic Solutions of' ~ , , . Compressible Turbulent Corner Flow," Trans. ASME,, Jl Engr. Power, V. 104, 1982, pp. 57'1-579. 12. Baker, A. J., Finite Element Computational Fluid Mechanics, McGraw- 13. Hill/Hemisphere, NY, 1983. Baker, A. J., "The CMC:3'DPNS; Computer Program For Prediction of ' 14. Three-Dimensional, Subsonic, Turbulent Aerodynamic Juncture Region Flow-Yolume I - Theoretical," NASA Technical Report CR=3645, 1982~. Manhardt, P. D., "The CMC:3DPNS Computer Program For Prediction of' Three-Dimensional, Siubsonic, Turbulent Aerodynamic Juncture Region Flow - Yollume II - User's Manual," NASA Technicai', Report CR-165997, 1982. 15'- Orzechowslciy, J. A.,, "The CMC:3D'PNS Computer Program For Prediction of' Three-Dimensional, Subsonic,, Turbulent Aerodynamic Juncture Region Flow - Volume I'Il: - Programmer's Manual," NASA Report CR-165998, 1982. 16. Baker, A. J., Orzechowski, J: A., and Stungis, G. E., "Prediction o1, Secondary Vortex Flowfields lhdWced By Multiple Free-Jets Issuing, in Close Proximity," Technical Paper AIAA-83-0289; 1983. 29

Figure i. Perspective View of the Multiple Jet Configuration. Figure 2. Smoke Tracer Visualization of Mul'tipie-Jert: Configuration Flowf~ield', a) Jet' Flows Non-operational, b) Jet Flows Operating. Figure 3'. Smoke Tracer Streakline Distribution For Estimation of Velocity Vector Tangent Angle Distribution 0. Figure 4. PNS Penaity Algorithm Symmetric Quarter Plane Solution Domain,. Four Multi-Jet Configuration. Dots denote Nodal Coordinate Distribution~ of' M= 19 x 19 Discretization. Figure 5. Geometric Specification For a Two-Dimensional Slot Jet. Figure 6. PNS Penalty Algorithm Solution Distributions, Half-Plane Symmetric Subsbnic Slot Jet, u? = Xm/s, k= 50, 0 < xl/Hf< 1.0,, a)~Axial Mean Veliacity, u 1„ b) Transverse Mean Velocity, u 2, c) Twrbulent Kinetic Energy k, dlsotropic Dissipation Function e. Figure 7. PNS Penalty Algprithm Symmetric Ha1f=Plane Solution Domain. Discretization, Discrete Approximate Circul'ar Jet, M= 13 x 25, Figure 8. PNS Algorithm Transverse Plane Mean Velocity Distributions u., ,. Symmetric-HJalf'Circufar' Free Jet, Laminar Flow, u? = 30'm/s, X = 10,, a)xl/D=0:25,uT =0.001~5, b) xl/D = 01.5, uT = 0.001!4, c) xl/D = 1.0, um = 0.001.18. Figure 9. PNS Algorithm Transverse Plane Mean Velocity Distributions uij, Symmetric-Ha1f' Circular Free Jet, Turbulent Flow,, u°r =, 30 m/s, ~= 110, kj = 0.01,vA/V = 10, a) xl/D = 0,25, uT'= 0.038, b)xl/D=0:5, uT =0.0'58, c)'xl/D=1.0, uT =0.044. Figure 10. PNS Algorithm Nparfield' Transverse Plane Mean Velocity Distributions, i4jcl),, Multiple Free Jet, u y= 12 m/s, X = 8, k9 = 0.005, vS/v = 35, a)'xl/At=0.15,uT =0.0E5, b)xl/ER=0:5, u~ = 0.058, c) x l/11 ' = 1.0, uT = 0.058, d) xl/Ft = 1.5, uR = 0.067. Figure 11. Composite of PNS Algorithm Transverse Plane Velocity Distribution uQ(xl), Four Mult iple Jet Geometry, u? , = 12 rnn/s, X= 8', xl/Pt = 1.5.

F0igwre12: PNS Algorithm Symmetric Qwart~er-Plane Vel'ocityD'ist'ribut'ions u-1='12m/s,~=8, o a) Turbulent Fl+ow, k~i= 0.0025, Nd/v = 9,, b) Laminar Flow, k® = 0, vt = 0.. Figure 13. Summary of PNS Algorithm, Axial Mean Velocity, u I Maximum as Funcxiionot JetIniRial Tur~bulentKinetikEnergyLevel Assumption. Figwre 14. Summary of PNS Algorithm Transverse Plane Mean Velbcity uL, Extremum as Function of Jet Initial Turbullent Kinetic Energy Level' Assumptiion. Figure 15. Perspective Composite Graph of Multiple Jet Initial Condition Specidications, a):Jet Velocity u?('xz,xl = 0),, bYMass Fraction Yotx,,,xl = OD. Figure 16. PNS Algorithm Species Mass Fraction, Distri'butions, Yfi, Four Jet Geometry, u= 12 m1s, l= 8, aDlxl/'>7l' = .0, Ym = 1.0,, b)lx1IAt = 1.0, Ym = 1.0, cl~xll'1t=S.Oy Ym=0:63; dDlxl/R = 11.0, Ym = 0.30. Figure 17. PNS Penalty Algorithm Predicted, Distribution of Extremum Velocity Vector Tangent Angles 0i andl 0e.

X3f R / rigure 1. Perspective View of the Multiple Jet Configuratil, i micJ PNS Vortex Ftowfield Solution. giZLGC&zOz

1 r . ~ -T Um r w ~ b) Figure 2. Smoke Trncer V'isualization of Mult:ipte-Jet Configuration Flowf'iield, a):Jet F1!o•x•s ?1orn-oper:3zion3l, b)lJieT Fiows Operatinv. Smoke Col urtara iDispersed- OA' . 40

Figure 3. Smoke Tracer St'reakline Distributiom For Vector Tangent Angle Distribution 0'. ~ d ~ ~. Estirnation of Velbcity ~' ~ ~ ~

1.0 '. TRANS. COOiRD I NATE - X 3 /R 2.0 ~. O hJ G0 Figure 4. PNS Penalty Algorithm Symmetric Quarter Plane Solution Dornair, ~ Four Multi-Jet Configuratilon. Dots denote Nodal Coordinate Distribution of M= ~ 19 x 19 Discretizatiiorn. ~ N N

I Vupper I ~ ~ gf -Ujiet --~ I lower ON. X1 F'igWre S. Geometric Specification For a Two-Dimensional Slot' J+et.

2.0 r _i- = 2:0 ~ N 0 I y I _ , 0.0 40 6.0 AXIAL VELOICIT'YZI r c 1.0 I i , ~ I _: 2.0 ` N K i , ~ . . 1 12:0 -t.6 -O:e 0.0 TRANSVERSE VELOCITY u7 ~ 181 ~ 2.0 2 ~ K 0.0 0.04 0.08 TURBULENT KINETIC ENEf2IGY k (c) . ~ ~ . O!D 12 D 2+4D ~ pISS1PATCON FUNCTION E 114 Figure 6. PNS Ppnalty Aligorithm Solution Distributions, Half'.-PLarne Symmetric Subsonic Slot Jet, ui = 30 m/s, J1 = 50, 0'< xl/Hf < 1.0, a) Axial Mean Velocity w 1, b) Transverse Mean V'elocity u 2, c) Turbulent Kinetic En,ergy k, d) Isotropic Dissipation Function e.

a 0 0 . Figure 7'. PNS Penalty Algorithm Symmetric Half-Plane Solution Domain Discretization, Discrete Approximate Circular Jet„ M, = 13 x 25.,

. . . ` : 1 111, VA11,,1A: I . i' 1'; j: ! j, II' ~ / • , . t • l+ ~.- ~ .1 rl. -f- . ``% ' I t1ll1!'t 1 1 111111 1 11 • % \ tI lll'111 J1'11'11l 1 • . . \ \ \ "111Ui ! !'1j/l / I / • ". " \ \ \\ll 1 l I 1'1'llll/ / / ` , . ~). " I \ % I, 11't t s 1iJ; llI/t I r If % % `\ t111i1 i J'1!111 1 i ( . , ~~~~~ , _ _ ~~ , ~Y Figure 8. PMS Algorithm Transverse Plane Mean, Velocity Distributions u.E , Symmetric-Half Circular Free Jet, Laminar Flow, u? = 30 m[s, a= 10, a)xl/'D=0i2S,um =0.001!S, b)lxl/D = 0.5, uR l= 0.0014, c) xl/D = 1.0, urp '= 0.0018. . • \ ~ \ ~ .~ 1~ : %~% ~ 1 I 1, ,~ j . I Ili I~ ,( ~ .~ .' I

V 1: i:1~1 : i i 1 i i 1'J i1 :; , ' X I & \ 1111 111 11 1 i'k,"i i : / % ~ ~ 1il1 ! J1 Wr1~'~ ~ ' ~~ ' `~ ;AI 1~ I 111! ~ ~ W" _. ~ ~ ••- - ..-- a) . , . . . .. , , ; , ~ • ` ` " x 1 I 11'111ii1 l1I1Jl1! I, : 'r . .. . , , , . ` >>11411111'111llill 1 ~ 11 ~. , r- - ~- _+. . D--~I ~ Figure 9. PNS Algorithm Transwerse Plane Mean Yslocity Distributions u Symmetric-Half Circufar Free Jet,, Turbulent. Flow, uli =, 30 mfs, X= 10, lC~= 0.01!, Va/V = 10; a) x l/'D = 0.25, uT = 0.03'g, b) x l/'D = 0.5, uT = 0.058, ~' c) xl/D = 1.0, uT'= 0.0!4'4.

~ ~ ~ l l 1I1IJJ11 1 l I I 1 [14 111111 / / . . . .... . . ; . ., , 4. O'~4 '_ -400 X3,/R'. I , , , , , , ~X~ . . ._= .. ` 1j.,/ 40 w. ~.~ ~~ , . . , ; i ~ - ~ ~ r b)' Figure 10. PNS Algorithm Nearfield Transverse Plane Mean Velocity Distributions, uqxl), Multiple Free Jet, u? = 12 m;/s, X= 8, kP= 0.005, vS/v = 35, alxl,/R=0.15,uT=0.085, b)~xl/R = 0.5, uT'= 0.05'8,. c) x l/R = 1.0, uT '= 0.058', dDxl/R= 1.5, uT 1 =0.067. '- . : : ,.: ' , :,~i~'j-, . ,

1 I I /,,.... , . I . I dill I I . !! 11 1'' ' ` i f_ l . -!'~~'~ _ ~1 % 6 r l . R3/R 1 i . 1 1 / ....ttl! ! ! d). x -Ii Figure 10. Ph1S Algorithm Nearfield Transverse Plane Mean Velbcity Distributions, ujxl), Multiple Free Jet, u~= 12 11 X= 8, ej = 0:005, vdLv= 35, a) xj/R = 0.15, uR = 0.085, b) x l/AC = 0.5, wT = 0.058, c) x l/R = 1.0, ui m = 0.05'8, d) xl/R' = 11.5, , wm = 0.067.

I 1 1.1111- - '' I rI I ...;.4 1 1 1 I t . % t 1i 1 r...,t 11 1( / i , , / / ! I,l1u, . , I I I I I ' j 1 111 ti,•., 1 l f, I I I i" ..... e 1' f 1 I ! • . . .......~t~ t t r 1 1 J ! 1 ...... ~ ~ . 1 I 1 1 j 4 •4 11 1 1 , . 1 I 1 I . •, tf ll'I / , • 1 ! I ii t I 1 t,.,,rs r 1 1 1 I tr.._,,„ ~ ~ 1 1 t t' t r....• .~ . Figure 11. Composite of ' PNS Algorithm Transverse Pllane Velocity Distribution u,(xl), Four Multiple Jet Geometry, ul = 12 mJs, A= 8„ xli/R = 1.5.

1 11 1 i i1/,rr'j/ 1/ ' 1 4 1 4 4 4 0.,4p 4f f f I 0 - fiJl I r ~ . I / I' I =, ~ l I ~I % f ~ ~. % d . %_ - f c.--''•-_ ~ - i IX2'R. t' 1 t t rrrrr•rr' r . - -{- X3/R a) t' t~~ r' r. U r r.r~~ r. r r' I~ .. .~~ p .. I t t t~. r t r r~ r~o.iI i. .. . . .. t~~ l' t t rtt.rr r r i . . . 1 /, f N` tllti 4 i . . . . t iiPlL'fyi i . . ~ Ifff // I 1 t0`S 0 0 . i Figure 12. PMS Algoriithrrn Symmetric Quarter-Plane Velocity Distributions ut url= 12m/s,a=8, : a) TwrbwlenA Flow„ k9= 0:0D2'S, va'ti' = 9, b) Laminar Flovw, h0 = 0, vt = 0-

- o- -o---o---t~-----~ } F 0:75 U 0 J w > J 0.50 Q x Q ~ .25' w ~ 0 0 0.5 1.0 AXIAL COORDINATE'- L / 11 . 1.5 Figure 13'. Summary of PNS Algorithm Axial Mean Velocity u 1 Maximum as Function ol' Jet Initial Turbulent Kinetic Energy Level Assumption.

0.5 160 AXIAL COORDINATE - x, / R 1.5 Figure 14- Summary of PNS AlgoriR'hm Transverse Plane Mean Veiiacity u Extremum,as Funcfion of Jei lnitial Turbulent Kinetic Energy Level Assumptior~:

v v N U N ^' fp N (~~ . ~ < .~ .r ~ ~ 1 I ~ 0= ~ ^ _ x .N. M x .~a + n ~y „xr II Q ~r v t,

Y Y Figure 16. PNS Algorithm Species Mass Fraction Distributions, yh, Four Jet. Geornetry, ul = 12 m/s, a= Si a)xli/R=.0, Ym= 1.0, b) x 1 LR = 1.0, Y'R' _, 1.0, c)xl/R=S:0,, Y'm=0.63, d)' x l/Et' = 11.0, Y'm = 0.3i0..

Y 2 Figure 16. ~,1~S Algorithm Species Mass Fraction Distributions, Yh, Four Jet Geornetry, ul = 12 m/s, = 8, a,llxl/'R = .0,, Ym = 11.0, bl xl/R = 1.0, Ym = 1.0, c) z1CR = 5.0, Y'm = 0.63, d) x 1 /R = 11.0, Ym = 0.30.

e a B a a AX IAL DISPLACEMENT XI/'R Figure 17. PNS Penalty Algorithrra Predictied! Velocity Vector Tangent Angles 01 and o!e. Distribution of Extremiwrn