Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 41, 1, pp. 19-32, Warsaw 2003 INFLUENCE OF THE COMPACT EXPLICIT FILTERING METHOD ON THE PERTURBATIONS GROWTH IN TEMPORAL SHEAR-LAYER FLOW Artur Tyliszczak1 Institute of Thermal Machinery, Technical University of Częstochowa e-mail: atyl@imc.pcz.czest.pl In this paper the effect of the explicit filtering of high frequency pheno- mena performed bymeans of the high order compact filters is presented. A 2D unsteady incompressible flow in which the perturbations super- imposed on the mean velocity profile evolve in time is considered. The dependence of the amplitudes of disturbances is analysed for various fil- tering procedures and comparedwith othermethods intended to remove high frequency waves. Key words: compact scheme, filtering, aliasing removal 1. Introduction The widespread use and continuous development of high-speed compu- ters turned numerical simulation a commonpractice for solving turbulent flow problems. This holds both for the practical industrial flows as well as for the fundamental studies where the physics of turbulence may be analysed either by direct numerical simulation (DNS) or by large eddy numerical simulation (LES). Both these methods require solution of time and length scales that differ bymany orders of magnitude, starting from the large vortex structures of size of the computational domain up to the Kolmogorow scale in the case of DNS and the so-called cut-off scale for the LES method. From this point of view, the classical approaches based on the finite differences, finite volume or finite element methods (Hirsh, 1990) are not suitable due to their limited 1The author won the second prize awarded at the biennial young researchers’ contest for the bestworkpresented at the 15thPolishConference onFluidMechanics, Augustów, September 2002 20 A.Tyliszczak possibilities of solving high frequency phenomena, or eventually, their appli- cation would require a huge number of computational nodes increasing the computational time and would also require more powerful computers. There- fore, for theapplication ofDNSandLES, themost appropriate seems tobe the spectral and pseudo-spectralmethods (Canuto et al., 1988) for which the high resolution is obtained at a relatively small number of the computational nodes. One drawback of thesemethods is their limitation to the simple geometry and dependence on the type of boundary conditions.An alternate approach,which ensures spectral-like resolution is the compact scheme introduced in the 1930s. The fundamental ideas as well as their derivation techniques are, among the others, given in Vichnevetsky and Bowles (1982) and in the seminar paper of Lele (1992). The implementation of the compact scheme is relatively simple and is not restricted to the uniform grid and to some type of the bounda- ry conditions. Moreover, it is generally more robust and less computationally costly than the spectral or pseudo-spectral methods. The succesfull application of the compact scheme itself or its combination with the spectral and classical discretization methods has already been pre- sented bymany authors (Lele, 1992; Boersma andLele, 1999; Cook andRiley, 1996; Metais et al., 1999; Tyliszczak, 2002) for various fluid flow problems. Theweak point of the compact schemes is that they do not share the so-called telescoping property, and therefore they are not conservative. It means that the systemrepresenting conservation laws at thepartial differential level, when discretized by means of the compact schemes, will not reflect its conservative form at the discretization level. Moreover, except for the boundary and near boundarynodes the compact schemes have purely central character whichme- ans that they do not have dissipative properties, as for example the upwind biased scheme does. In effect, high frequency waves, arising either due to the round-of errors or the aliasing errors coming from nonlinear interaction, will grow and affect a low frequency wave what usually leads to instability, and the solution diverge. Therefore, for the computational stability it is necessary to remove or damp such errors what in practice may be done by high-pass explicit filtering which eliminate high frequency waves regardless of their so- urce or in the case of the aliasing errors, theymay beminimized by using the Arakawa (Arakawa, 1996; Peyret, 2002) skew-symmetric formof the nonlinear terms. Within this paper, based on the two-dimensional incompressible tem- poral shear-layer test case we compare the influence of both these methods on the growth of the disturbances in view of their amplitudes and the time of occurrence of their maximum. The flow problem chosen for the analysis assumes periodicity in one space direction allowing one in this way to use also Influence of the compact explicit filtering method... 21 the pseudo-spectral approximation, which will be applied for comparison. In this case, we additionally consider themethod for removing the aliasing errors by applying the 2/3 law (Canuto et al., 1988). 2. Governing equations and flow configuration TheNavier-Stokes equations for the incompressible flow together with the continuity equation are given in a nondimensional form as ∂ui ∂t + ∂uiuj ∂xj =− 1 ρ ∂p ∂xi + 1 Re ∂τij ∂xj (2.1) ∂ui ∂xi =0 where τij = ∂ui ∂xj + ∂uj ∂xi − 2 3 δij ∂uk ∂xk (2.2) where Re is the Reynolds number, ui stands for the velocity component whi- le ρ and p denote the density and pressure, respectively. The Navier-Stokes equations (2.1)1 are discretized on the uniform collocated grid with III-rd order Adams-Bashforthmethod for time integration. The continuity equation is enforced according to the projection method in connection with the algori- thm of Ku et al. (1987) allowing one to fullfil the continuity equation in every computational node including the boundary nodes as well. The flow problem considered in this work together with the prescribed boundary conditions and symbolically sketchedmean velocity profile is shown inFig.1. Thewalls at the upper and lower boundarymovewith the speed ±1, respectively. The mean velocity profile was taken the same as in the work of Moser andRogers (1993) which is used for comparison and verification of the obtained results. It is given as u(y)=Uerf(y √ π) (2.3) where U is the velocity outside the mixing layer. The initial disturbance su- perimposed on themean velocity profilewas introduced in two different ways. First, we wanted to analyse the evolution of the streamwise velocity ener- gy spectrum in dependence of the method used to remove high frequency 22 A.Tyliszczak Fig. 1. Temporal shear-layer flow phenomena. In this case, the initial velocity field was superimposed with the disturbances defined as u′(y)= N/4∑ k=1 exp(−πy2)[A1 cos(kx+φ)+A2 sin(kx+φ)] (2.4) where φ is the random phase shift and the amplitudes A1,A2 are equal to A1 =A2 = 1 20 + 1 100 ( ran− 1 2 ) (2.5) where ran denotes random number in range [0,1]. The one-dimensional energy spectrum of the streamwise velocity compo- nent is defined as E1(k)= 1 L L/2∫ −L/2 |û(k,y)|2 dy (2.6) where û(k,y) is the k-th Fourier mode. In thenextpart of the computationsweanalysed the evolution of theparti- cularmodes explicitly defined at the initial state. In this casewe disturbed the initial vorticity field ω(x,y) = ω(x,y)+ω′(x,y), where ω(x,y) = ∂u(y)/∂y. Influence of the compact explicit filtering method... 23 The vorticity disturbance is defined as ω′(x,y)= ∑ α=1,1 2 A0αReal [ φα(y)e i(kαx) ] (2.7) where the symbols α=1, 1 2 are given in the form α= kαλ 2π (2.8) where kα and λdenote thewave number andwavelength of themost unstable mode and its subharmonics. The length of the most unstable disturbance is equal to λ = 2.32π (Moser and Rogers, 1993). The symbol A0α appearing in Eq. (2.7) denotes the amplitude of a given mode, while φα(y) represents the vorticity eigenfunction obtained from the solution of the Orr-Sommerfeld equation for the wave numbers defined according to the parameter α. The length of the computational domain Lwas equal four times the length of the themost unstablemode, allowing it to form four vortices corresponding to the most unstablemode, whichwas then subject to pairing due to the presence of the subharmonic modes (Moser and Rogers, 1993; Lesieur et al., 1988). The amplitudes of a given mode of the initial perturbation or in the evolved field are measured by the integrated RMS of the velocity modes. Thus they are defined as Aα = [+L/2∫ −L/2 2ûi(α)û ∗ i(α) dy ]1 2 (2.9) where ûi(α) is the Fourier coefficient of the velocity component ui and û ∗ i(α) is its conjugate value. 3. Space discretization method In the computation performed, in the periodic directionwe have used both the compact scheme and, also for comparison, the pseudo-spectral method (Canuto et al., 1988) based on the Fourier series expansion. The approxima- tion of the first derivative bymeans of the compact scheme, with the assumed uniformly spaced mesh and the nodes indexed by i, are given below. The se- condderivativeswere computedby successive application of thefirstderivative approximation. 24 A.Tyliszczak Fornode i the relationbetween thevalues of a function f and its derivative f ′ is defined [?] by a linear combination of both the function f and also f ′ as βf ′i−2+αf ′ i−1+f ′ i +αf ′ i+1+βf ′ i+2 = (3.1) = c fi+3−fi−3 6h + b fi+2−fi−2 4h +a fi+1−fi−1 2h where h is themesh size. The values of the coefficients α and β together with the coefficients a, b, and c are obtained by expanding the function f and its derivative f ′ into Taylor series around node i, and next bymatching various orders of these expansions.Note thatby setting α=β=0the explicit formula for the first derivative is obtained which formally can provide a sixth-order scheme. As shown by Lele (1992), the highest order compact approximation which can be obtained from Eq. (3.1) is of tenth order. In our computations we used the sixth-order scheme with the following values of coefficients α= 1 3 β=0 a= 14 9 b= 1 9 c=0 (3.2) Thevalues of derivatives are obtainedby solving the linear systemof equations given as Af ′ =Bf (3.3) where thematrices A and B consist of the coefficients defined in Eq. (3.2). In the case of the periodic direction, thematrices A and B are of a periodic type. For the non-periodic boundaries the crucial issue of the compact schemes is the correct approximation at the boundary and near boundary nodes. While formula given by Eq. (3.1) with the set of coefficients (3.2) can be applied for the inner points starting from the node i = 3 up to i = N − 2, the approximations on i = 1,2 and i = N,N − 1 need another treatment. We note that awrongchoice of theapproximationon thesenodes causes theoverall approximation to violate the condition of stability based on the eigenvalues of the spatial operator – there exist eigenvalues on the right side in the complex plane of the Fourier footprint. This aspect is discussed in detail in Carpenter et al. (1993). According to Gustafsson (1975), who showed that the formal overall accuracy of the scheme can be atmost one order higher than the order of the approximation at the boundary nodes, for our sixth order inner scheme it is advisable to use at the boundary nodes an approximation of at least fifth order. One of the possibilities to construct fifth order closure formulas at the Influence of the compact explicit filtering method... 25 boundary nodes is a onesided (Lele, 1992) compact approach which can be applied to the first and second grid nodes. However, in this case the scheme violates the stability condition and can not be used. The highest possible order of the approximation obtained with compact formulas which produce a correct eigenvalue spectrum is of the form (3-4-6-4-3), where the first and last numbers stand for the order of the approximation on the boundarynodes, second and before last numbers for near boundary nodes, and the central number represents the order of the inner approximation. Thus, the overall formal approximation is of fourth order. The equations for the derivatives at the boundary nodes are given as f ′1+2f ′ 2 = 1 h ( − 15 6 f1+2f2+ 1 2 f3 ) (3.4) f ′N +2f ′ N−1 = 1 h (15 6 fN −2fN−1− 1 2 fN−2 ) and for the second and before last node the equation is f ′i−1+4f ′ i +f ′ i+1 =3 fi+1−fi−1 h (3.5) 4. Filtering procedure For the filtering operation we have chosen a compact filter defined (Lele, 1992) as βf̂i−2+αf̂i−1+ f̂i+αf̂i+1+βf̂i+2 = (4.1) = afi+ d 2 (fi+3+fi−3)+ c 2 (fi+2+fi−2)+ b 2 (fi+1+fi−1) where f̂ stands for the filtered value at a given node.Writing the above equ- ation for each computational node forms a system of linear equations, and their solution gives the vector of the filtered variables. The transfer function G(ω) of this filter is given as G(ω)= a+ bcosω+ ccos2ω+dcos3ω 1+2αcosω+2β cos2ω (4.2) where ω=2πk/N is the so-called scaledwavenumber,which for 0¬ k¬N/2 belongs to the range [0,π]. The Nyquist wave number corresponds to ω=π. 26 A.Tyliszczak Thenecessary condition for the filter is to remove all the phenomenaoccurring at ω=π as they can not be solved on a givenmesh. Therefore the coefficients appearing inEq. (4.1) are determinedbymatching various orders of theTaylor series expansions with the additional requirement for the transfer function G(π) = 0. In the computations we have used the filters of the IV -th and VI-th order. The first one is defined by the following set of coefficients a= 1 8 (5+6α) b= 1 2 (1+2α) (4.3) c=−1 8 (1−2α) where α varies in range α = 0.4-0.45. The VI-th order filter is desi- gned with the additional constrains for the transfer function in the form (d2G/dω2)(π) = (d4G/dω4)(π) = 0. These conditions cause effective dam- ping of the wave numbers also close to ω = π. The coefficients for this filter are defined as α=0 β= 3 10 a= 1 2 b= 3 4 c= 3 10 d= 1 20 (4.4) Fig. 2. Transfer function G(ω) for the IV -th order filter with α=0.4 and for the VI-th order filter The transfer functions for the IV -th and VI-th order filter are shown in Fig.2 where, as one may see, both filters have the property G(π) = 0. More- Influence of the compact explicit filtering method... 27 over, the sixth-order filter will additionally damp the phenomena occurring at the high wave numbers. In the case of non-periodic boundaries the filter given byEq. (4.1)must be complemented by additional formulas for the boundary nodes. We have used an explicit filter of the form f̂1 = 15 16 f1+ 1 16 (4f2−6f3+4f4−f5) (4.5) f̂2 = 3 4 f2+ 1 16 (f1+6f3−4f4+f5) with similar expressions for the nodes N−1 and N. The filter defined in Eqs (4.5) is of the IV -th order and has the property G(π)= 0. 5. Results and discussion The computational mesh consists of 256× 256 uniformly spaced nodes. The time step used in the calculation was equal to 0.01, theReynolds number based on the vorticity thickness and the velocity outside themixing layer was equal to 250. The amplitudes of themost unstablemode and its subharmonic mode computed according to Eq. (2.9) were equal to 0.04 and 0.03, respec- tively. The computations were performed for two discretization methods, first using the compact scheme in both directions (method denoted wherein as CD-CD) and also using the compact scheme in the vertical direction together with the pseudo-spectral methods in the horizontal direction (method deno- tedwherein asCD-PS). In the results presented below, the filtering procedure was performed every five time steps. The influence of the filtering procedure interval is not analysed in this work, but it must be stressed that it has no influence on the general trends obtained in the results presented. However, we note thatwithout thefilteringorwithoutusinganothermethods (e.g.Arakawa skew-symmetric form of the non-linear terms, 2/3 law) for removing the high frequency waves, the computations are unstable. In Fig.3 we present sample results for the computations performedwith the CD-CDmethod showing the evolution of the vorticity. Depending on the discretization method in a given direction, the high frequency waves have been removed in the following manner. In case of the CD-PS discretization, in the horizontal periodic direction the computations were performed with the VI-th order periodic filter, the 2/3 law and also 28 A.Tyliszczak Fig. 3. The vorticity izolines for time instants t=0,12,22,30 Fig. 4. Evolution of the energy spectrum for computations performed with the VI-th order periodic filter with theArakawa form of the nonlinear terms. In the vertical directionwe ap- plied the VI-th order filter. The computations with the CD-PS discretization method were performed for the disturbances introduced randomly as defined in Eq. (2.4). The energy spectra obtained for these computations are shown in Fig.4 to Fig.6. As one may see, the explicit filtering, see Fig.4, has the best properties in the sense of removing the high frequency waves. The spectrum is smooth and there is no characteristic kink which appears at the high wave numbers in the case of the two remainingmethods. The reason for this is that neither applying the 2/3 law nor Arakawa form of the nonlinear terms takes into account the viscous terms of the Navier-Stokes equations, which in fact may also produce the unresolved scale. From this point of view the explicit filtering, which does not distinguish between the source of the high frequency Influence of the compact explicit filtering method... 29 phenomena, seems to be the most appropriate, however in comparison with the application of the 2/3 law, the filtering procedure is computationally a bit more costly. Fig. 5. Evolution of the energy spectrum for computations performed with the application of the 2/3 law Fig. 6. Evolution of the energy spectrum for computations performed with the Arakawa form of the nonlinear terms For the computations performed with disturbances introduced based on the eigenfunction of the Orr-Sommerfeld equation both CD-CD and CD-PS discretizationmethods give practically the same results, regardless of the pro- cedure used for removing the disturbances in the periodic direction. The re- sults presented in the following, concern theCD-CDdiscretizationmethod for which in the horizontal direction we used the VI-th order filter denoted as FVIx while in the vertical direction we applied the Arakawa form of the non- linear terms denoted as AR and the IV -th and VI-th order filters denoted as FIVy and FVIy . The obtained results presenting the time evolution of the 30 A.Tyliszczak amplitude of disturbances are shown in Fig.7. The detailed data concerning the occurrence times of the first local maximum of the most unstable mode (denoted as t1) and its subharmonic mode (denoted as t1/2) together with the values of the amplitudes measured at the corresponding times are given in Table 1. These results are in good agreement with the data given in Mo- ser and Rogers (1993) where the same flow problem has been solved with the pseudo-spectral approximation in both directions based on theFourier/Jacobi series expansions. Fig. 7. Growth of the amplitude disturbances for the discretizationmethod CD-CD for variousmethods used to eliminate high frequency waves Table 1. The amplitudes A1 and A1/2 and the occurrence time of their maxima for the CD-CD discretization method A1 t1 A1/2 t1/2 AR−FVIx 0.707 12.50 1.727 22.50 FIVy −FVIx 0.591 12.25 1.744 22.50 FVIy −FVIx 0.547 12.50 1.760 22.25 Influence of the compact explicit filtering method... 31 In this paper there was no data concerning the values of the amplitudes, while thepresented times t1 and t1/2were equal to 11.9and 21.5, respectively. The results shown in Fig.7 confirm good properties of the filtering procedure. In comparison to the results obtained with theArakawa form of the nonlinear terms, the amplitude of the high frequency disturbance A1, which in this case is the most unstable mode, is decreased. Furthermore, applying the VI-th order filter, which has the property to completely remove not only the waves corresponding to Nyquist frequency but also the waves close to it, see Fig.2, allows for the most efficient damping of the high frequency phenomena. The important thing is that none of the filters considerably influence neither the time of the local amplitudemaxima nor the low frequency waves. The ampli- tude of the subharmonic disturbance is almost the same for all the method applied. This work was supported by State Committee for Scientific Research, Poland, under grant No. 8T10B00919. References 1. ArakawaA., 1966,Computational design for long-termnumerical integration of the equations of fluid motion: Two dimensional incompressible flow. Part 1, Journal of Computational Physics, 1, 119-143 2. Boersma B.J., Lele S.K., 1999, Large eddy simulation of compressible jet, Annual Research Briefs, Stanford University Center for Turbulence Research 3. Canuto C., Hussaini M.Y., Quarteroni A., Zang T.A., 1988, Spectral Methods in Fluid Dynamics, Springer-Verlag 4. Carpenter M.H., Gottlieb D., Abarbanel S., 1993, The stability of nu- merical boundary treatements for compact high-order finite-difference schemes, Journal of Computational Physics, 108, 272-295 5. Cook W.C., Riley J.J., 1996, Direct numerical simulation of a turbulent reactive plume on a parallel computer, Journal of Computational Physics, 129, 263-283 6. Gustafsson B., 1975, The convergence rate for difference approximations to mixed initial boundary value problems,Math. Comput., 29, 396-406 7. HirshCh., 1990,Numerical Computation of Internal andExternal Flows, John Wiley & Sons, Chichester 32 A.Tyliszczak 8. KuH.C., HirshR.S., TaylorT.D., 1987, Pseudospectralmethods for solu- tion of the incompressible Navier-Stokes equations,Computers and Fluids, 15, 195-214 9. Lele S.K., 1992, Compact finite difference with spectral-like resolution, Jour- nal of Computational Physics, 103, 16-42 10. Lesieur M., Staquet Ch., Le Roy P., Comte P., 1988, Themixing layer and its coherence examined from the point of view of two-dimensional turbu- lence, Journal of Fluid Mechanics, 192, 511-534 11. Metais O., Lesieur M., Comte P., 1999,Transition, Turbulence and Com- bastion Modelling, Chapter – Large-Eddy Simulations of Incompressible and Compressible Turbulence, KLUWERAcademic Publisher 12. MoserR.D.,RogersM.M., 1993,The three-dimensional evolutionof aplane mixing layer: pairing and transition to turbulence, Journal of FluidMechanics, 247, 275-320 13. PeyretR., 2002,SpectralMethods for Incompressible Viscous Flow, Springer- Verlag NewYork 14. TyliszczakA., 2002,Applicationofpreconditioningmethods for compressible flows, PhD thesis, PolitechnikaCzęstochowska/vonKármán Institute for Fluid Dynamics 15. Vichnevetsky R., Bowles J.B., 1982, Fourier analysis of numerical appro- ximations of hyperbolic equations, SIAM, Philadelphia Wpływ jawnej filtracji metodą kompaktową na wzrost zaburzeń w nieustalonym przepływie z warstwą ścinającą Streszczenie W pracy przedstawiono efektywność jawnej filtracji zjawisk wysokoczęstotliwo- ściowych przy zastosowaniu filtracji kompaktowej. Analizie poddano nieściśliwy dwu- wymiarowyprzepływ zwarstwą ścinającąw którym rozpatrywano ewolucję zaburzeń nałożonych na pole prędkości średniej. Przedstawiono porównania wartości amplitud zaburzeń obliczonych przy zastosowaniu różnych procedur filtracji oraz dla innych metod stosowanych do eliminowania zjawisk wysokoczestotliwościowych. Manuscript received November 18, 2002; accepted for print November 27, 2002