- Email: [email protected]
Direct numerical simulation of turbulent spots
Computers & Fluids 30 (2001) 533±541
www.elsevier.com/locate/comp¯uid
Direct numerical simulation of turbulent spots Arup Das 1, Joseph Mathew * Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560 012, India Accepted 10 November 2000
Abstract A group of high-order ®nite-dierence schemes for incompressible ¯ow was implemented to simulate the evolution of turbulent spots in channel ¯ows. The long-time accuracy of these schemes was tested by comparing the evolution of small disturbances to a plane channel ¯ow against the growth rate predicted by linear theory. When the perturbation is the unstable eigenfunction at a Reynolds number of 7500, the solution grows only if there are a comparatively large number of (equispaced) grid points across the channel. Fifth-order upwind biasing of convection terms is found to be worse than second-order central dierencing. But, for a decaying mode at a Reynolds number of 1000, about a fourth of the points suce to obtain the correct decay rate. We show that this is due to the comparatively high gradients in the unstable eigenfunction near the walls. So, high-wave-number dissipation of the high-order upwind biasing degrades the solution especially. But for a well-resolved calculation, the weak dissipation does not degrade solutions even over the very long times O 100 computed in these tests. Some new solutions of spot evolution in Couette ¯ows with pressure gradients are presented. The approach to self-similarity at long times can be seen readily in contour plots. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Turbulent spots; Direct numerical simulations; Stability
1. Introduction An isolated region of turbulence that is convected in an otherwise laminar ¯ow has been termed a spot. Spots have been observed in wall-bounded shear ¯ows such as boundary layers, (free surface) ®lm ¯ows and channel ¯ows. (see Ref. [9], for a review). In these ¯ows, transition begins in the region where spots appear. Spots grow as they travel downstream, and coalesce to give rise *
Corresponding author. Tel.: +91-80-309-2480; fax: +91-80-360-0683/+91-80-360-0134. E-mail address: [email protected] (J. Mathew). 1 Present address: Department of Theoretical and Applied Mechanics, University of Illinois at Urbana-Champaign, Urbana, IL 61820 USA. 0045-7930/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 0 0 4 - 4
534
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
to the fully turbulent region. Although several properties of spots have been documented, it is still not known whether sucient organization exists within turbulent spots to construct a simple model of their dynamics. It would be useful to explain the dierences in shapes and spreading of spots, observed in dierent background ¯ows, in a uni®ed way. The search for simple organization within spots can be facilitated by conducting direct numerical simulations (DNS), which can provide complete ®eld data ± especially vorticity ®elds that are not easily obtained in experiments. Often, resolution and geometry requirements have limited DNS to studies in homogeneous turbulence or other simple ¯ows. But, as a naturally isolated region of turbulence, the spot presents a particularly suitable ¯ow for DNS studies. We have applied variants of a method presented by Rai and Moin [8] to simulate spot evolution in channel ¯ows driven by shear and by pressure gradients. The method involves high-order upwind-biased ®nite dierences which gives high-order dissipation. Good results were demonstrated in the original presentation [8], but the method has since then fallen out of favour because its high-wave-number dissipation is unsuitable for Large Eddy Simulations (LES) [6]. Rai and Moin [8] have already noted the bene®cial characteristics. We expect some speci®c bene®ts. As a ®nite-dierence method, we expect it to be extended easily to calculations in practical geometries of small curvature such as wing±body corners. As a practical matter, high-order dissipation which is eective only at the smallest scales can be helpful in stabilizing DNS. For economy of overall eort, it is common to begin DNS studies with grids on which the solutions are slightly underresolved, examine the solutions, and then employ ®ner grids. The merit of high-order dissipation is that it stabilizes the calculations to provide some preliminary results on small grids and then provides better solutions on ®ner grids. In this paper we focus attention on a development experience. To understand the long-time ®delity of the method, knowing that it has weak dissipation, we applied the method to compute the evolution of disturbances to an unstable Poiseuille ¯ow. It turns out that the unstable eigenfunction has a pathological, rich structure close to the channel wall which requires very ®ne grids. The evolution on coarse grids showed decay of perturbation energy which can be mistakenly attributed to numerical dissipation. However, the numerical method is able to compute the evolution of decaying modes for very long times accurately, even on such grids, because the eigenfunction is smooth. We conclude with some aspects of the solutions to spot simulations. 2. Numerical method DNS have often employed spectral methods, but increasingly, it is becoming clear that ®nitedierence and ®nite volume methods can be used successfully when the algorithms are constructed with sucient care. Thus Rai and Moin [8] presented a general, ®nite-dierence scheme using high-order upwind-biased dierence formulae. High-order formulae are used to compete with spectral methods and keep phase errors small over a larger wave number range. Upwind biasing, rather than upwinding, permits compact stencils while retaining high-order dissipation. DNS implies an attempt to resolve all signi®cant scales of motion, down to the scales where dissipation occurs. So the weak, high-order dissipation should be acceptable as demonstrated in Ref. [8]. But, for LES, as very little energy dissipation should occur at the smallest resolved scales, this method is unsuitable [6].
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
535
We integrated the Navier±Stokes equations by taking the ¯ow to be at constant density and viscosity. The numerical method used here involves only minor variations to that presented in Ref. [8]. As a splitting method, the ®rst step is an explicit integration of the momentum equations to ®nd an intermediate velocity ®eld, with a given pressure ®eld. The second, corrector step ®nds the pressure distribution and the correction to the intermediate velocity ®eld that makes the resultant ®eld divergence free. The corrector step requires the solution of the coupled system of equations for pressure. In the present problems involving channel ¯ows, the solutions are periodic in the streamwise and spanwise directions. So we use Fourier expansions to reduce the problem and then use a direct, banded-matrix solver. The discretization equations are not repeated here as Ref. [2] is widely accessible. While developing the codes, the ®rst-order Euler and third-order Runge±Kutta methods were used for time stepping. Convection terms were represented with second-order central, third-order upwind-biased and ®fth-order upwind-biased dierence formulae while second-order and sixthorder central dierences were used for viscous terms. The results of spot evolution presented here were obtained using Euler time stepping, ®fth-order upwind biasing of convection terms and second-order central dierencing of viscous terms. The numerical divergence and gradient operations on the velocity and pressure ®elds that are done in the corrector step were done with fourthorder central dierences. The main dierences between the present implementations and that in Ref. [8] are: (1) we use uniform grids, and (2) we use linear interpolation to ®nd velocity crossproduct terms rather than the fourth-order cubic interpolation. However, for the subject of this paper, the use of uniform grids has had the more signi®cant eect than the order of the scheme. All quantities are non-dimensionalized using the channel halfwidth as the length scale. The velocity scale is the velocity of the moving plate in all Couette ¯ows, and the centerline velocity for Poiseuille ¯ow. The Reynolds number Re is always constructed from these two scales and the ¯uid viscosity. 2.1. Linear stability test for long-time accuracy The initial tests of the codes for steady state solutions to the 2-D and 3-D lid-driven cavity at low Reynolds numbers were satisfactory and presented no surprises. Solutions to 3-D cavity ¯ow are few, but ours compared closely with those of Deshpande et al. [2]. The comparisons are in Ref. [1]. A more appropriate test of a numerical method which has been designed to compute turbulent ¯ows accurately must evaluate the ®delity of the computations against a known unsteady ¯ow. The growth of prescribed small amplitude disturbances to a Poiseuille ¯ow is a natural choice as solutions can be found accurately by solving Orr±Sommerfeld eigenvalue problems. Moreover, though obtained numerically, these solutions are well established. As Rai and Moin [8] had reported good comparisons with this test problem as well, this test was also expected to be routine. The temporal stability of the plane Poiseuille ¯ow has been studied by Orszag [7]. All modes are damped when the Reynolds number is below 5772. We considered the mode of wave number unity and Re 7500 for which the complex frequency x 0:24989154 i0:00223498. The eigenfunction was obtained using a spectral code and superposed on the parabolic Poiseuille ¯ow. The perturbation amplitude is 0.0001. Fig. 1(a) shows the growth of the energy of the perturbation velocity ®eld on three grids: 32 65, 32 129 and 32 257. In each case, there are 32 equispaced grid points in the
536
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
streamwise direction and an increasingly ®ner equispaced grid normal to the walls. This was computed with ®fth-order upwind-biased convection terms and second-order central dierenced viscous terms. The result of linear theory, which is the exponential growth of the perturbation is
Fig. 1. Evolution of perturbation energy with time. Graphs (a) and (b) are of growing mode at Re 7500; graph (c) is of the decaying mode at Re 1000. The scheme used is ®fth-order upwind-biased convection and second-order central dierenced viscous terms for results in graphs (a) and (c); all terms are with second-order central dierences in (b). The solid line in each graph is the reference (Orr±Sommerfeld) solution. Curves in graphs (a) and (b) are from the dierent grids: (- - -) 32 65; (± ± ±) 32 129; ( ) 32 257. The dashed curve in (c) is from a 32 65 grid.
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
537
the straight line. It is clear that even 257 points are not sucient to simulate time-accurate evolution as even the slopes at long times are noticeably dierent. At lower resolutions, the solutions are decaying. Also, the initial evolution is quite dierent from that expected from linear theory. The corresponding curves in Fig. 1(b) are from a calculation that used second-order central dierencing of all terms and the long-time results are better. The dierences at short times remain. Fig. 1(c) is from a similar calculation, but the initial perturbation is the damped eigenfunction of wave number unity at Re 1000. Now, even the initial dierences are minor, and a 32 65 grid and the ®fth-order upwind-biased scheme is able to capture the correct decay rate. It is tempting to dismiss the upwind-biased scheme as being too dissipative by comparing Fig. 1(a) and (b), whereas Fig. 1(c) shows that the method can remain accurate over long times. Simulating the growing mode turns out to be a specially dicult problem. There are two aspects to the diculties, and both are related to the solution, i.e., the structure of the eigenfunctions. First, the growing mode is not adequately resolved until a comparatively large number of grid points are used. Secondly, the ®ne structure exaggerates the high-wave-number dissipation of the upwind-biased method compared to the dissipation-free central-dierence method. The eigenfunction has a very rich structure near the walls and unless those variations are resolved, it is not possible to simulate its behaviour accurately. Fig. 2 shows the structure of the eigenfunctions W y of the growing (Re 7500) and decaying (Re 1000) modes close to a wall. The third derivative of W y, which appears as a viscous term, has been plotted to see the sharper gradients of the growing mode. Signi®cantly more points are required in a uniform-grid computation to resolve these gradients. Rai and Moin [8] were able to get good solutions on a 32 65 grid by clustering wall-normal grid points. The spacing increased in geometric progression from the wall. It turns out that the choice of clustering is also important: even when Chebyshev points were used, Malik et al. [5] showed that about 256 points across the channel are needed to reproduce an acceptable growth rate. They used a second-order ®nite-dierence scheme. Both schemes provide clustering near the wall, but one turns out to be better at resolving the particular
Fig. 2. Structure of eigenfunctions of growing and decaying modes. Real and imaginary parts of eigenfunctions are shown. (± ±, ± ±) Re 1000; (±±, - - -) Re 7500.
538
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
structure of this eigenmode. Finally, it is worth emphasizing that although a large number of points appear to be necessary to obtain the correct decay rates, there is no inconsistency in the present schemes as the solutions improve monotonically with increasing resolution. In all our simulations the grid spacing is uniform as that is the more appropriate spacing for turbulence simulations: small scale eddies occur everywhere and there are high frequency motions that mitigate against clustering. A turbulent fully developed channel ¯ow was also simulated and the mean velocity pro®le was found to follow the law of the wall and log-law approximately (see Ref. [1] for details). 3. Spot simulations In this section, we present some results of the spot simulations. Turbulent spots have been studied most often in laboratory ¯ows in boundary layers [9]. Spots have been also simulated in a boundary layer and a plane Poiseuille ¯ow [3]. Simulation [4] has preceded experiment [10] for zero pressure-gradient Couette ¯ow. But, analyses of these simulation data have not yet provided a clear understanding of the mechanisms at work in the evolution of spots. A way to search for the mechanisms would be to examine the behaviours in background ¯ows with dierent vorticity distributions. We have conducted spot simulations for a variety of Couette ¯ows with pressure gradients which gives dierent background vorticity distributions. To the best of our knowledge, these are the only simulations of spots in Couette ¯ows with pressure gradients. In a Cartesian frame with the x-axis aligned with the stream and y-axis normal to the parallel plates, a Couette ¯ow is speci®ed by the streamwise velocity component u y 12U 1 y 12C 1
y 2 ;
1
where U 1 is the speed of the upper plate at y 1, and C Re dp=dx. Dierent pressure gradients dp=dx are prescribed by assigning dierent values to C. Plane Poiseuille ¯ow is obtained by setting U 0 and C 2. The initial, localized perturbation to the velocity ®eld is given by a cross¯ow streamfunction w x; y; z 3 1
y 2 2 xz exp x2
z2 ;
from which cross¯ow velocity components v wz and w wy are obtained. A slightly dierent streamfunction was used to initialize the spot in the plane Poiseuille ¯ow. These dierent choices were necessary for comparing with previous work. This smooth disturbance evolves into disordered ¯ow which ®lls the region between the walls, but is otherwise localized. Simulations were conducted on grids of 128 33 128 points over a region of size 24p 2 16p in the streamwise, transverse and spanwise directions, and the time step was 0.04. The Reynolds number was 750. As in other DNS of turbulent ¯ows, the adequacy of the grid is not checked by grid re®nement in the usual way. Spot shape and propagation speed were found to be nearly the same on a smaller grid (128 33 64). Fourier expansion coecients of the velocity distribution along the streamwise coordinate begin decaying at intermediate wave numbers, indicating adequate resolution. Expansion coecients with respect to the spanwise coordinate suggest that slightly ®ner grids may still be necessary. As the calculations proceed without diculty, the spectrum would always show a decay towards high wave numbers, because the algo-
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
539
Fig. 3. Contours of the vorticity magnitude on midplane y 0 at time t 80. (a) Poiseuille ¯ow (b) Couette ¯ow at zero pressure gradient.
rithm provides for dissipation at high wave numbers. But when the decay begins at intermediate wave numbers, the eect of numerical dissipation is small. Grid resolution is deemed adequate in this sense. Due to practical limitations of computer resources, further re®nement must wait; re®nement will be done to test ideas formed from the data obtained from the present grids. On an R10000 SGI serial computer, each calculation requires approximately 40 h of CPU time. Fig. 3(a) and (b) shows spots at time t 80 in Poiseuille ¯ow and in a Couette ¯ow at zero pressure gradient. These are contours of the magnitude of the vorticity on midplanes y 0. Flow is from left to right. In each case, the shape of the spot is dierent. The characteristic arrowhead shape of the spot in the direction of ¯ow that has been observed in channel ¯ows [3] is evident in Fig. 3(a). At zero pressure gradient the spot must have fore-and-aft symmetry (Fig. 3(b)) and is roughly elliptical in planform. At very long times, the shape is expected to become circular [4]. Fig. 4 shows stages in the evolution at dierent times t 20, 40, 60 and 80. Here, the pressure gradient parameter in formula (1), C 1. The adverse gradient is then large enough to create a back¯ow. Due to fore±aft and top±bottom symmetry, there is no substantive dierence between a ¯ow with C 1 and one with C 1. But, in the present reference frame, C 1 implies a back¯ow. We know of no other simulation or experiment under these conditions with pressure gradients. Even these pictures suggest that a self-preserving state is likely as the same set of vortical structures seem to be present at the dierent times. We have not yet identi®ed substructures. The approach to self-preservation was con®rmed by noting that the area of the spot in a plane divided by the square of the length tends to a constant with time. Details are available [1]. In all three cases, there appears to be some degree of organization. It is strongly symmetric, as the Navier±Stokes equations preserves the symmetry of the initial data; special care has not been taken to obtain this
540
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
Fig. 4. Vorticity magnitude on midplane y 0 in adverse pressure gradient ¯ow C respectively.
1. (a)±(d) t 20, 40, 60 and 80,
Fig. 5. Spot propagation speeds of leading and trailing edges along centerline y 0, z 0 at dierent pressure gradients. ( ) leading edge; () trailing edge.
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
541
symmetry. Random, computing errors have not grown to destroy the symmetry. There appears to be a simple linear relationship between spot propagation rate and the applied pressure gradient. Fig. 5 shows propagation speeds of the front and rear edges of spots in the dierent ¯ows. 4. Conclusions We have examined the behaviour of a ®nite-dierence scheme for DNS of turbulent spots. The high-order upwind-biased scheme has the additional desirable property of stabilizing preliminary turbulence calculations that are usually on coarse grids. The solutions then improve as the grid is re®ned. This is an economical approach to DNS. The linear stability test demonstrated the accuracy of the method over long times when the grid is suciently ®ne to describe the dominant eigenfunction. However, this proves to be an uneconomical test: the growing mode has very large gradients near the walls which requires very many points for proper resolution. As the stability calculation is only a test, it would be preferable to validate the scheme with less computing eort. Spots in Couette ¯ows over a range pressure gradients were simulated with this method. The data show an organization of ¯ow which may not be very complicated and a tendency to reach a selfpreserving state. These results suggest that the spot is a rich turbulence phenomenon awaiting DNS studies. Acknowledgements We thank AR&DB (Aerodynamics Panel) for ®nancial support for a part of the work reported here. We acknowledge the use of computational facilities of the SERC and the AR&DB CFD Centre, IISc. for this work. References [1] Das A. A direct simulation study of turbulent spots in Couette ¯ows with pressure gradients. MSc (Engng) Thesis, Department of Aerospace Engineering, Indian Institute of Science, Bangalore, 1998. [2] Deshpande MD, Balaji S, Milton SG. Laminar ¯ow®elds in a three-dimensional cavity. PD CF-94-12, National Aerospace Laboratory, Bangalore, 1994. [3] Henningson D, Spalart P, Kim J. Numerical simulations of turbulent spots in plane Poiseuille and boundary layer ¯ows. Phys Fluids 1987;30(10):2914±7. [4] Lundbladh A, Johansson AV. Direct simulation of turbulent spots in plane Couette ¯ow. J Fluid Mech 1991;229:499±516. [5] Malik MR, Zang TA, Hussaini MY. A spectral collocation method for the Navier±Stokes equation. J Comp Phys 1985;61:64±88. [6] Moin P, Mahesh K. Direct numerical simulation: a tool in turbulence research. Ann Rev Fluid Mech 1998;30:539± 78. [7] Orszag SA. Accurate solution of the Orr±Sommerfeld stability equation. J Fluid Mech 1971;50:689±703. [8] Rai MM, Moin P. Direct simulations of turbulent ¯ow using ®nite-dierence schemes. J Comp Phys 1991;96:15±53. [9] Riley JJ, Gad-el-Hak M. The dynamics of turbulent spots. In: Davis SH, Lumley JL, editors. Frontiers in ¯uid mechanics, Berlin: Springer; 1984. p. 123±55. [10] Tillmark N, Alfredsson PH. Experiments on transition in plane Couette ¯ow. J Fluid Mech 1992;235:89±102.
www.elsevier.com/locate/comp¯uid
Direct numerical simulation of turbulent spots Arup Das 1, Joseph Mathew * Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560 012, India Accepted 10 November 2000
Abstract A group of high-order ®nite-dierence schemes for incompressible ¯ow was implemented to simulate the evolution of turbulent spots in channel ¯ows. The long-time accuracy of these schemes was tested by comparing the evolution of small disturbances to a plane channel ¯ow against the growth rate predicted by linear theory. When the perturbation is the unstable eigenfunction at a Reynolds number of 7500, the solution grows only if there are a comparatively large number of (equispaced) grid points across the channel. Fifth-order upwind biasing of convection terms is found to be worse than second-order central dierencing. But, for a decaying mode at a Reynolds number of 1000, about a fourth of the points suce to obtain the correct decay rate. We show that this is due to the comparatively high gradients in the unstable eigenfunction near the walls. So, high-wave-number dissipation of the high-order upwind biasing degrades the solution especially. But for a well-resolved calculation, the weak dissipation does not degrade solutions even over the very long times O 100 computed in these tests. Some new solutions of spot evolution in Couette ¯ows with pressure gradients are presented. The approach to self-similarity at long times can be seen readily in contour plots. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Turbulent spots; Direct numerical simulations; Stability
1. Introduction An isolated region of turbulence that is convected in an otherwise laminar ¯ow has been termed a spot. Spots have been observed in wall-bounded shear ¯ows such as boundary layers, (free surface) ®lm ¯ows and channel ¯ows. (see Ref. [9], for a review). In these ¯ows, transition begins in the region where spots appear. Spots grow as they travel downstream, and coalesce to give rise *
Corresponding author. Tel.: +91-80-309-2480; fax: +91-80-360-0683/+91-80-360-0134. E-mail address: [email protected] (J. Mathew). 1 Present address: Department of Theoretical and Applied Mechanics, University of Illinois at Urbana-Champaign, Urbana, IL 61820 USA. 0045-7930/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 3 0 ( 0 1 ) 0 0 0 0 4 - 4
534
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
to the fully turbulent region. Although several properties of spots have been documented, it is still not known whether sucient organization exists within turbulent spots to construct a simple model of their dynamics. It would be useful to explain the dierences in shapes and spreading of spots, observed in dierent background ¯ows, in a uni®ed way. The search for simple organization within spots can be facilitated by conducting direct numerical simulations (DNS), which can provide complete ®eld data ± especially vorticity ®elds that are not easily obtained in experiments. Often, resolution and geometry requirements have limited DNS to studies in homogeneous turbulence or other simple ¯ows. But, as a naturally isolated region of turbulence, the spot presents a particularly suitable ¯ow for DNS studies. We have applied variants of a method presented by Rai and Moin [8] to simulate spot evolution in channel ¯ows driven by shear and by pressure gradients. The method involves high-order upwind-biased ®nite dierences which gives high-order dissipation. Good results were demonstrated in the original presentation [8], but the method has since then fallen out of favour because its high-wave-number dissipation is unsuitable for Large Eddy Simulations (LES) [6]. Rai and Moin [8] have already noted the bene®cial characteristics. We expect some speci®c bene®ts. As a ®nite-dierence method, we expect it to be extended easily to calculations in practical geometries of small curvature such as wing±body corners. As a practical matter, high-order dissipation which is eective only at the smallest scales can be helpful in stabilizing DNS. For economy of overall eort, it is common to begin DNS studies with grids on which the solutions are slightly underresolved, examine the solutions, and then employ ®ner grids. The merit of high-order dissipation is that it stabilizes the calculations to provide some preliminary results on small grids and then provides better solutions on ®ner grids. In this paper we focus attention on a development experience. To understand the long-time ®delity of the method, knowing that it has weak dissipation, we applied the method to compute the evolution of disturbances to an unstable Poiseuille ¯ow. It turns out that the unstable eigenfunction has a pathological, rich structure close to the channel wall which requires very ®ne grids. The evolution on coarse grids showed decay of perturbation energy which can be mistakenly attributed to numerical dissipation. However, the numerical method is able to compute the evolution of decaying modes for very long times accurately, even on such grids, because the eigenfunction is smooth. We conclude with some aspects of the solutions to spot simulations. 2. Numerical method DNS have often employed spectral methods, but increasingly, it is becoming clear that ®nitedierence and ®nite volume methods can be used successfully when the algorithms are constructed with sucient care. Thus Rai and Moin [8] presented a general, ®nite-dierence scheme using high-order upwind-biased dierence formulae. High-order formulae are used to compete with spectral methods and keep phase errors small over a larger wave number range. Upwind biasing, rather than upwinding, permits compact stencils while retaining high-order dissipation. DNS implies an attempt to resolve all signi®cant scales of motion, down to the scales where dissipation occurs. So the weak, high-order dissipation should be acceptable as demonstrated in Ref. [8]. But, for LES, as very little energy dissipation should occur at the smallest resolved scales, this method is unsuitable [6].
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
535
We integrated the Navier±Stokes equations by taking the ¯ow to be at constant density and viscosity. The numerical method used here involves only minor variations to that presented in Ref. [8]. As a splitting method, the ®rst step is an explicit integration of the momentum equations to ®nd an intermediate velocity ®eld, with a given pressure ®eld. The second, corrector step ®nds the pressure distribution and the correction to the intermediate velocity ®eld that makes the resultant ®eld divergence free. The corrector step requires the solution of the coupled system of equations for pressure. In the present problems involving channel ¯ows, the solutions are periodic in the streamwise and spanwise directions. So we use Fourier expansions to reduce the problem and then use a direct, banded-matrix solver. The discretization equations are not repeated here as Ref. [2] is widely accessible. While developing the codes, the ®rst-order Euler and third-order Runge±Kutta methods were used for time stepping. Convection terms were represented with second-order central, third-order upwind-biased and ®fth-order upwind-biased dierence formulae while second-order and sixthorder central dierences were used for viscous terms. The results of spot evolution presented here were obtained using Euler time stepping, ®fth-order upwind biasing of convection terms and second-order central dierencing of viscous terms. The numerical divergence and gradient operations on the velocity and pressure ®elds that are done in the corrector step were done with fourthorder central dierences. The main dierences between the present implementations and that in Ref. [8] are: (1) we use uniform grids, and (2) we use linear interpolation to ®nd velocity crossproduct terms rather than the fourth-order cubic interpolation. However, for the subject of this paper, the use of uniform grids has had the more signi®cant eect than the order of the scheme. All quantities are non-dimensionalized using the channel halfwidth as the length scale. The velocity scale is the velocity of the moving plate in all Couette ¯ows, and the centerline velocity for Poiseuille ¯ow. The Reynolds number Re is always constructed from these two scales and the ¯uid viscosity. 2.1. Linear stability test for long-time accuracy The initial tests of the codes for steady state solutions to the 2-D and 3-D lid-driven cavity at low Reynolds numbers were satisfactory and presented no surprises. Solutions to 3-D cavity ¯ow are few, but ours compared closely with those of Deshpande et al. [2]. The comparisons are in Ref. [1]. A more appropriate test of a numerical method which has been designed to compute turbulent ¯ows accurately must evaluate the ®delity of the computations against a known unsteady ¯ow. The growth of prescribed small amplitude disturbances to a Poiseuille ¯ow is a natural choice as solutions can be found accurately by solving Orr±Sommerfeld eigenvalue problems. Moreover, though obtained numerically, these solutions are well established. As Rai and Moin [8] had reported good comparisons with this test problem as well, this test was also expected to be routine. The temporal stability of the plane Poiseuille ¯ow has been studied by Orszag [7]. All modes are damped when the Reynolds number is below 5772. We considered the mode of wave number unity and Re 7500 for which the complex frequency x 0:24989154 i0:00223498. The eigenfunction was obtained using a spectral code and superposed on the parabolic Poiseuille ¯ow. The perturbation amplitude is 0.0001. Fig. 1(a) shows the growth of the energy of the perturbation velocity ®eld on three grids: 32 65, 32 129 and 32 257. In each case, there are 32 equispaced grid points in the
536
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
streamwise direction and an increasingly ®ner equispaced grid normal to the walls. This was computed with ®fth-order upwind-biased convection terms and second-order central dierenced viscous terms. The result of linear theory, which is the exponential growth of the perturbation is
Fig. 1. Evolution of perturbation energy with time. Graphs (a) and (b) are of growing mode at Re 7500; graph (c) is of the decaying mode at Re 1000. The scheme used is ®fth-order upwind-biased convection and second-order central dierenced viscous terms for results in graphs (a) and (c); all terms are with second-order central dierences in (b). The solid line in each graph is the reference (Orr±Sommerfeld) solution. Curves in graphs (a) and (b) are from the dierent grids: (- - -) 32 65; (± ± ±) 32 129; ( ) 32 257. The dashed curve in (c) is from a 32 65 grid.
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
537
the straight line. It is clear that even 257 points are not sucient to simulate time-accurate evolution as even the slopes at long times are noticeably dierent. At lower resolutions, the solutions are decaying. Also, the initial evolution is quite dierent from that expected from linear theory. The corresponding curves in Fig. 1(b) are from a calculation that used second-order central dierencing of all terms and the long-time results are better. The dierences at short times remain. Fig. 1(c) is from a similar calculation, but the initial perturbation is the damped eigenfunction of wave number unity at Re 1000. Now, even the initial dierences are minor, and a 32 65 grid and the ®fth-order upwind-biased scheme is able to capture the correct decay rate. It is tempting to dismiss the upwind-biased scheme as being too dissipative by comparing Fig. 1(a) and (b), whereas Fig. 1(c) shows that the method can remain accurate over long times. Simulating the growing mode turns out to be a specially dicult problem. There are two aspects to the diculties, and both are related to the solution, i.e., the structure of the eigenfunctions. First, the growing mode is not adequately resolved until a comparatively large number of grid points are used. Secondly, the ®ne structure exaggerates the high-wave-number dissipation of the upwind-biased method compared to the dissipation-free central-dierence method. The eigenfunction has a very rich structure near the walls and unless those variations are resolved, it is not possible to simulate its behaviour accurately. Fig. 2 shows the structure of the eigenfunctions W y of the growing (Re 7500) and decaying (Re 1000) modes close to a wall. The third derivative of W y, which appears as a viscous term, has been plotted to see the sharper gradients of the growing mode. Signi®cantly more points are required in a uniform-grid computation to resolve these gradients. Rai and Moin [8] were able to get good solutions on a 32 65 grid by clustering wall-normal grid points. The spacing increased in geometric progression from the wall. It turns out that the choice of clustering is also important: even when Chebyshev points were used, Malik et al. [5] showed that about 256 points across the channel are needed to reproduce an acceptable growth rate. They used a second-order ®nite-dierence scheme. Both schemes provide clustering near the wall, but one turns out to be better at resolving the particular
Fig. 2. Structure of eigenfunctions of growing and decaying modes. Real and imaginary parts of eigenfunctions are shown. (± ±, ± ±) Re 1000; (±±, - - -) Re 7500.
538
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
structure of this eigenmode. Finally, it is worth emphasizing that although a large number of points appear to be necessary to obtain the correct decay rates, there is no inconsistency in the present schemes as the solutions improve monotonically with increasing resolution. In all our simulations the grid spacing is uniform as that is the more appropriate spacing for turbulence simulations: small scale eddies occur everywhere and there are high frequency motions that mitigate against clustering. A turbulent fully developed channel ¯ow was also simulated and the mean velocity pro®le was found to follow the law of the wall and log-law approximately (see Ref. [1] for details). 3. Spot simulations In this section, we present some results of the spot simulations. Turbulent spots have been studied most often in laboratory ¯ows in boundary layers [9]. Spots have been also simulated in a boundary layer and a plane Poiseuille ¯ow [3]. Simulation [4] has preceded experiment [10] for zero pressure-gradient Couette ¯ow. But, analyses of these simulation data have not yet provided a clear understanding of the mechanisms at work in the evolution of spots. A way to search for the mechanisms would be to examine the behaviours in background ¯ows with dierent vorticity distributions. We have conducted spot simulations for a variety of Couette ¯ows with pressure gradients which gives dierent background vorticity distributions. To the best of our knowledge, these are the only simulations of spots in Couette ¯ows with pressure gradients. In a Cartesian frame with the x-axis aligned with the stream and y-axis normal to the parallel plates, a Couette ¯ow is speci®ed by the streamwise velocity component u y 12U 1 y 12C 1
y 2 ;
1
where U 1 is the speed of the upper plate at y 1, and C Re dp=dx. Dierent pressure gradients dp=dx are prescribed by assigning dierent values to C. Plane Poiseuille ¯ow is obtained by setting U 0 and C 2. The initial, localized perturbation to the velocity ®eld is given by a cross¯ow streamfunction w x; y; z 3 1
y 2 2 xz exp x2
z2 ;
from which cross¯ow velocity components v wz and w wy are obtained. A slightly dierent streamfunction was used to initialize the spot in the plane Poiseuille ¯ow. These dierent choices were necessary for comparing with previous work. This smooth disturbance evolves into disordered ¯ow which ®lls the region between the walls, but is otherwise localized. Simulations were conducted on grids of 128 33 128 points over a region of size 24p 2 16p in the streamwise, transverse and spanwise directions, and the time step was 0.04. The Reynolds number was 750. As in other DNS of turbulent ¯ows, the adequacy of the grid is not checked by grid re®nement in the usual way. Spot shape and propagation speed were found to be nearly the same on a smaller grid (128 33 64). Fourier expansion coecients of the velocity distribution along the streamwise coordinate begin decaying at intermediate wave numbers, indicating adequate resolution. Expansion coecients with respect to the spanwise coordinate suggest that slightly ®ner grids may still be necessary. As the calculations proceed without diculty, the spectrum would always show a decay towards high wave numbers, because the algo-
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
539
Fig. 3. Contours of the vorticity magnitude on midplane y 0 at time t 80. (a) Poiseuille ¯ow (b) Couette ¯ow at zero pressure gradient.
rithm provides for dissipation at high wave numbers. But when the decay begins at intermediate wave numbers, the eect of numerical dissipation is small. Grid resolution is deemed adequate in this sense. Due to practical limitations of computer resources, further re®nement must wait; re®nement will be done to test ideas formed from the data obtained from the present grids. On an R10000 SGI serial computer, each calculation requires approximately 40 h of CPU time. Fig. 3(a) and (b) shows spots at time t 80 in Poiseuille ¯ow and in a Couette ¯ow at zero pressure gradient. These are contours of the magnitude of the vorticity on midplanes y 0. Flow is from left to right. In each case, the shape of the spot is dierent. The characteristic arrowhead shape of the spot in the direction of ¯ow that has been observed in channel ¯ows [3] is evident in Fig. 3(a). At zero pressure gradient the spot must have fore-and-aft symmetry (Fig. 3(b)) and is roughly elliptical in planform. At very long times, the shape is expected to become circular [4]. Fig. 4 shows stages in the evolution at dierent times t 20, 40, 60 and 80. Here, the pressure gradient parameter in formula (1), C 1. The adverse gradient is then large enough to create a back¯ow. Due to fore±aft and top±bottom symmetry, there is no substantive dierence between a ¯ow with C 1 and one with C 1. But, in the present reference frame, C 1 implies a back¯ow. We know of no other simulation or experiment under these conditions with pressure gradients. Even these pictures suggest that a self-preserving state is likely as the same set of vortical structures seem to be present at the dierent times. We have not yet identi®ed substructures. The approach to self-preservation was con®rmed by noting that the area of the spot in a plane divided by the square of the length tends to a constant with time. Details are available [1]. In all three cases, there appears to be some degree of organization. It is strongly symmetric, as the Navier±Stokes equations preserves the symmetry of the initial data; special care has not been taken to obtain this
540
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
Fig. 4. Vorticity magnitude on midplane y 0 in adverse pressure gradient ¯ow C respectively.
1. (a)±(d) t 20, 40, 60 and 80,
Fig. 5. Spot propagation speeds of leading and trailing edges along centerline y 0, z 0 at dierent pressure gradients. ( ) leading edge; () trailing edge.
A. Das, J. Mathew / Computers & Fluids 30 (2001) 533±541
541
symmetry. Random, computing errors have not grown to destroy the symmetry. There appears to be a simple linear relationship between spot propagation rate and the applied pressure gradient. Fig. 5 shows propagation speeds of the front and rear edges of spots in the dierent ¯ows. 4. Conclusions We have examined the behaviour of a ®nite-dierence scheme for DNS of turbulent spots. The high-order upwind-biased scheme has the additional desirable property of stabilizing preliminary turbulence calculations that are usually on coarse grids. The solutions then improve as the grid is re®ned. This is an economical approach to DNS. The linear stability test demonstrated the accuracy of the method over long times when the grid is suciently ®ne to describe the dominant eigenfunction. However, this proves to be an uneconomical test: the growing mode has very large gradients near the walls which requires very many points for proper resolution. As the stability calculation is only a test, it would be preferable to validate the scheme with less computing eort. Spots in Couette ¯ows over a range pressure gradients were simulated with this method. The data show an organization of ¯ow which may not be very complicated and a tendency to reach a selfpreserving state. These results suggest that the spot is a rich turbulence phenomenon awaiting DNS studies. Acknowledgements We thank AR&DB (Aerodynamics Panel) for ®nancial support for a part of the work reported here. We acknowledge the use of computational facilities of the SERC and the AR&DB CFD Centre, IISc. for this work. References [1] Das A. A direct simulation study of turbulent spots in Couette ¯ows with pressure gradients. MSc (Engng) Thesis, Department of Aerospace Engineering, Indian Institute of Science, Bangalore, 1998. [2] Deshpande MD, Balaji S, Milton SG. Laminar ¯ow®elds in a three-dimensional cavity. PD CF-94-12, National Aerospace Laboratory, Bangalore, 1994. [3] Henningson D, Spalart P, Kim J. Numerical simulations of turbulent spots in plane Poiseuille and boundary layer ¯ows. Phys Fluids 1987;30(10):2914±7. [4] Lundbladh A, Johansson AV. Direct simulation of turbulent spots in plane Couette ¯ow. J Fluid Mech 1991;229:499±516. [5] Malik MR, Zang TA, Hussaini MY. A spectral collocation method for the Navier±Stokes equation. J Comp Phys 1985;61:64±88. [6] Moin P, Mahesh K. Direct numerical simulation: a tool in turbulence research. Ann Rev Fluid Mech 1998;30:539± 78. [7] Orszag SA. Accurate solution of the Orr±Sommerfeld stability equation. J Fluid Mech 1971;50:689±703. [8] Rai MM, Moin P. Direct simulations of turbulent ¯ow using ®nite-dierence schemes. J Comp Phys 1991;96:15±53. [9] Riley JJ, Gad-el-Hak M. The dynamics of turbulent spots. In: Davis SH, Lumley JL, editors. Frontiers in ¯uid mechanics, Berlin: Springer; 1984. p. 123±55. [10] Tillmark N, Alfredsson PH. Experiments on transition in plane Couette ¯ow. J Fluid Mech 1992;235:89±102.