New explicit two-dimensional higher order filters

Computers & Fluids 39 (2010) 1848–1863

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

New explicit two-dimensional higher order filters Tapan K. Sengupta *, Yogesh G. Bhumkar Department of Aerospace Engineering, Indian Institute of Technology Kanpur, Kanpur 208016, India

a r t i c l e

i n f o

Article history: Received 15 December 2009 Received in revised form 4 April 2010 Accepted 16 June 2010 Available online 23 June 2010 Keywords: Two-dimensional filter Spatial filters Stabilization effects Rotary oscillation Accelerated flows Transitional flows q-Waves

a b s t r a c t In the present work, a new class of explicit two-dimensional filter is proposed that is significantly different from the conventional one-dimensional Padè type filters. Comparison of performance between oneand two-dimensional filters are made in studying propagation problems and actual viscous flow problems. In this preliminary work, the focus is kept solely on developing central filters with its real transfer function. Many applications related to Navier–Stokes equation have been shown to demonstrate that these two-dimensional filters help in numerical stabilization; control aliasing and control over spurious upstream propagating disturbances. Furthermore, its potential for large eddy simulation (LES) is suggested by band-limiting the solution via filters in post-processing mode. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Discrete computations introduce implicit filters – as shown in [12,14,28,29] using Fourier transform technique. On the other hand, quite often explicit filters are used to perform multiple functions in computing. Main purpose of such filtering has been to provide numerical stabilization at high wavenumbers [4–6,30] arising due to mesh non-uniformities and boundary conditions. Filtering also alleviate nonlinear instabilities and aliasing – as noted in [14,16]. Explicit filters used in a post-processing mode can bandlimit computed unknowns, to achieve similar results as is practiced in traditional LES by filtering the governing equation along with sub-grid scale models, described in [7,11,24,27] and other references contained therein. In recent times, explicit spatial filters have been used on the unknowns obtained by computing unfiltered governing equations without using sub-grid scale model as in MILES or ILES approach [9,10,22]. Bogey and Bailly [1] have investigated on the choice of quantity to be filtered explicitly in solving aeroacoustic problems. All the explicit filters in the literature are one-dimensional, even when applied for multi-dimensional problems. Here, we present explicit two-dimensional filters intended to perform both the roles discussed above, with applications shown for unsteady flows and instability studies. This is aided due to advances in understanding functions of filters based on spectral analyses reported in [4,6,13,17,22]. In [4,6], the transfer function of the

* Corresponding author. Tel.: +91 512 2597945; fax: +91 512 2597561. E-mail address: [email protected] (T.K. Sengupta). 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.06.014

central, one-dimensional filters have been obtained in the spectral plane explaining the control of high wavenumber instabilities by these Padè filters. Structure of Padè filters [4,31] makes use of easy solution of tridiagonal matrix equation for the filtered quantities, expressed as quantities with caret evaluated in terms of the unfiltered variables u on the right-hand side of the equation given by,

a1 u^j1 þ u^j þ a1 u^ jþ1 ¼

M X an ðujþn þ ujn Þ 2 n¼0

ð1Þ

with M as the stencil size defining the order of the filter used. For example, M = 1 provides a second order filter and successive increase of M by one increases the order of the filter by two. In Eq. (1), a1 is a free parameter that is restricted between 0.5 and +0.5. In addition to this filter in the interior, boundary filters were proposed in [5,6,31] and the transfer function for interior and near-boundary points reported separately without any coupling between these two sets of points. Effects of boundary filters on numerical instability were also obtained in isolation. The transfer function of the filters in the full domain can be provided by express^ j g ¼ ½Bfuj g and the nodal value obtained in the ing Eq. (1) as ½Afu spectral (k) plane following the methodology described in [13] by,

PN bjl eikðxl xj Þ T j ðkÞ ¼ Pl¼1 N ikðxl xj Þ l¼1 ajl e

ð2Þ

The role of Tj(k) in stabilizing a numerically unstable method can be explained by noting that the original amplification factor at the jth b j ðkÞj ¼ jGj ðkÞjjT j ðkÞj. One notices that when a node, Gj(k) alters to j G central filter is used, as in (1), then the role of Tj(k) is to reduce the unknown differentially for different values of k. It has also been

1849

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

noted in [6,22] that the intensity of filtering can be reduced in two ways for one-dimensional filters: either by increasing the order of the filter (i.e. increasing the value of M in (1)) or increasing the value of a1 towards its limiting value of 0.5. Often one-sided filters have been used in solving non-periodic problems, as developed in [6] and analyzed in [22]. However, one-sided filters also make the transfer function complex, thereby providing additional phase shift and causing numerical dispersion relation to change apart from changing dissipation also [22]. This is one of the serious aspects of filtering that can affect computations adversely. At the same time, proper analysis also allows one to use filters constructively. For example, for most of the numerical methods in imposing the numerical dispersion relation (that is different from the physical dispersion relation) show an extreme form of dispersion error for a range of kh, for which waves move spuriously in an opposite direction to that of the physical direction. This has been rigorously shown in [20] using one-dimensional convection equation. This has been predicted earlier in [29] for specific numerical methods in a heuristic fashion and these spurious waves have been called as q-waves. In this context, the authors in [22] have introduced a new spatial one-sided or upwind composite filter and it has been shown that such filters not only allow one to add a controlled amount of dissipation, but also change the basic dispersion relation of the numerical scheme in such a way that for a small range of Nc values, one prevents formation of q-waves. Despite the advantages of such well-designed filters, it is preferable to use central filters, even when very high order filters are

a

1 0.98

Transfer function

needed. While higher order filters have large stencil size and cannot be used at points close to the boundary, an alternative was proposed in [6] through the use of least order centered (LOC) filters. It is noted from Eq. (1) that one can use a second order filter at all the interior points and hence one can progressively keep increasing the order of the filter, as one moves inwards from the boundary. This is the essence of using LOC filters that does not additionally alter the dispersion relation beyond that is incorporated by the basic numerical method. There is another aspect of filters that has remained unattended so far. This is related to avoiding the directional artifact introduced in filtering when one-dimensional filters are applied sequentially in a multi-dimensional problem. This is attempted in the present work where we propose and apply two-dimensional filters for some typical two-dimensional signal and flow problems. To our knowledge, this has not been done before and the present work reports two-dimensional filters of different orders, while utilizing LOC concept. The present paper is formatted in the following manner. In the next section, we develop a new class of two-dimensional (2D) filters of different orders. These filters are tested on a standard test case of wave propagating obliquely with respect to a co-ordinate axis. In Section 3, we study few problems requiring the solution of two-dimensional Navier–Stokes equation, which cannot be solved without using filters. This will help one to test the efficacy of the proposed filters in comparison to one-dimensional (1D) filters. The paper closes with a summary and some conclusion.

kjhj=0.0

0.96

2 order filter, α2=0.2495 th 4 order filter, α2=0.2495 th 6 order filter, α2=0.2495 2nd order filter, α2=0.2 th 4 order filter, α2=0.2 6th order filter, α2=0.2 nd

