https://ntrs.nasa.gov/search.jsp?R=19910011768 2019-04-12T12:55:40+00:00Z N 9 l - 2 PRESSiTRE fiGORITjBM FOR EWIBTIG PLOW CmIKATIONS WITB m&P DF mTBOD M. S. Anand Allison Gas Turbine Indianapolis, Indiana 46206 S.B. Pope Cornell Univer§ity Ithaca, New York 14853 H.C. Mongia Allison Gas Turbine Indianapolis, Indiana 46206 % $, 6 !t r ABSTRACT An algorithm to determine the mean pressure field for elliptic flow calculations with the probability density function (pdf) method is developed and applied. The pdf method is a most promising approach for the computations of turbulent reacting flows. Previous computations of elliptic flows with the method have been performed in conjunction with conventional finite-volume based calculations that provided the mean pressure field. The algorithm developed and described here permits the mean pressure field to be determined within the pdf calculations. The pdf method incorporating the pressure algorithm is applied to the flow past a backward-facing step. The results are in good agreement with data for the reattachment length, mean velocities, and turbulence quantities including triple correlations. INTRODUCTION The understanding and modeling of turbulent reacting flows is of great import- ance due to the occurrence of such flows in most practical combustion devices. The probability density function (pdf) method offers several advantages over conven- tional turbulence models for the computations of turbulent reacting flows (ref. 1). In this method a modeled transport equation for the joint pdf of velocities and scalars is solved using the Monte Carlo technique. The primary advantages of the pdf method are that, in the pdf transport equation, terms representing convection and reaction appear in closed form and need not be modeled. Hence, conventional turbulence models such as the k-E or second-moment closure models which are based on gradient-diffusion assumptions are entirely avoided. Also, any arbitrarily com- plex finite-rate reactions that are beyond the reach of conventional turbulent re- action models can be treated without approximation or restrictive assumptions. Until recently, most of the flows studied using the pdf method have been two- dimensional (2-D parabolic or 1-D time-dependent) in nature. These studies (refs. 2-7 to name a few) have demonstrated the potential advantages of the pdf method es- pecially for reacting flows. There are no theoretical limitations for applying the pdf method for complex (2-D or 3-D) elliptic flows. However, a suitable algorithm for extracting the mean pressure field within the pdf solution algorithm is needed. A pdf method for elliptic flows was first demonstrated by Anand, et al. (ref. 8) for a 2-D flow. The method was applied to 3-D in-cylinder flows by Haworth and El Tahry (Refs. 9, 10). In this method, the Monte Carlo (MC) calculations for the pdf are combined with conventional finite-volume (FV) calculations (ref. 11) for mean fields. The mean pressure field and the turbulent time scale are supplied to the MC calculations by the FV calculations. In turn, the MC calculatio~~ssu pply to the FV calculations the Reynolds stresses extracted from the pdf solution, thereby avoiding the use of conventional turbulence models. This combined algorithm ex- ploits the advantages of both methods to yield an economical algorithm for pdf eal- culations of elliptic flows. This method was intended to demonstrate the feasi- bility of pdf calculations for such flows, However, for more complex flows, especially variable density and reacting flows with fast or finite-rate multistep chemistry, the coupling between the two methods becomes quite complex and it becomes more difficult to fully exploit the advantages of the pdf method. Therefore, it is desirable to perform pdf calculations inde- pendently, without recourse to FV calculations. This paper describes an algorithm by which pdf calculations are performed independently and economically (in terms of computer resources required) for complex elliptic flows. While the present work focusses on the determination of the mean pressure field within the pdf calcula- tions, parallel work is in progress (ref. 12) that deals with the determination of the turbulent time scale. For the present work, the turbulent time scale is sup- plied to the calculations as described later in this paper. The relevant equations for determining the mean pressure are presented, the solution strategy and the numerical scheme are developed and suitably implemented in the pdf calculations. The method is applied to a sample elliptic flow, namely the recirculating flow be- hind a backward-facing step, and compared against data of reference 13. THE PDF METHOD A very brief description of the pdf method is presented here. The reader is referred to reference 1 and other studies listed in the introduction for more de- tails. The joint pdf, f(V, $; x, t), is the probability density of the simultaneous event Y(x, t) = Y and @(x t) = 9, where Y is the velocity vector, $ is a set of scalars, and y and 9 are independent variables in the velocity-scalar space. Starting from the usual conservation equations for mass (continuity), momentum and the scalar quantities, the transport equation for the joint pdf can be derived as described in reference 1. In this equation, the terms involving convection, reac- tion, body forces and the mean pressure gradient (including the variable density effects in those terms) appear in closed form and need not be modeled. Hence, the use of conventional turbulence and reaction models is avoided. However, the terms in the pdf equation representing the effects of viscous dissipation, fluctuating pressure gradient and molecular diffusion need to be modeled. A Lagrangian view- point is adopted in modeling and solving the joint pdf equation. The modeled pdf transport equation is solved by the Monte Carlo technique. In the Monte Carlo solution technique, notional particles representing fluid particles are distributed throughout the solution domain. Each particle is attri- buted with values for the velocity components, spatial coordinates, and scalar quantities. Each of these values evolves according to its corresponding Lagrangian evolution equation which incorporates modeled terms where needed. The solution of these evolution equations constitutes the solution of the pdf transport equation r e f 1). Means of any functions of the independent variables are determined by appropriately summing over the values of all the particles in small spatial sub- volumes and fitting splines over the whole domain to these sums. In the context of the scope of the present study, the scalar quantities are ig- nored in the presentation except to note that the density p is a function of the scalar variables. During an interval of time dt, the change in particle position x>k is given by where 4$t is the particle velocity. With the use of the simplified Langevin model (refs. 2, 3) the increment in y* in the interval dt is given by: where angled brackets indicate mean quantities, tilde denotes density-weighted (Favre) means, <P>(x*) is the mean pressure, p$t is the density of the particle, ZU( x>k) is the Eulerian mean velocity, k is the turbulent kinetic energy, w(x*) is the turbulent frequency, Co is a universal constant and W(t) is an isotropic Wiener process. The last two terms in Eq 2 jointly represent the effects of viscous dissipation and fluctuating pressure gradient. The value Co = 2.1 was determined b~ Anand and Pope (ref. 2) and has been used in the present study. The quantities U and k are easily extracted from the pdf solution at time t. In previous studies, -. the mean pressure gradient has been determined through boundary-layer assumptions or has been supplied by the accompanying finite-volume calculations. Its determina- tion within the pdf calculations for elliptic flows is the topic of this paper. The determination of w is also discussed. PRESSURE ALGORITHM Overview The pressure algorithm is constructed for steady-state, constant- or variable- density, high-Reynolds-number flows. The equation for mean pressure is the Poisson equation (Eq. 4) derived from the mean momentum equation (Eq. 3): Since the right-hand side of the equation can be evaluated from the pdf solution, Eq. 4 can be solved for the mean pressure cP> with appropriate boundary conditions. Starting from arbitrary initial conditions, if an iterative or psuedo-time marching algorithm is used to reach the steady-state for variable-density flows, then, be- cause the transient terms are absent from Eq. 4, the mean pressure given by Eq. 4 is not correct until the steady-state is achieved. This is true even for constant- density cases unless the initial velocity field satisfies the continuity equation and Eq. 4 is solved in a coupled manner with the mean momentum equation which is not solved directly in the pdf method. A consequence is that the mean continuity equation will not be satisfied. Given a field that does not satisfy Eq. 5, a velocity cor- rection All can be obtained by + By requiring that {<pUi> <p>AUi) be divergence free, we obtain a Poisson equation for @: These observations suggest the following basic algorithm: every so-many steps 1. solve Eq. 7 for (3 2. add the velocity correction obtained from Eq. 6 3. solve Eq. 4 for <P> The spirit of the algorithm is similar to that of the SIMPLER algorithm used in conventional FV calculations (ref. 11). However, the solution strategy and the im- plementation in the pdf calculations are markedly different. Boundary Conditions for tP> The mean pressure <P> is uniquely determined (to within an arbitrary constant) by - Eqs. 3 and 4. That is, given Ri* - <pUiUj>, <P> can be determined from the Poisson equation (Eq. 4) wit& the Neumann boundary conditions provided by Eq. 3: where is the outward pointing normal and n is a coordinate in that direction. The numerical solution of a Poisson equation with Neumann conditions everywhere is somewhat perilous. Unless suitably treated, the coefficient matrix is singular; and a solution exists only if the boundary conditions are precisely consistent with the source. To avoid these difficulties, and to make the same boundary conditions applicable to @ (see below), a Dirichlet condition is used at the outflow bound- ary. There the flow is approximately parallel (to the x-axis) and fully developed, and so the lateral momentum equation becomes (approximately) where v" denotes the lateral velocity fluctuation. Hence the Dirichlet condition used at the outflow boundary is The arbitrary constant Po in the pressure solution is chosen each time so as to match a known pressure at any given poiiit in the solution domain. Boundary Conditions for @ At walls, planes of symmetry, and inflow boundaries, the normal component of the velocity correction should be zero. Hence the Neumann condition is appropriate. At the outflow boundary, specifying the Dirichlet condition, Q = 0, results in there being no tangential velocity correction, but allows a correction to the normal velocity. This is necessary if the mass flow rate in differs from the mass flow rate out. NUMERICAL SOLUTION In the context that the Monte Carlo calculations are grid-free and the mean fields are represented by spline-fits, it is advantageous to solve the Poisson equations (Eqs. 4 and 7) in a way that is completely consistent with the spline representation of the mean fields. The following discussion is in the context of 2-D flows, but its extension to 3-D flows is straightforward. Let g(x,y) be any function within the solution domain. Then its spline rep- resentation is where aa(x)(a = 1, M) are the basis functions in x, bn(y)(13 = 1, N) are those in y, and G~~ are the spline coefficients. Let Qyfi be the spline R?: coefficients of -<pUi>, and let be those of -<pU.U.>. With l.7 the substitution of the appropriate spline representatlons, the Poisson equations (Eq. 7 and Eq. 4) become and a13 II C C P (a "b n + a bR)= a a a 13 where pfaind paB are the spline coefficients of (4 and <P>, and primes denote differentiation. Equations 13 and 14, along with the appropriate boundary conditions can each be represented by A.s = h, (15) - where s is the solution (spline coeffici-ents of <P> or @) and A is the co- efficient matrix. Since both <P> and Q satisfy the same types of equations and boundary conditions, their coefficient matrices are identical. Further, since is a constant matrix (independent of R(t..8 )a nd Qi), its J J. decomposition need be performed once only, and then Eq. can be solved many times by back substi- tution. The source vector h is different for <P> and @ and changes from iteration to iteration as the steady-state is approached. The structure of matrix A is shown in Figure 1, with x's indicating the only entries in & when cubic B-splines (refs. 1, 14) are used. A special band solver, which exploits the structure of the matrix & was developed to minimize the storage and computational requirements. The solver uses EU decomposition and Gaussian elimination with partial pivoting. As pointed out earlier, the LU decomposition is performed only once and the Poisson equation is solved repeatedly by back sub- stitution. The work to perform back substitution is of order m(4M + 5) operations: just 34,000 for N = M = 20 used typically (M = N = 21 in the present study). The majority of the work is done in calculating b. This involves forming splines for <pUi> and <pUiUj> with appropriate boundary conditions. These splines need to be very accurate especially since their derivatives and second derivatives are to be evaluated (see Eqs. 4 and 7). Also, the boundary conditions used in forming these splines have a significant influence on the solution of the Poisson equations. Both the issues need careful consideration. In general, the splines in the pdf calculations are formed using the cross-vali- dation procedure (refs. 1, 14), by which splines are fit to two independent sets of sums (or samples) from each of the spatial bins such that the resulting spline is the best fit for both sets of samples. This procedure considerably reduces the er- ror in the spline fit with respect to the actual function, compared to the error with a simple least-square spline fit with a single set of samples. For spline-fits .... 1 XXXXXX ....... XX.XX XX........XXXXXX............................................................................. . XX..........XXXX....... XXXX.XX.... .XX ...XX X X.. XXXXXX... XXXX..X... XXX ....... XX...... ..XXXX....... XXXXXX ..... .X.XXX.XX... . ...XXX.X XXX . . .......XXXXXXX.....XX..X..XXXX....X.X.X.X..XX.....X..XX..X.XX......X.XXXX..XX.......X..XXX.X. .......X....X.X................................................................................................................................................................................................... ...... .................................................................. .................................................................. ....... ... ...... .. ... ..................................................... .. .................3....X....X..... .XX ... XXX XXX X ..X ...X....X......X X. . XXX.XXX ..X...X......XXX.X.XXX.XXX. ..X............................................................................................ . ............................X...XX. XXX XX X. .......X..X X..XXX XX..X......X...XX.XXXXX..X..................................................................................... ........................,...........X...X...X . XX XXXX ..X...X...X....XX . ..XX .XXXX....X.....XX..XX..XX................................................................. ........ ...............................................................................................................................................................................................................................................X....... .......XX...... .......XXX ......... .. XXXX...........X.XX............ XX.......... .. XXX... ........XXXX...X...... .XX.XX..XX .......X. .XXX..X. . . . .X XXXX XX.........X XXX.XX........X.XX. XX........ .X. ..X .XXX..... . . ..XX XX...XX XX....X..X.XXX. ..XXX XXX.X.X... X..XX.XX...X....... . X.XX.XX.X....Xx.......x .X.XX..X.X.XX.x........X.XX..X.X.XX..........X.X..X.X..XX........x.X.XX.........x..........x..XXX..X...................X..XX......xx...........X...X...xx.....................xx.................................... ................. .......................................... .................... ...... ...... .... .. ........ ...............................XX X....X.xx.......x xx........~ ..............................................................................................................................X....XxX .. XxX.. xX x..x.x...X.. XXX..XXXx.XX x....x.....X....XXX..XXXx.XXx....x..................... ........................................................................ X.XX XXX.......x-xxxx.x......x.xxxx.x...... ........................................................................... .X.XX XX.x......x.xxxx,x.....x..xxxx.x. . .................................................................... ......Xx xxxx..-~...x.xxxxx.......x.xxxx x ................................... XXX....xxx....xxx.... .......................................................................... X XXXx.X......x.xXxX.X.....x..xXxX.X..... . ............................................................................ .x xxXx.x....-x..xxxx..x....x.x.xxx..x . ....................................... Xxx....xxx....xxx Xxx....xxx....xxx - e A M = 7, N = 8, (56 x 56). - TE89-4904 Figure 1. Structure of the coefficient matrix A (MN x MN) for sample values of M and N. There are 5 bands, each with a width of 5 and separated by a width of M -5. The total bandwidth is 4M + 5. Table 1: Boundary conditions for splines needed to evaluate h Quantity Inflow boundary Outflow bourldary Wall Wall 4 normal to x normal to x normal to x normal to y 0'" R >" qx = qx = 0 q = 0 - q = 0 q = 0 - q = 0 >*" * subscripts x and y denote differentiation with respect to x and y, respectively - indicates that no boundary condition is specified in two-dimensions, four independent samples are needed. However, even with cross- validation, although the first derivatives are acceptable, the second derivatives evaluated from the cubic spline fits can have errors significant enough to adversely affect the solution of Eq. 4 and hence the convergence of the whole pressure algorithm. Hence for Eq. 4, splines are fit for Rij as well as to its appropriate first derivative from which the required second derivative is evaluated by a single differentiation. From the same raw sums for Rij four different cross-validated splines are formed for Rij using four different sets of values for M and N (the number of basis functions). The first derivatives evaluated from each of the four splines form the samples, although not entirely independent, for the spline-fit for the first derivative. This procedure significantly reduces the errors in the second derivatives evaluated and results in a stable solution of Eq. 4. The four splines were formed in the present study with values of M (=N) equal to 17, 19, 21 and 25. It is also necessary to specify a consistent set of boundary conditions for the spline fits. The quantities splined for solving Eqs. 4 and 7 and the boundary con- ditions used are listed in Table 1. Note that the density has been assumed con- stant. The table should be applicable to variable-density flows also, with the conventional means replaced by Favre means. Additionally, boundary conditions for <p> will be needed--a zero normal-gradient condition at all boundaries is sug- gested. The development and the implementation of the algorithm to solve for the pres- sure field within the Monte Carlo solution algorithm were evolutionary processes. Primary concerns were to obtain an accurate pressure solution and also keep the computational times low. The algorithm can be summarized as follows: I. form the matrix and perform its decomposition 2. initially set <P> = 0 everywhere 3. march the pdf calculations by a step in psuedo-time 4. perform velocity correction; i.e. o form splines of <p> and <pUi> o solve Eq. 7 by back substitution for o correct particle velocities by Eq. 6 5. on the first step and subsequently, every S steps update pressure; i.e. solve for the pressure correction A<P> from Eq. 4 a 2 ~ 2 2 old V A<P> = ii - V <P> , ax. ax 1 j and reset the mean pressure by where r is a relaxation parameter. 6. if steady-state has not been reached, go to step 3 7. stop Two parameters S and r related to pressure updates, are introduced. For the case studied, neither parameter influenced the results significantly. The pressure updates are considerably more expensive than velocity corrections. The pressure updates need not be performed frequently since the velocity corrections, which cause a similar effect as the pressure corrections, are made every step. Also, the velocity corrections are quite small from step to step (especially if continuity is satisfied initially) so that the pressure updates are not necessary every step. It is better to update the pressure after the velocity field has evolved appreciably, say every 10 to 20 steps. For the present study, values between 0.5 to 0.8 were used for the relaxation factor without any apparent effect on the results or the number of steps needed for convergence. Although convergence is not difficult to achieve, it is difficult to monitor. Several parameters were considered as indicators of convergence. But, due to the stochastic nature of the calculations and due to the fact that the solution evolves slowly from step to step, a robust and reliable parameter has not been found. In the present study; the convergence was inferred by comparing the results at the 200th step with those at the 220th and the 250th step and concluding that the re- sults were unchanged within statistical errors. TEST CASE AND INITIAL CONDITIONS The pressure algorithm and the pdf calculations were tested for the flow past a backward facing step with data from reference 13. The schematic of the flow is shown in Figure 2. The step height is H and the centerline mean axial velocity in the channel is Uref. The fluid used in the experiments was water, and the Reynolds number based on H and Uref was approximately 12,000. The inlet velocity pdf was prescribed to be Gaussian for each of the velocity components with the means and variances being those of fully developed turbulent Figure 2. Schematic of the flow considered. Backward-facing step studied experimentally by Pronchick and Kline (1983). flow in the inlet channel as calculated by the k-E model. The inlet turbulent kinetic energS is distributed into the three components according to k = <ut2>= >. 2<vf2>= 2<wt The turbulence levels a t the inlet have little effect on the final results since most of the shear and turbulence production occur in the re- circulation zone within the solution domain. The initial pdf in the solution domain was set equal to the inlet pdf and the initial mean pressure was zero everywhere. The solutio~id omain extended from the step (X/H = 0) to x/H = 15 and between the two walls in the y-direction. The turbulent frequency (w = ~ / k ,w here E is the dissipation rate of k) is calculated at each step from the current value of k by where the turbulent length scale % is supplied from the k-E solution of ref. 8 (see Eq. 19) and is held constant throughout the calculations: The choice of supplying Q was preferred over the choices of supplying w or E from the k-E solution since the profiles of Q were nearly uniform across the flow (approximately 0.2 W3) except near the walls where the values decrease. This way, the w used in the calculations (Eq. 18) would be more consistent with the current values of k. As mentioned earlier, work i s in progress (ref. 12) that deals with the determination of w within the pdf calculations. RESULTS AND DISCUSSION The number of particles used in the simulation was 60,000. The number of spatial bins used was 1500 (50 in x X 30 in y). The computations were performed for 250 steps with pressure updates every 20 steps and the pressure relaxation fac- tor of 0.8. The time interval for each step was 0.3 H/Uref. This time interval was less than one-tenth of the turbulent time scale (inverse of w) over most of the flow domain except for small regions near the wall where it was between one- fifth and one-tenth of the turbulent time scale. The calculations reached a steady- state around the 200th step. As noted earlier, a more precise way of determining convergence needs to be identified. The final results shown here (except for the mean pressure) are caIcufated from spline-fits to sums formed over 50 steps after steady-state has been reached. This contributes significantly to reducing the statistical error in the splined results especially for the triple and higher correlations which are more prone to such error. This procedure, as an alternative to increasing the number of particles in the simulation, allows the computer storage requirement to be kept at a modest level. The computed reattachment length of 6.65H is in good agreement with the experi- mentally observed reattachment length (xRE) of 6.75H. The results for mean velocities and various turbulent quantities are shown in Figures 3 - 11. The re- sults are shown for different axial stations normalized with respect to the experi- mental reattachment length: Overall, the results are in excellent agreement with data both in magnitude and in qualitative trends. It is noteworthy that the turbulence quantities, especially the third moments (which cannot be calculated with the k-E model, and are modeled in second-moment closures) are in good agreement with data. The slight disagreement between some of the computed turbulent quantities and the data in the vicinity of y/H=2.0 could be due to the inadequacy of the prescription of the turbulent frequ- ency in the present study. It is expected that the use of the stochastic dissipa- tion model (ref. 12) will remedy this deficiency. -.5 .0 1.0 .0 1.0 .0 1.0 .0 1.0 TE90-960 <U>w ref Figure 3. Computed mean streamwise velocity profiles (lines) compared against data (symbols).