0.94 0.92 0.9 0

0.5

1

1.5

2

2.5

3

2

2.5

3

kihi

b

1

Transfer function

0.8

kihi=kjhj

0.6

2nd order filter, α2=0.2495 th 4 order filter, α2=0.2495 6th order filter, α2=0.2495 2nd order filter, α2=0.2 th 4 order filter, α2=0.2 6th order filter, α2=0.2

0.4 0.2 0 0

0.5

1

1.5

kihi Fig. 1. Transfer function variation of 2D second, fourth and sixth order filter at the indicated values of filtering coefficient a2 plotted as a function of kihi along the lines: (a) kihi axis; (b) kjhj = kihi. (c) Transfer function contours in (kihi, kjhj)-plane for (i) a 1D fourth order filter with a1 = 0.20 and (ii) a 2D fourth order filter with a2 = 0.20 are shown. Frame (iii) shows a non-aliased region and aliased region in (k–k0 )-plane for product evaluation. Aliased region is shown with hatches. Note that 2D filter controls aliasing more effectively by attenuating aliased components near kihi = kjhj = p and does not attenuate low wavenumber components.

1850

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

side giving rise to a penta-diagonal matrix. To solve this linear algebraic equation, one must have the diagonal dominance of the matrix on the left-hand side of Eq. (3) and that is ensured by ja2j 6 1/4. Transfer function of the filter is defined in the k-plane, obtained by taking Fourier transform of both sides in (3). Thus, the transfer function of the 2D second order filter for the (i, j)th node is given by,

2. Two-dimensional filters Despite discussed advantages, 1D filters also show strong directional behavior of the solution, as noted in [21,22] for two-dimensional fluid flow problems. Filtering in azimuthal or wall-normal direction causes unphysical vortex smearing in the respective direction. This is proposed to be removed with the help of twodimensional filters as much as possible and is described below. Stencil for the two-dimensional filters is proposed to have a general form,

^ i1;j þ u ^iþ1;j þ u ^ i;j1 þ u ^ i;jþ1 Þ ¼ ^ i;j þ a2 ðu u

M X n¼0

an ðuin;j þ ui;jn Þ 2

T ij ðki hi ; kj hj Þ ¼

a0 þ 2a1 ¼ 1 þ 4a2

3

a0 ¼ 2a1

ð6Þ

obtained from Tij(kihi = kjhj = p) = 0. The reason for this normalization of the filter is discussed later.

(ii)

2D 4th order filter, α2 = 0.20 0.05 0.2

3

0.05 0.2

0.

0.

kjhj

kjhj

0.8 0.9 0.95

1

0.

2

0.7

1

0. 99

0. 9

0.5

7

0.5

2

ð5Þ

For the second relation, we set the condition on transfer function Tij to attenuate completely at the Nyquist limit (kihi = kjhj = p), which gives the relation as,

1D 4th order filter, α1 = 0.20

(i)

9

0.

8

95

9

0.999

0.999

0 0

1

kihi

2

ð4Þ

In the above expression hi and hj are the uniform grid spacings in the physical plane in i and j directions, respectively. In order to obtain the coefficients a0 and a1 in terms of the filtering parameter a2, one equates the Taylor series expansion on both sides of the Eq. (3) about the (i, j)th point. The first relation ^ i;j , that gives is obtained by equating the coefficients of ui,j with u the consistency condition as,

ð3Þ

where M is the order of the filter and is equal to one for a 2D second order filter and every increase of M by one increases the order of the filter by two. For the second order, 2D filter, i and j varies from 2 to (Ni  1) and 2 to (Nj  1), respectively, in i and j directions. For the fourth order filter, M = 2 and Eq. (3) applies from 3 to (Ni  2) and 3 to (Nj  2), respectively. Similarly, higher order filters can be defined. Development of this filter is helped by the analysis method of 1D central filters given by (1) and its upwind extension described ^ , on the left-hand side are in [22]. In Eq. (3), filtered quantities u evaluated in terms of the unfiltered variables u, on the right-hand side. Filtering coefficient a2 and other coefficients an’s on the righthand side, determine the transfer function. Application of the filter amounts to solving a linear algebraic equation with the left-hand

c

ao þ a1 ½cosðki hi Þ þ cosðkj hj Þ 1 þ 2a2 ½cosðki hi Þ þ cosðkj hj Þ

3

0 0

1

kihi

k km Aliased region

(iii)

Non-aliased region km k’ Fig. 1 (continued)

2

3

1851

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

Unfiltered solution

5

6

0.004 0.002

4

t=7 1E-12

0

t=0

3

0

1

2

0.004

3

x

4

5

-0.2

t=0 0.002

3

0

6

d

t=7 3.98403

0 11.3371

0

1

2

3

x

4

5

6

0.004

t=7

0.002

t=7

y

t=0

0

1E-12

0.002

0

3 2

6

6th order composite 2D filter α2=0.24925

4

y

0

t=7 Max : 0.967 Min : -0.663

5

t=7

y

4

y

0

1

2nd order 1D filter α1=0.4985

5 0.002

t=0

0.002

3 2

t=7 Max : 12.701 Min : -3.411

1 0

4

t=7

2 t=7 Max 0.970 Min :-0.664

1

6

0.002

0

2

b

5

t=7 -0.1

0.002

y

y

0

0

0.004

t=7

y

6

th

6 order composite 1D filter α1=0.4985

c

y

a

0

1

2

3

x

4

5

t=7 Max : 0.968 Min : -0.663

1 6

0

0

1

2

3

x

4

5

6

Fig. 2. Propagation of a wave-packet moving along a line inclined at angle 45° to the co-ordinate axes shown when (a) the solution is not filtered; (b) solution is filtered after every time step using 1D second order filter in both the directions; (c) solution is filtered after every time step using 1D composite sixth order filter in both the directions; (d) solution is filtered after every time step using 2D sixth order composite filter.

Solving Eqs. (5) and (6), one gets a0 and a1 in terms of a2 as,

1 a0 ¼ þ 2a2 2 1 a1 ¼ þ a2 4

ð7Þ

For the fourth order filter, one has M = 2 and the transfer function is obtained as, T ij ðki hi ; kj hj Þ ¼

ao þ a1 ½cosðki hi Þ þ cosðkj hj Þ þ a2 ½cosð2ki hi Þ þ cosð2kj hj Þ 1 þ 2a2 ½cosðki hi Þ þ cosðkj hj Þ ð8Þ

And here the consistency condition is given by,

a0 þ 2a1 þ 2a2 ¼ 1 þ 4a2

ð9Þ

To evaluate three unknown an’s, additionally we equate the coefficients of the second derivatives on either side of (3) that gives,

a1 þ 4a2 ¼ 2a2

ð10Þ

Above two equations are augmented by the vanishing condition of the transfer function (8) at the Nyquist limit:

a0  2a1 þ 2a2 ¼ 0

ð11Þ

Solution of (9)–(11) yields the following

a0 ¼

ð5 þ 12a2 Þ ; 8

a1 ¼

1 þ 4a2 ; 4

a2 ¼

1 þ 4a2 16

ð12Þ

A similar exercise for the sixth order filter yields the following coefficients for the right-hand side of (3) as,

20a2 þ 11 68a2 þ 15 ; a1 ¼ ; 16 64 12a2  3 1  4a2 ¼ ; and a3 ¼ 32 64

a0 ¼

a2 ð13Þ

Variations of Tij in the spectral (kihi, kjhj)-plane are plotted in Fig. 1a and b, along specific directions for two values of a2. Two-dimensional filters proposed in (3) are directional and to test it, we plot the transfer functions for the second, fourth and sixth order filters along the kihi-axis and along the kihi = kjhj line. Note that the value of the transfer function along kjhj-axis will be same as obtained for the variation along the kihi-axis. Results are plotted as a function of kihi, from zero to its Nyquist limit in these figures. The transfer function is purely real and in Fig. 1a results are shown for a2 = 0.2495 and 0.20 along the kihi-axis. The transfer function for a2 = 0.2495 varies very little between 1.0 and 0.999, with some difference among different ordered filters noted in the fourth decimal place for intermediate wavenumbers only. Whereas for a2 = 0.20, there is significant variation of Tij among different ordered filters, where it varies from 1 to 0.90 along the kihi-axis, again only in the intermediate range of wavenumbers. In contrast, along the kihi = kjhj line, variation of Tij is from a value of one (at kihi = kjhj = 0) to zero at the Nyquist limit of (kihi = kjhj = p) for both values of a2. However, one cannot

1852

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

(a) Experimental flow visualisation [25] Re = 150, A= 2.0, ff/fo = 1.5

(c) 4th order 2D filter α2 = 0.2495 t = 49.06

t = 49.71

t = 50.36

nd

(b) 2 order 2D filter α2 = 0.2495

th

(d) 6 order 2D filter α2 = 0.2495

t = 49.06

t = 49.06

t = 49.71

t = 49.71

t = 50.36

t = 50.36

Fig. 3. Flow past a circular cylinder executing rotary oscillation for Re = 150; A = 2 and ff/f0 = 1.5. (a) Experimental visualization pictures (top left) from [25] are compared with computed results using: (b) second order 2D filter; (c) fourth order 2D filter and (d) sixth order 2D filter. These cases are computed using the indicated filter coefficient a2. The negative vorticity contours are shown by dotted lines and positive vorticity contours by solid lines.

Table 1 Quantitative comparison of the normalized drag coefficient. Cases A = 2, A = 2, A = 2, A = 2,

ff/f0 = 1.5 ff/f0 = 2.0 ff/f0 = 3.0 ff/f0 = 4.0

Experimental value of CD/CD0 [25]

Numerical value of CD/CD0

1.171 1.024 0.878 0.902

1.245 1.010 0.855 0.864

distinguish among the three different ordered cases plotted for a2 = 0.2495 and all of them drop from one to zero for kihi P 2.5. For a2 = 0.20, the transfer functions are distinct for different filters with different orders. The advantage of using these two-dimensional filters given by Eq. (3), over conventional one-dimensional filters is explained with

the help of Fig. 1c. Frame (i) shows variation of transfer function in the (kihi, kjhj)-plane for one-dimensional fourth order filter with a1 = 0.20 applied sequentially. The transfer function variation for a two-dimensional fourth order filter with a2 = 0.20 is shown in frame (ii) of Fig. 1c. The restriction on transfer function to attain zero value for kihi = p or kjhj = p, imposed for one-dimensional filter coefficient leads to excessive filtering of the solution as shown in frame (i) for high values of either kihi or kjhj. Transfer function attains a value of 0.95 at kihi ’ 1.2 in frame (i), while same value of transfer function is attained at kihi ’ 2.15 in frame (ii), for the two-dimensional fourth order filter with a2 = 0.20. This shows that there is a unphysical attenuation of low wavenumber components for a one-dimensional filter as compared to two-dimensional filters. An ideal filtering procedure should only remove high wavenumber components from the solution which are responsible for numerical instabilities and aliasing. Thus, one-dimensional filtering

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

15

1853

Unfiltered case nd 1D 2 order filter, α1=0.4990 nd 2D 2 order filter, α2=0.2495 th 2D 4 order filter, α2=0.2450 th 2D 6 order filter, α2=0.19

Re = 150, A= 5, ff /f0 = 0.5

12

Cd

9

6

3

0

20

30

40

50

60

time Fig. 4. Computed drag coefficient variation with time showing unfiltered and filtered solutions for the case of flow past a circular cylinder executing rotary oscillation for Re = 150; A = 5 and ff/f0 = 0.5. The unfiltered solution blows up after t = 47. Filter type and filter coefficients are as indicated in the legend.

is inefficient from this perspective as compared to two-dimensional filtering. In an actual governing equation involving product terms, one-dimensional filter removes physically important low wavenumbers – which the two-dimensional filters will not do as per the plotted transfer function shown in frame (ii) of Fig. 1c. The restriction on transfer function chosen for the two-dimensional filter to attain zero value at kihi = kjhj = p, imposed while deriving the filter coefficient has two important numerical consequences. The first one is to avoid excessive filtering which is a problem with conventional one-dimensional filter in evaluating product terms. The second aspect is related to aliasing of numerical solution. Aliasing error results from evaluation of product terms in a finite grid. Aliasing error can either be present while numerically evaluating nonlinear product terms or even while evaluating linear terms in the transformed plane where products have to be evaluated due to grid transformation terms appearing in transformed plane equations [14]. Calculations performed with a high resolution discretization schemes such as OUCS3 scheme [13], are more susceptible to aliasing error as compared to conventional low resolution CD2 spatial discretization schemes, as explained in [14] and briefly explained with frame (iii) of Fig. 1c. Consider a product term f(n)x(n) to be evaluated. The Fourier spectral representations of these two terms f(n) and x(n) are given as,

f ðnÞ ¼

Z

km

km

0

0

0

Fðk Þeik n dk ;

and xðnÞ ¼

Z

km

XðkÞeikn dk

km

where km is maximum resolved wavenumber for a given grid. Aliasing occurs when the term (k + k0 ) is greater than km [14] as shown hatched in frame (iii) of Fig. 1c in the spectral plane. The high wavenumber components near kihi = kjhj = p are folded back to resolved wavenumbers. The maximum problem occurs when both the Fourier amplitudes are near to the Nyquist limit. Two-dimensional filters automatically help in dealiasing the solution due to the nature of the transfer function shown in frame (ii) of Fig. 1c. We note specifically that the filtering does not attenuate the solution by a large amount when either kihi or kjhj has small values as is the case of using one-dimensional filters applied sequentially so that the transfer function is as shown in frame (i) of Fig. 1c. Furthermore the

amount of filtering near kihi = p or kjhj = p can be explicitly controlled via the choice of a2. 2.1. Implementation procedure for the two-dimensional filter The left-hand side of the stencil proposed for the two-dimensional filters by Eq. (3) represents a Laplacian operator. Thus, while applying the filter, the boundary values of the domain remain unchanged and the filter is applied only for interior points. When LOC approach [6] is used, no near-boundary filters are needed. Problems of numerical instabilities near inflow and excessive damping near outflow, while using high order filters due to one-sided stencils for near-boundary nodes are avoided here by the LOC approach. To solve Navier–Stokes equation, we have used the stream function–vorticity formulation. Success of the method depends upon capturing the vorticity dynamics correctly by solving accurately the vorticity transport equation (VTE), as given in Eq. (18). Vorticity is continuously generated at the body due to no-slip condition. Numerical instabilities and aliasing problems mostly arise due to adverse numerical properties at high wavenumbers for velocity and vorticity fields. These can be controlled effectively only when the problematic higher wavenumbers acting together are removed by explicit filtering. The 2D filter is applied by following the steps given below for VTE. (1) Consider that the numerical solution is available in the domain at the nth time step. (2) The right-hand side of Eq. (3) is calculated from the known vorticity field for all the points inside the domain. (3) Eq. (3) is solved in an iterative manner, using the Bi-CGSTAB variant of the conjugate gradient method [26] to obtain the filtered vorticity data. (4) As the maximum modulus of the residue of Eq. (3) attains a value less than 106 in the domain, filtered vorticity field obtained is taken as converged. Otherwise this operation is repeated. (5) The filtered vorticity field is used to time advance the unknown vorticity field at (n + 1)th time step by solving VTE (Eq. (18)).

1854

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

t = 47

Flow direction

Re=150, A=5, ff/f0=0.5 Unfiltered case t = 47

Re=150, A=5, ff/f0=0.5 nd

1D 2 order filter α1=0.4990 t = 47

Re=150, A=5, ff/f0=0.5 nd

2D 2 order filter α2=0.2495 Fig. 5. Flow field past a circular cylinder executing rotary oscillation for Re = 150; A = 5 and ff/f0 = 0.5 at t = 47. (a) Computed vorticity contour plot for unfiltered case; (b) computed vorticity contour plot for the second order 1D filter with a1 = 0.499; (c) computed vorticity contour plot for the second order 2D filter with a2 = 0.2495.

Depending on the level of filtering needed to control high wavenumber components, above sequence is either repeated at every time step or after user-specified number of time steps. Thus, along with the choice of the filter coefficient a2 and the order of the filter, frequency of filtering is an additional parameter that can be used to control filter performance. As we would also show that one does not need to filter the variables over the complete domain. This actually helps in bringing down the computational cost enormously – a quantitative estimate of it is also provided later. 2.2. Performance comparison between 1D and 2D filters To study this aspect, we have considered the propagation of a wave-packet moving along a line that is inclined at h = 45° to the co-ordinate axes. The wave-packet is initially located at (xo = 4, yo = 4) and whose functional form is given by,

f ¼ e500½ðxxo Þ

2

þðyyo Þ2 

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos 100 ðx  xo Þ2 þ ðy  yo Þ2

ð14Þ

Thus the solution is sought for the governing equation

@u @u @u þ cx þ cy ¼0 @t @x @y

ð15Þ

where cx = c cosh and cy = c sinh and we have chosen cx = cy = 0.1. Time step is fixed by choosing Ncx = 0.1, where N cx ¼ cx Dt=Dx. The domain considered is 0 6 (x, y) 6 6, with 1201 equi-spaced points. We have shown the results of calculations, with and without the filters, in Fig. 2. The top left frame (a) shows unfiltered solution at t = 0 and 7. In frame (b), results shown are obtained using a second order 1D filter applied after every time step in both x and y directions. In this case, value of a1 is chosen as 0.4985 and for the solution at t = 7, one notices accumulation of error near the origin. The spurious numerical solution in the vicinity of the origin is shown in an inset within the frame and one can note very large contour values in Fig. 2b. While the wave-packet propagates along h = 45°, very high resolved wavenumber components travel upstream, for which the numerical group velocity is negative for the chosen spatial and temporal discretization schemes (OUCS3 and RK4 schemes, respectively).

1855

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

1

Transfer function

0.8

0.6 nd

1D 2 order filter, nd 2D 2 order filter, nd 2D 2 order filter, 2D 2nd order filter,

0.4

0.2

α1=0.4990 0 α2=0.2495, θ = 45 0 α2=0.2495, θ = 30 α2=0.2495, θ = 150

0 0

0.5

1

1.5

2

2.5

3

kh Fig. 6. Comparison of transfer function between 1D and 2D filters for the indicated filter coefficients. For the 2D filter, transfer function varies with h-angle made by a line with either kihi or kjhj axes. For a1 = 2a2, transfer functions are identical for 1D and 2D filter along h = 45°.

a t = 1.92

t = 1.92

4 order 2D composite filter, α2=0.249 th

t = 1.92

2 order 2D filter, α2=0.249 nd

t = 1.92

6 order 2D composite filter, α2=0.249 th

Fig. 7. (a) Flow past NACA 0015 aerofoil at 90° angle of attack for the accelerated oncoming flow for Re = 5200 (as defined in [2]). Computed results are compared with experimental results [2] shown in top left frame at t = 1.92. Computed results are obtained using second, fourth and sixth order filters with indicated filter coefficients. Negative vorticity contours are shown by dotted lines and positive vorticity contours by solid lines. (b) Flow past NACA 0015 aerofoil at 90° angle of attack for the accelerated oncoming flow for Re = 5200 (as defined in [2]). Computed results are compared with experimental results [2] shown in top left frame at t = 2.42. Computed results are obtained using second, fourth and sixth order filters with indicated filter coefficients. Negative vorticity contours are shown by dotted lines and positive vorticity contours by solid lines.

1856

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

b

t = 2.42

t = 2.42

4th order 2D composite filter, α2=0.249

t = 2.42

t = 2.42

2nd order 2D filter, α2=0.249

6th order 2D composite filter, α2=0.249 Fig. 7 (continued)

The numerical properties of the chosen discretization method are given in [20,22]. These upstream propagating waves reflect from the computational domain boundary and contaminate the solution strongly. It is important to explain the reason as to why the unfiltered solution did not produce any instability, while the 1D filter with the chosen value of a1 causes an instability. As shown in [20], the error dynamics is governed by the following propagation equation,

h @e @e cN i @uN þc ¼ c 1  @t @x c @x  Z Z  0 V gN  cN 0 0  ik A0 ½jGjt=Dt eik ðxcN tÞ dk dk k Z LnjGj A0 ½jGjt=Dt eikðxcN tÞ dk  Dt

cantly less amount of accumulated error at the origin at t = 7, as seen in Fig. 2c. This is due to the fact that the sixth order 1D filter modifies the numerical amplification factor less, as compared to the second order 1D filter. Finally, results for another case are reported where a 2D sixth order composite filter has been used with a2 = 0.24925. This value of filter coefficient for the 2D filter corresponds identically to the chosen sixth order 1D filter coefficient. Here, the solution contours near the origin, shown in the inset of Fig. 2d indicate negligibly small value of error (of the order of 1012). This establishes utility and superiority of 2D filters over the 1D filters in directional propagation of signals.

3. Solutions of Navier–Stokes equation

ð16Þ

It is evident that the forcing of error is contributed by three sources on the right-hand side. If one chooses a neutrally stable scheme, then jGj = 1 and the last term on the right-hand side does not contribute to error. For the unfiltered case, the chosen numerical parameters give rise to exact neutral stability. However, when the solution is filtered [22], modified numerical amplification facb ¼ jGjjT ij jÞ does not remain neutral anymore. This seed of ertor ðj Gj ror at high wavenumbers causes upstream propagating waves that is evident in Fig. 2b. Next, when the second order 1D filter is replaced by a composite sixth order 1D filter [22], we have signifi-

Applications of two-dimensional filters are shown here by computing three different unsteady flow cases: (a) flow past a cylinder executing rotary oscillation [18,19,25]; (b) accelerated flow past a NACA-0015 aerofoil at 90° angle of attack [2,21] and (c) transitional flow past a SHM-1 Honda aerofoil [3]. These flows are solved in stream function (w)–vorticity (x) formulation, that ensures solenoidality of velocity and vorticity field simultaneously. The Navier–Stokes equations in (w–x) formulation are given by,

    @ h22 @w @ h11 @w ¼ h11 h22 x þ @n h11 @n @ g h22 @ g

ð17Þ

1857

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

10000 105

30000 25000

Amplitude

α2=0.249

8000 10

= filter 1.92 2nd order t2D 4th order 2D filter 6th order 2D filter

6000 103

time = 1.92

4

2

10 4000

20000

101 2000

15000

100

0

2nd order 2D filter 4th order 2D composite filter 6th order 2D composite filter 0

0.2 0.1

0.4 0.2

0.6 0.3

0.8 0.4

0.5 1

10000 5000 0

1

2

3

kh 105

30000

10

t = 2.42

103

25000

Amplitude

α2=0.249

4

102

20000

2nd order 2D filter 4th order 2D composite filter 6th order 2D composite filter

101 100

15000

0

0.1

0.2

0.3

0.4

0.5

10000 5000 0

1

2

3

kh Fig. 8. Comparison of FFT of vorticity data at t = 1.92 and t = 2.42, for the computed cases using second, fourth and sixth 2D filters with indicated filter coefficient.

@x @x @x þ h22 u þ h11 v @t @n @g      1 @ h22 @ x @ h11 @ x ¼ þ Re @n h11 @n @ g h22 @ g

hðtÞ ¼ h0 cosðxf tÞ

h11 h22

ð18Þ

where h11 and h22 are the scale factors used in mapping the physical (x, y)-plane to a computational (n, g)-plane [21–23]. The vorticity transport equation is time-advanced using RK4 scheme while the Laplacian operator is discretized in self-adjoint form by second order central differencing and the convection terms are discretized by OUCS3 scheme [13], for all the numerical results presented in the following. 3.1. Flow past a cylinder executing rotary oscillation Rotary oscillation of cylinder has been studied by many researchers for the purpose of understanding basic fluid flow past a bluff body with the rotary oscillation giving rise to drag reduction [8,18,19,25]. All these studies except [19] are for low Reynolds numbers, for which the unsteady two-dimensional flow is laminar. In the present study also, the flow is simulated for Re = 150, for which good quality experimental results are available in [25]. The rotary oscillation of the cylinder is defined by,

ð19Þ

where xf = 2pff and A = h0xfd/(2U1), define the circular frequency and amplitude of the imposed rotary oscillation. In [25], experimental results are provided for Re ¼ U 1m d ¼ 150 with different amplitudes (A) and forcing frequencies (ff) of oscillation. Here d is the diameter of the cylinder, U1 is the oncoming free stream velocity and m is the kinematic viscosity of the ambient fluid. Here, the forcing frequency ff is normalized by f0, where f0 is the shedding frequency of the stationary cylinder at same Re. Briefly, the rotary oscillation of the cylinder causes unsteady separation of the flow from the cylinder surface and the resulting small scale vortices interrupt the formation of regular Karman– Bènard vortices in the wake. Thus, the imposed articulation interferes with the natural shedding. This interference, termed as aerodynamic tripping [18], creates small-scale vortices at the cylinder surface. This leads to reduced wake-width, and thereby brings early drag crisis and drag reduction. In this flow, large range of wavenumbers are created by surface excitation, as compared to the stationary cylinder case at the same Reynolds number. This flow therefore provides a case to test the efficacy of a filter in capturing small-scale physical vortices, while eliminating spurious higher wavenumber errors. In [19], the compact scheme OUCS3 was used to solve this problem at the higher Reynolds number of

1858

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

a Flow Direction

t = 5.0

th

4 order 1D composite filter α1 = 0.49

6

Re = 10.3X10 Grid Size = 511X300 Angle of attack = 00

t = 5.10

t = 5.20

t = 5.30

Fig. 9. (a) Computed vorticity contours of flow past SHM-1 aerofoil [3] at 0° angle of attack for Re = 10.3  106 shown at the indicated time instants, obtained using 1D fourth order composite filter with a1 = 0.49. (b) Computed vorticity contours of flow past SHM-1 aerofoil [3] at 0° angle of attack for Re = 10.3  106 shown at the indicated time instants, obtained using 2D fourth order composite filter with a2 = 0.245.

Re = 15,000, to control the flow using a strategy encoded by genetic algorithm. In Fig. 3, computed vorticity contours are compared for the case of A = 2, ff/f0 = 1.5 and Re = 150, with the corresponding experimental visualization pictures [25] to establish the correctness of the numerical method used with the two-dimensional filters. This is one of the low amplitude cases reported [25] that can also be computed without the need of any filter. A (150  450) grid with 150 equi-angular points in the azimuthal direction and 450 points in the radial direction has been used. Wall resolution in the radial direction is fine with the second point 0.00095d distance away from the cylinder surface and the calculations are performed with Dt = 5  105. Results are shown for almost a full cycle of oscillation in Fig. 3, with frames shown at a non-dimensional time interval of 0.65. Computational results are shown here using three filters whose transfer functions were shown in Fig. 1a for a2 = 0.2495. In this case filters are applied in the full domain for the vorticity data. Despite very large value of a2, one notices distinct differences between the results with second order filter and the results obtained using fourth and sixth order filters. The latter

computed cases exhibit excellent match with the experimental results. The second order filter causes the shed vortices to diffuse and align wider with respect to center-line, while the fourth and sixth order filter cases show correct shape and alignment of the shed vortices. Thus, this exercise shows that for this flow, the second order filter has larger influence with possible alteration of the flow field, while the fourth and sixth order filters do not alter the flow field for this value of a2 and play the roles of numerical stabilization and removal of aliasing error. The qualitative comparison of vorticity contours with the experimental flow visualization pictures shown in Fig. 3 highlights the need for the proper choice of filtering coefficient a2 along with the order of the two-dimensional filter. This will result in least alteration of the flow field while numerically stabilizing the flow and removing aliasing error. In addition to this qualitative comparison, we have obtained a quantitative comparison for the cases of Re = 150, A = 2 and different forcing frequencies ff/f0 = 1.5,2,3,4, as shown in Table 1. All the computations are performed by applying a fourth order two-dimensional filter with a filtering coefficient a2 = 0.2495. Mean drag coefficient CD of a particular rotary

1859

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

b

Flow Direction

4th order 2D composite filter α2 = 0.245

t = 5.0

6

Re = 10.3X10 Grid Size = 511X300 Angle of attack = 00

t = 5.10

t = 5.20

t = 5.30

Fig. 9 (continued)

oscillation case is normalized by mean drag coefficient CD0 of a stationary cylinder at the same Reynolds number. Comparison of this normalized drag coefficient obtained experimentally [25] and numerically shows a good quantitative comparison. In Fig. 4, drag coefficient variations with time are shown for a higher amplitude case of Re = 150, A = 5 and ff/f0 = 0.5. A grid with (150  450) points has been taken with the outer boundary located at a distance of 20d from the surface of the cylinder. These computed cases display high mean drag coefficient with large excursion about the mean with time. For the unfiltered case, the solution blows up a little before t = 48. In this figure, effect of different order filters on numerical stabilization property is studied. In one of the filtered cases, conventional 1D second order filter [4] is used in the azimuthal and radial direction over the full computational domain taking a1 = 0.499. When a sixth order composite filter was employed for the same filter coefficient, the numerical solution blew up (not shown). In the other three displayed filtered cases, we have applied the proposed two-dimensional second, fourth and sixth order filters, respectively, with different values of a2. Filter coefficient value has to be reduced with the increase in order of the 2D filters for stable numerical solution. From the figure, 1D and 2D second order filters with corresponding filter coef-

ficients have computed solutions without showing numerical instability for the chosen a1 and a2 indicated in the figure. For fourth order filter, value of a2 was reduced to 0.2450 and to 0.19 for the sixth order filter to compute without numerical instabilities. The qualitative changes in the flow after t = 50 result in different levels of drag coefficient variation. Effectiveness of any filter can be furthermore ascertained by looking at the detailed flow structures. This is shown in Fig. 5, where results obtained with second order 1D and 2D filters are compared along with the unfiltered solution case. The difference between 1D and 2D filters can be explained by looking at the instantaneous vorticity contours for the above case plotted at t = 47, in this figure. Vorticity contours for unfiltered case show presence of large coherent vortices in the wake with sharp vorticity gradient. Experimental data given in Fig. 10 of [25] for this case also indicate the presence of a coherent vortex in the near-wake with high spatial gradient of vorticity. Unfiltered solution blew up due to aliasing associated with this large gradient that creates a large bandwidth disturbance field. For the case of filtered numerical solution (by applying second order 1D filter after every time step), the stabilization results in weakening of the core of this vortex. Filtering is performed in both the radial and the azimuthal

1860

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

t = 5.0 2D filter α2=0.245 1D filter α1=0.490

105

Amplitude

104

1D 4th order filter 2D 4th order filter

3

10

2

10

1

10

0

0.5

1

1.5

2

2.5

3

2

2.5

3

2

2.5

3

kh 5

10

t = 5.1

4

Amplitude

10

3

10

102 101 100 0

0.5

1

1.5

kh 5

10

t = 5.2

4

Amplitude

10

103 2

10

101

0

0.5

1

1.5

kh Fig. 10. Comparison of FFT of vorticity data for flow past SHM-1 aerofoil [3] obtained using 1D fourth order composite filter with a1 = 0.49 and 2D fourth order composite filter with a2 = 0.245, respectively.

directions over the full domain. In addition to attenuation of the vorticity values at the core, vorticity also is smeared or diffused as compared to the unfiltered case. However, for a 2D second order filter, attenuation of the vorticity value at the core and the diffusion of vorticity is significantly less as compared to the 1D filter case. Difference between 1D and 2D filters for the same flow field with equivalent filter coefficients is explained from the transfer function shown in Fig. 6. In Fig. 6, transfer functions for the second order 1D and 2D filters are provided. For a 2D filter, h corresponds to the angle made with kihi axis. As we have chosen a1 = 2a2, transfer functions for 1D and 2D second order filters for h = 45° are identical. However for a 2D filter, transfer function plotted at other angles (h = 15°, 30°) are significantly higher as compared to the h = 45° angle case, across whole wavenumber range. Difference between these transfer functions increases sharply near kihi = p. Variation of transfer function in different directions results in less attenuated and diffused vortices in the wake for the 2D filter as compared to 1D filter. Thus, application of 2D filter is advantageous

over 1D filter in not altering the flow field. The filtering action along h = 45° can be further reduced by increasing a2. 3.2. Accelerated flow past a NACA-0015 aerofoil This problem is selected to identify the order of the most effective 2D filter by filtering the solution of Navier–Stokes equation for a difficult unsteady flow. This flow cannot be computed without a filter [21]. Here, a NACA-0015 aerofoil is kept at 90° angle of attack in the oncoming accelerated flow. Freymuth [2] has studied uniform accelerating flow past NACA-0015 aerofoil by flow visualization technique for different angles of attack and Reynolds numbers. In [2], Reynolds number was defined in terms of the uniform acceleration and the chord of the aerofoil with visualization made for constant Reynolds number cases. In the convention adopted here by defining Reynolds number based on instantaneous convection speed, the studied flow in [2] has linearly varying Reynolds number with time. Such time-dependent flows continually change their energy spectrum. Out of all the cases considered in [2], here we

1861

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

(i)

Time requirement for Streamfunction Equation

C.P.U.time(s)

0.4 0.3 0.2 Mean value = 0.022 s

0.1 0 5

5.0005

5.001

5.0015

5.002

5.0025

Non-dimensional time of Navier-Stoke sequation solver

(ii)

Time requirement per non-dimensional time step of Navier-Stokes equation

C.P.U.time(s)

0.4 0.3 0.2 Mean value = 0.2815 s

0.1 0

5

5.0005

5.001

5.0015

5.002

5.0025

Non-dimensional time of Navier-Stoke sequation solver Time requirement for 2D 4th order filter

(iii) C.P.U.time(s)

0.4 0.3 0.2 Mean value = 0.0129 s

0.1 0

5

5.0005

5.001

5.0015

5.002

5.0025

Non-dimensional time of Navier-Stoke sequation solver Fig. 11. Various time requirements for computing (i) stream function equation; (ii) one step of time increment of Navier–Stokes equation and (iii) time required for 2D filtering operation.

have computed the flow for Re = 5200 (in the notation of [2]) and an angle of attack of 90°. Various cases of accelerated flows past aerofoils at large angles of attack have been reported in [21], where other computational details can be found. The governing two-dimensional Navier–Stokes equations, Eqs. (17) and (18), are solved on a numerically generated orthogonal O-grid [14], with 257 points in azimuthal direction and 301 points in the wall-normal direction. After every time step, vorticity values for the first 20 g = constant lines are filtered using three 2D filters of different orders keeping a2 same (=0.249). For second order filter, 2D filtering is straightforward, while LOC or composite filtering [6,22] is required for fourth and sixth order cases. For this flow, vorticity generated on the aerofoil surface convects downstream in the wake and to control numerical instability, it is established here that one needs to apply 2D filters in a region very close to the body only. In Fig. 7a and b, experimental visualization pictures [2] at t = 1.92 and 2.42 are compared with computed vorticity contours obtained using three 2D filters. At these time instants, we have

good match between the flow structure obtained via visualization and the computed vorticity contours. Asymmetry in vortex structures observed in the flow visualization pictures is also noted in the numerical solution. Wavy appearance of the link connecting the vorticity on the surface of the aerofoil to the vortex core is observed in the experimental and the numerical solution obtained using the fourth order filter only. The use of 2D filter allows one to compute the flow field without excessive attenuation of flow structures and imposed directionality of computed flow structures, as obtained using 1D filters in [21]. Out of three 2D filters used for this example, the fourth order 2D filter is able to capture the flow structures with finer details. For example, undulations in the lower vortex structure at t = 1.92 and in the top vortex structure at t = 2.42 are seen clearly for the fourth order 2D filter case – as seen in the visualization picture. Why the fourth order filter is better than second or sixth order filters is explained with the help of FFT plots shown in Fig. 8 at t = 1.92 and t = 2.42. The FFT plots are for the filtered vorticity data on a line close to the aerofoil obtained using three 2D filters. Inset in the top frame for t = 1.92

1862

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

shows that the Fourth order filter is able to capture the flow structures, in the wavenumber range from kh = 0.15 to 0.5, better as compared to the second and sixth order filters. Similar observations can be made for the data shown in Fig. 8 at t = 2.42. The results of rotary oscillation cases shown in Figs. 3–6, established that fourth and sixth order filters are superior to second order filter. For the accelerated flow case shown in Figs. 7 and 8, one notes further that the fourth order filter is superior to the sixth order filter. This helps one to conclude that there is no need to seek filters of order higher than fourth order filter for higher accuracy. 3.3. Equilibrium flow past a SHM-1 Honda aerofoil Transitional flow past natural laminar flow (NLF) aerofoil at a high Reynolds number is studied here by applying the 2D fourth order filter to the solution of Navier–Stokes equation after every time step. The flow experiences bypass transition and the flow cannot be computed without an appropriate filter. Here, an orthogonal grid [14] around the SHM-1 aerofoil is generated for which sectional aerodynamic data are available in [3]. The grid consists of 511 points in the azimuthal direction and 300 points in the wallnormal direction, with the outer boundary located at ten chord distance from the aerofoil surface. The aerofoil is held at zero degree angle of attack and the Reynolds number based on the chord and free stream speed is 10.3 million, for which the flow is laminar up to the mid-chord [3]. Beyond that the flow downstream experiences physical instability and the purpose of present computations is to calculate it without the use of any models for transition and/or turbulence. Here, we have computed the flow field using 1D fourth order composite filter with a1 = 0.49 and a 2D composite fourth order filter with a2 = 0.245. Both the filters are applied to first 20 glines only for the vorticity data. This is the least number of lines which require filtering to avoid numerical instability. In Fig. 9a and b, we have plotted vorticity contours at the indicated times obtained using 1D and 2D fourth order filters. The top frames in these figures show the flow field over the complete aerofoil, while the bottom frames show zoomed view of the flow near trailing edge only where the flow experiences unsteady separation. Aft of the leading edge of the aerofoil, flow accelerates due to favorable pressure gradient and as it passes over the maximum thickness point, it experiences adverse pressure gradient. This adverse pressure gradient leads to formation of unsteady separation bubbles on the aerofoil surface. Presence of high wavenumbers in the solution leads to spurious unphysical upstream propagating disturbances [22] that must be avoided by using filters. In Fig. 9a and b, one can observe formation of small scale bubbles on aft portion of both the upper and the lower surfaces of the aerofoil. As these vortices move downstream, these cause secondary vortex-induced instabilities that results in formation of newer bubbles on the aerofoil surface. Use of 1D and 2D filters removes high wavenumber components that can cause aliasing and unphysical spurious upstream propagating q-waves [22]. In Fig. 10, FFT of vorticity data are plotted shown corresponding to the cases in Fig. 9a and b. These figures show that filtering action due to 1D and 2D filter starts showing differences from a non-dimensional wavenumber kh = 1.5 upwards, with 1D filter showing finite value of the transform at the Nyquist limit. At t = 5.0, the transform even turns upward at this limit, indicating aliasing problem. Moreover, for the transform at the other two time instants, one-dimensional filter attenuates moderate to large wavenumbers as compared to twodimensional filter. 3.4. Algorithmic cost estimate for the proposed 2D filters Explicit filtering is a post-processing operation. Estimation of the algorithmic cost of two-dimensional filter in terms of addi-

tional requirement of computing time, is attempted next. For this purpose, we have considered a case shown in Fig. 9b of flow past SHM-1 aerofoil when a 2D fourth order filter is applied after every time step. For the filter, we have used a2 = 0.245 and it is applied on first 20 g-lines only for vorticity. Frame (i) in Fig. 11 shows computational time required in seconds for solving stream function Eq. (17) for 500 non-dimensional time steps of the Navier–Stokes equation solver. The average time required to calculate the stream function equation is 0.022 s. The average time required for every time increment of Navier–Stokes equation is 0.2815 s, as shown in frame (ii). Compared to these values, applying 2D filter needs 0.0129 s. This is almost 22 times lower compared to time required for every time increment of Navier– Stokes equation. In addition to it, one has the explicit control over the frequency of filtering. Reduction in computational cost can be achieved by applying 2D filter infrequently. Thus 2D filters can be used effectively with minimum additional cost of computation. This makes an alternate case for using filters for LES and DES of 2D and 3D flows.

4. Summary and conclusions In the present work, a new class of explicit 2D filters has been proposed that is different from the conventional one-dimensional filters. Filters are used to (i) stabilize numerical methods [4– 6,22,31]; (ii) remove aliasing problem [22]; (iii) band-limit unknowns, as needed in LES [9,10,22] and also (iv) remove spurious upstream propagating high wavenumber components, identified as Gibbs’ phenomenon [15] and q-waves [22]. It has been reported that application of 1D filters in different directions causes distortions of flow structures by smearing and diffusion [21]. These are the additional motives for introducing 2D filters, that fulfills above mentioned roles while doing it expeditiously. In Fig. 1a and b, transfer functions for the proposed 2D second, fourth and sixth order filters are shown along h = 0° and 45° angles with respect to kihi axis. Presented central filters have purely real transfer function and hence do not interfere with the dispersion relation of the basic numerical method. Transfer function plots show that it remains close to one when the 2D filter coefficient, a2 is near its possible highest value of 0.25, except along h = 45° where the transfer function is fixed to zero at the Nyquist limit. This helps to control aliasing problem with least attenuation of low wavenumber components by 2D filters, as shown in Fig. 1c. Filtering by 1D filter in two directions successively for the case of skewed wave propagation can cause numerical errors to contaminate the solution near one of the boundary segment of the computational domain, as shown in Fig. 2. For the same problem, use of 2D filters avoids such build up of numerical instability. Comparison of propagation of vortices in the wake of a cylinder executing rotary oscillation in Fig. 3, for different filtered cases with experimental results [25] helps us identify that fourth and sixth order 2D composite filters interfere least while stabilizing calculations. Quantitative comparison of normalized drag coefficient for different rotary oscillation cases given in [25] shows excellent match with the numerically obtained results using a 2D fourth order filter as shown in Table 1. Figs. 4–6 show that the problem of vortex elongation, smearing in one filtering direction and attenuation associated with 1D filters to be less for 2D filters. Apart from removing directionality via filtering with 2D filters, computational effort also can be brought down significantly. This is due to the fact that 2D filters are used only over a small region of the domain close to the body. This has been studied for aerofoil cases shown in Figs. 7–10. Vorticity contour plots in Fig. 7a and b show that the use of a fourth order 2D composite filter with a2 = 0.249 captures vortical structures more accurately as com-

T.K. Sengupta, Y.G. Bhumkar / Computers & Fluids 39 (2010) 1848–1863

pared to the second and the sixth order 2D filters. This also suggests that there is no need to go for higher than fourth order 2D filter. Application of fourth order 2D filter for flow past a NLF aerofoil (SHM-1), shows that the proposed 2D filters can be applied for high Reynolds number flows. The algorithmic cost in terms of computing time required for proposed 2D filter is explained with the help of Fig. 11. Application of 2D filter needs almost 22 times less time as compared to that required for time advancing the Navier–Stokes equation. With additional choice on frequency of filtering, one can further bring down the computing cost. One would like to develop upwind 2D filters, as in [22] for 1D filters. Also, there is a need to extend this for 3D flow field using similar 3D filters. This aspect will help one to use it for LES and DES. References [1] Bogey C, Bailly C. On the application of explicit spatial filtering to the variables or fluxes of linear equations. J Comput Phys 2007;225:1211–7. [2] Freymuth P. The vortex patterns of dynamic separation, a parametric and comparative study. Prog Aerosp Sci 1985;22:108–61. [3] Fujino M, Yoshizaki Y, Kawamura Y. Natural-laminar-flow airfoil development for a lightweight business jet. J Aircraft 2003;40(4):609–15. [4] Gaitonde DV, Shang JS, Young JL. Practical aspects of higher-order accurate finite volume schemes for wave propagation phenomena. Int J Numer Methods Eng 1999;45:1849–69. [5] Gaitonde DV, Visbal MR. Further development of a Navier–Stokes solution procedure based on higher-order formulas. In: 37th aerospace sciences meeting and exhibition, AIAA 99-0557, Reno, NV; 1999. [6] Gaitonde DV, Visbal MR. Padè-type higher-order boundary filters for the Navier–Stokes equations. AIAA J 2000;38(11):2103–12. [7] Najjar FM, Tafti DK. Study of discrete test filters and numerical approximations for the dynamic subgrid-scale stress model in turbulent channel flow simulations. Phys Fluids 1996;8(4):1076–88. [8] Protas B, Wesfreid JE. Drag force in the open-loop control of the cylinder wake in the laminar regime. Phys Fluids 2002;14:810–26. [9] Rizzetta DP, Visbal MR, Blaisdell GA. A time-implicit high-order compact differencing and filtering scheme for large-eddy simulation. Int J Numer Methods Fluids 2003;42:665–93. [10] Rizzetta DP, Visbal MR, Morgan PE. A high-order compact finite-difference scheme for large-eddy simulation of active flow control. In: 46th aerospace sciences meeting and exhibition, AIAA 2008-526, Reno, NV; 2008.

1863

[11] Sagaut P. Large eddy simulation for incompressible flows. Springer-Verlag; 2002. [12] Sengupta TK, Nair MT. Upwind schemes and large eddy simulation. Int J Numer Methods Fluids 1999;31:879–89. [13] Sengupta TK, Ganeriwal G, De S. Analysis of central and upwind compact schemes. J Comput Phys 2003;192(2):677–94. [14] Sengupta TK. Fundamentals of computational fluid dynamics. Hyderabad, India: Univ. Press; 2004. [15] Sengupta TK, Ganeriwal G, Dipankar A. High accuracy compact schemes and Gibbs’ phenomenon. J Sci Comput 2004;21(3):253–68. [16] Sengupta TK, Vikas V, Johri A. An improved method for calculating flow past flapping and hovering airfoils. Theor Comput Fluid Dyn 2005;19(6):417–40. [17] Sengupta TK, Sircar SK, Dipankar A. High accuracy schemes for DNS and acoustics. J Sci Comput 2006;26(2):151–93. [18] Sengupta TK, Gaurav Kumar. Bluff-body flow control by aerodynamic tripping. Presented at ASME PVT 2006/ICPVT-11 conference on flow-induced vibration, Vancouver, Canada; July 2006. [19] Sengupta TK, Deb K, Talla SB. Control of flow using genetic algorithm for a circular cylinder executing rotary oscillation. Comput Fluids 2007;36:500–78. [20] Sengupta TK, Dipankar A, Sagaut P. Error dynamics: beyond von Neumann analysis. J Comput Phys 2007;226:1211–8. [21] Sengupta TK, Lim TT, Sajjan SV, Ganesh S, Soria J. Accelerated flow past a symmetric aerofoil: experiments and computations. J Fluid Mech 2007;591(1):255–88. [22] Sengupta TK, Bhumkar Y, Lakshmanan V. Design and analysis of a new filter for LES and DES. Comput Struct 2009;87:735–50. [23] Sengupta TK, Lakshmananan V, Vijay VVSN. A new combined stable and dispersion relation preserving compact scheme for non-periodic problems. J Comput Phys 2009;228(8):3048–71. [24] Stoltz S, Adams NA, Kleiser L. An approximate deconvolution model for largeeddy simulation with application to incompressible wall-bounded flows. Phys Fluids 2001;13(4):915–97. [25] Thiria B, Goujon-Durand S, Wesfreid JE. Wake of a cylinder performing rotary oscillations. J Fluid Mech 2006;560:123–47. [26] Van Der Vorst HA. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J Sci Stat Comput 1992;13(2):631–44. [27] Vasilyev OV, Lund TS, Moin P. A general class of commutative filters for LES in complex geometries. J Comput Phys 1998;146:82–104. [28] Vichnevetsky R. Numerical filtering for partial differential equation. Numerical applications memorandum. Rutgers University, NAM 156; 1974. [29] Vichnevetsky R, Bowles JB. Fourier analysis of numerical approximations of hyperbolic equations. Philadelphia, USA: SIAM Stud Appl Math 5, SIAM; 1982. [30] Visbal MR, Gaitonde DV. High-order accurate methods for complex unsteady subsonic flows. AIAA J 1999;37(10):1231–9. [31] Visbal MR, Gaitonde DV. On the use of higher-order finite-difference schemes on curvilinear and deforming meshes. J Comput Phys 2002;181:155–85.