An adaptive Non-Linear Frequency Domain method for viscous flows

Computers & Fluids 75 (2013) 140–154

Contents lists available at SciVerse 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

An adaptive Non-Linear Frequency Domain method for viscous flows A. Mosahebi ⇑, S. Nadarajah McGill University, CFD Lab, Department of Mechanical Engineering, Montreal, Quebec, Canada

a r t i c l e

i n f o

Article history: Received 11 June 2011 Received in revised form 5 September 2012 Accepted 16 December 2012 Available online 5 February 2013 Keywords: Adaptive Nonlinear Frequency Domain Method (adaptive NLFD) Vortex shedding 2D cylinder

a b s t r a c t An adaptive Nonlinear Frequency Domain method (adaptive NLFD) has been developed and validated for the Navier–Stokes equations. Although the compute time for periodic flows is reduced by using the NLFD approach over classical time marching schemes, the adaptive NLFD approach leads to further reductions in the overall computational effort. The adaptive NLFD method is validated for the laminar vortex shedding behind a stationary two-dimensional cylinder. A novel adaptation strategy that allows both mode augmentation and reduction, enables the adaptive approach to gain a factor of two reduction in memory and computational cost at an equivalent solution accuracy to the non-adaptive approach. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Periodic phenomena widely appear in physical fluid flow problems where the flow contains repeating patterns in time. For instance, fluid flow around helicopter rotor blades, wind turbines, and gas turbines are some well-known cases associated with periodic compressible flows. In the incompressible regime, there is current interest in developing Micro Air Vehicles (MAVs) that produce lift by a periodic flapping mechanism. Numerical methods based on Computational Fluid Dynamics (CFD) are widely used for the modeling and simulation of fluid flow problems. In the case of periodic problems, usually one is most interested in the periodic solution when the initial transient behaviors vanish. However, because of the hyperbolic nature of the governing equations, simulation of the unsteady part of the solution is inevitable. Therefore, usual time marching schemes start solving the equations from an initial condition and march in time until the periodic condition is achieved. The time step should be selected in such a way that all of the complicated flow physics involved in the transient part are appropriately captured. For many cases, the time associated with the initial transient portion is much larger than the time period. Therefore, most of the computational resources are spent on resolving the initial transient before reaching the desired periodic state. As an example, the fluid flow through an axial turbine is a physical problem with multiple timescales. Yao et al. [1] performed a numerical simulation on a 1 12 stage turbine using the unsteady Reynolds Averaged Navier–Stokes (RANS) model. It was shown ⇑ Corresponding author. E-mail addresses: [email protected] (A. Mosahebi), siva. [email protected] (S. Nadarajah). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.12.016

that, around 2500 iterations are needed to reach the periodic steady state condition. However, only 80 time steps are used to resolve the fundamental period of the oscillation when the solution converged, which is approximately 3.2% of the total computational effort. The vast industrial demands for periodic flow solvers has resulted in new numerical methods, which are developed for this particular category of fluid flow problems. These methods are based on the same general concept, that the solution is periodic over time. This is unlike the time-accurate approaches which admit transient decays as components of the solution. Instead, these methods directly calculate the solution when it is in a periodic state. Given the assumption of periodicity, trigonometric interpolants such as Fourier series, can be used to accurately and efficiently represent the solution. However, the application of these methods to a non-linear system of equations presents some difficulties. Linearized frequency domain and deterministic stress methods [2] are examples of periodic methods that are widely used in industry. In these methods, the conservation equations are linearized and the solution is constructed from the superposition of the steady state mean flow and unsteady perturbations which are periodic in time. The mean flow components are obtained using steady state assumptions and conventional CFD techniques, then any desired number of temporal frequencies of the solution are obtained independently. For this part, the governing equations and the associated boundary conditions are linearized about the mean flow solution. The cost of these methods is proportional to the cost of the steady state solution and product of the number of temporal modes which are computationally cheaper than the mean solution. Unfortunately, for systems that contain strong nonlinearities such as boundary layers, shock waves, and

141

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

vortex shedding, these linear methods do not provide accurate results. In fact, strong nonlinearities in a system violate the basic assumption of these linearized methods which consider the temporal modes independent of each other [3–5]. In 1998, Ning and He [6] used the harmonic approach for numerical simulation of a non-linear system of equations resulted through the introduction of ‘‘unsteady perturbations’’ or ‘‘deterministic stress’’ components. They include the coupling effects of a linearized system of wave equations on a time-averaged mean flow equation. The authors illustrated that the developed method properly addressed the non-linearities within the unsteady viscous flows through turbomachinery blades, which was not observed by the previously linearized frequency method [7,6,8]. In September 2000, Hall et al. [3] used a harmonic balance technique to obtain an efficient periodic solution of a fully non-linear system of equations, where all modes are non-linearly coupled. A pseudo-spectral approach for representing non-linearities is employed in the temporal domain. Unlike Hall’s approach, where the equations are solved in the time domain, McMullen et al. [9,10] developed the Non-Linear Frequency Domain (NLFD) method which solves the equations in the Frequency domain. In this method, the fluxes are evaluated in the time domain and then transformed to the frequency domain via a Fast Fourier Transform technique (FFT). The state vector is then updated in the frequency domain and transformed back to the time domain where the surface and farfield boundary conditions are updated. The residual is a nonlinear function of the state vector in both time and frequency domains; therefore, regardless of the approach, iterative methods can be employed to drive the residual to negligible values. However, since the time derivative would be in pseudo-time rather than real time, various techniques may be used to accelerate the convergence to a steady state solution. Both Hall and McMullen validated their methods against the Euler and Navier–Stokes equations for a number of unsteady periodic problems and their methods successfully accounted for strong nonlinearities [11–13]. The cost associated with spectral methods like McMullen’s NLFD method is proportional to the cost of the steady state solution multiplied by the number of desired temporal modes. For inviscid flow, McMullen has shown that to accurately model an oscillating airfoil pitching about its quarter chord, a temporal resolution of only 1 mode (or equivalently three time samples per period) is needed using the NLFD method [14] versus the 45 time samples needed with a backward difference formulation of the time derivative [15]. These results demonstrate the potential of the method to provide significant reduction in the computational cost for the analysis and design of more realistic problems such as helicopter rotors, turbomachinery, and other unsteady devices operating in the transonic regime. Nadarajah et al. demonstrated a spectral-type convergence rate for the Euler equations [16], further extended the method to three-dimensional viscous flows [17], developed the adjoint equations for both two- and three-dimensional flows, and performed optimum shape design in unsteady viscous flows [16,17]. Although a significant improvement in execution time can be obtained using the NLFD approach, this method has one primary disadvantage. The need for a large number of modes places a restriction on the hardware resources as large amounts of memory are required. The NLFD method solves all modes simultaneously in comparison with usual time marching methods that only need to keep the solutions in one or two previous time levels. To alleviate the memory cost, one possible approach would be to perform a temporal domain decomposition among a number of nodes or blades of a supercomputer, where message passing can be achieved via MPI or OpenMP. However, this approach may still place a large strain on a system that might already require spatial domain decomposition for very large grids. Therefore, there is a

need to develop the necessary algorithms to reduce the memory requirements of the NLFD scheme and subsequently improve the overall computational effort. In many physical problems, regions with high level of nonlinearities are restricted to the wake of an object in external flows, the boundary layer after separation, and around shock waves. Far from these regions, a small number of modes may be sufficient to accurately represent the flow field. It would be inefficient use of memory to allocate the same number of modes to all of the cells in the domain, and instead, the allocation should be done according to the level of local unsteadiness. The primary concept of the adaptive NLFD approach is to augment the Fourier series in each control volume independent of the others based on the level of the normalized wave amplitude of the highest mode. An adaptive harmonic balance approach was initially proposed by Maple et al. [18–20]. The authors demonstrated their approach on a quasi-one-dimensional flow and achieved a 86% reduction in computational effort in comparison with a non-adaptive harmonic balance method [18]. An extension of the adaptive concept to two-dimensional viscous flows was performed by the present author [21]. In this work, a novel adaptive strategy for the adaptive NLFD method is explored for periodic viscous flows. We propose to employ both mode augmentation and reduction as the solution develops. Unlike the previous adaptive approaches [18,21], in which the number of modes could only be increased, the primary advantage of having both mode augmentation and reduction is the ability to recover equivalent solutions to the non-adaptive approach with the least number of modes in the computational domain. Through the proposed adaptive approach, the final converged modal distribution pattern would be independent of the adaptation intervals, employed temporal solver, and flow initial conditions. The paper presents a thorough analysis of the viability of the adaptive approach for laminar viscous flows. Results are presented and validated for a range of Reynolds numbers for the laminar vortex shedding behind a two-dimensional cylinder. For this test case, since the dominant frequency is unknown a priori, a gradientbased-variable-time period technique of McMullen [14] is applied to find the exact time period as part of the solution. A complete analysis and comparison of the computational effort of the approach as well as the effect of Reynolds number is investigated and presented. 2. Governing equations and their discretization The Cartesian coordinates and velocity components are denoted by x1, x2, and u1, u2. Einstein notation simplifies the presentation of the equations, where summation over i = 1–2 is implied by a repeated index i. The two-dimensional Navier–Stokes equations then take the form,

@w @ðfi  fv i Þ ¼ 0 in D; þ @t @xi

ð1Þ

where the state vector w, inviscid flux vector f, and viscous flux vector fv are described respectively by

w¼

8 9 q > > > > > < qu > = 1

> qu2 > > > > > : ; qE

;

fi ¼

8 > > > <

qui

9 > > > =

qu1 ui þ pdi1 ; > qu2 ui þ pdi2 > > > > > : ; qEui þ pui

8 > > > > <

9 > > > = sij dj1 > and f v i ¼ : s d > > ij j2 > > > > > > : uj sij þ k @T ; @xi 0

ð2Þ In these definitions, q is the density, E is the total energy, and dij is the Kronecker delta function. The pressure is determined by the equation of state

142

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

Fig. 1. Simplified dataflow diagram of the time advancement scheme illustrating the pseudo-spectral approach used in calculating the non-linear spatial operator R.

  1 p ¼ ðc  1Þq E  ðui ui Þ ; 2

ð3Þ

where c is the ratio of specific heats. The viscous stresses may be written as

sij ¼ l

    @ui @uj @uk þ þk dij ; @xj @xi @xk

ð4Þ

2.1. Spatial discretization The fluxes are represented in discrete form for each computational cell using a central second-order discretization. Eq. (6) can then be written for each computational cell as,

V i;j

dwi;j þ RðwÞi;j ¼ 0 in D: dt

ð7Þ

The residual can be represented as, where l and k ¼  23 l are the first and second coefficients of viscosity. These equations can be non-dimensionalized and written in integral form for a control volume as follows:

RðwÞi;j ¼ hiþ1;j  hi1;j þ hi;jþ1  hi;j1 ;

pffiffiffi X cM1 X~ ~ @w X~ ~ fc  ~ fd  fv  s ¼ 0 in D: þ V s @t Re1 cv cv cv

quantity is calculated at the cell faces. The values of the flow variables are stored as conservative variables at the cell centers and can be regarded as cell averages. Accordingly, the convective flux at the cell face is computed by taking the average of the flux contributions from each cell across the cell face.

2

2

In this equation, ~ s is the normal surface vector and fc, fv, and fd are the convective, viscous, and artificial dissipation fluxes, respectively. To eliminate odd–even decoupling of the solution and overshoots before and after shock waves, artificial viscosity is added to the convective terms. The artificial dissipation scheme used in the present study is a blended first- and third-order flux, first introduced by Jameson et al. [22]. A finite volume scheme is derived by applying the above equation directly to control volumes to give a set of ordinary differential equations of the form

2

The derivation of the NLFD method starts with Eq. (6) and assumes that the vector of flow variables w and local residual R can be represented by separate Fourier series: m X

2pk

^ k ei T t ; w

R¼

k¼m

ð6Þ

where V is the cell volume, and the residual R(w) is evaluated by summing the fluxes through the cell faces. The simulations performed in the present study are restricted to rigid mesh and therefore the cell volumes remain constant in time.

2

2

2.2. Temporal discretization

w¼ dðwÞ þ RðwÞ ¼ 0 in D; V dt

ð8Þ

2

where hiþ1;j ¼ fciþ1;j  fdiþ1;j  fv iþ1;j . The  12 notation indicates that the 2

ð5Þ

2

m X

b k ei2pT kt ; R

ð9Þ

k¼m

pffiffiffiffiffiffiffi where i ¼ 1. In the case of real problems, the Fourier coefficients for negative wave numbers are simply the complex conjugates of the coefficients for the positive wave numbers and can be eliminated in the computation. The Fourier representations are then substituted into the semi-discrete form of the governing equations as described in Eq. (6) to yield,

Table 1 Grid study at Re = 100. St

C ps

C pb

C Dt

C Dp

C Df

Grid study 128  64 192  96 256  128 384  192

0.18984 0.17046 0.16704 0.16455

1.08323 1.08549 1.08649 1.08711

0.80008 0.75530 0.75268 0.74338

1.44380 1.40090 1.39807 1.39000

1.09164 1.05620 1.05360 1.04682

0.35217 0.34469 0.34448 0.34318

Experimental results Williamson [28] McMullen [14] Henderson [27] Park et al. [26]

0.16434 0.159 – 0.165

– – – 1.04

– 0.697 0.73901 0.73

– 1.323 1.34500 1.33

– 0.98 1.00476 0.99

– 0.342 0.34524 0.34

Error (%) (I) (II)

1.51 0.13

0.06 –

1.25 0.59

0.58 3.35

0.64 4.19

0.37 0.60

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154 Table 2 NLFD results at Re = 100 for various coefficients. Case name

St

C ps

C pb

C Dt

C Dp

C Df

m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13

0.17585 0.16756 0.16642 0.16680 0.16699 0.16703 0.16704 0.16704 0.16704 0.16704 0.16704 0.16704 0.16704

1.08666 1.08647 1.08645 1.08649 1.08649 1.08649 1.08649 1.08649 1.08649 1.08649 1.08649 1.08649 1.08649

0.75329 0.74355 0.74781 0.75193 0.75259 0.75266 0.75268 0.75268 0.75268 0.75268 0.75268 0.75268 0.75268

1.39895 1.39064 1.39380 1.39739 1.39798 1.39805 1.39807 1.39807 1.39807 1.39807 1.39807 1.39807 1.39807

1.05418 1.04758 1.05016 1.05305 1.05352 1.05358 1.05359 1.05360 1.05360 1.05360 1.05360 1.05360 1.05360

0.34478 0.34306 0.34364 0.34434 0.34446 0.34447 0.34447 0.34448 0.34448 0.34448 0.34448 0.34448 0.34448

143

wavenumbers. The unsteady residual c Rk can then be calculated b by adding R k to the spectral representation of the temporal deriv^ k (Eq. (11)). The iteration is advanced and w ^ k is upative, i 2pT k V w ^ k is then transformed dated (Eq. (12)). Using an inverse FFT, w back to the physical space resulting in an updated state vector, w(t), sampled at evenly distributed intervals over the time period. The wall and far-field boundary conditions are imposed within the time domain. A modified five-stage Runge–Kutta time integration scheme is applied to march the solution to a steady state solution. Local time stepping is employed to accelerate the convergence. A diagram detailing the transformations used by the pseudospectral approach at each stage of the modified multistage Runge–Kutta scheme is provided in Fig. 1. 2.3. Adaptive NLFD scheme

Fig. 2. Error of the Strouhal number versus time step for the non-adaptive NLFD method.

V

" # m m X d X 2pk b k ei2pT kt ¼ 0 in D: ^ k ei T t þ R w dt k¼m k¼m

ð10Þ

As a result,

2pk c b k ¼ 0 in D; ^k þ R Rk ¼ i Vw T

ð11Þ

forms the new residual in the frequency domain for each wavenumber and must be solved iteratively using a pseudo time marching concept,

^k dw V  þc Rk ¼ 0 in D: dt

ð12Þ

The solver attempts to find a solution, w, that drives this system of equations to zero for all wavenumbers. The nonlinearity of the unsteady residual and, hence, bulk of the computing effort stems from the spatial operator, R. The pseudo-spectral approach begins by initializing the state vector, w(t), at all time instances. At the first iteration, w(t) will take on the value of the initial condition and for subsequent iterations, the values are based on the previous iterations. At each of these time instances, the operator R(w(t)) can be computed by summing the convective, viscous, and artificial dissipation fluxes as mentioned in the previous subsection. A FFT is then used to transform the state vector and spatial operator to b k are known for all ^ k and R the frequency domain where w

The primary concept of the adaptive NLFD approach is to augment the Fourier series in each control volume independent of the others, based on the level of the Normalized Wave Amplitude (NWA) of the highest mode. This allows for a large reduction in the total degrees of freedom (DOF) required in the temporal domain to accurately resolve the solution at a desired level of accuracy as well as it allows for very high number of modes to be introduced in regions of the flow where a larger amount of unsteadiness is present. Here unsteadiness denotes higher frequencies in the temporal direction of the flow solution which is synonymous to p-adaptation in the spatial direction, where higher-order polynomials are applied in regions of high gradients in the flow solution. The DOF is defined as the sum of the total number of time steps at each control volume. In the case of the nonadaptive approach, the number of time steps is fixed and equal for all control volumes, while this number varies for the adaptive approach. The procedure to augment the Fourier series is as follows. First, the initial solution is set to the freestream conditions and each control volume is assigned three time steps per period or the fundamental frequency. This represents the minimum number of modes that can be specified. As the solution progresses, one additional mode is added to each control volume if the ratio of the normalized wave amplitude of the highest mode is above a Reference Normalized Wave Amplitude (RNWA). The reference normalized wave amplitude is a user specified threshold to maintain a minimum desired level of normalized wave amplitude in the highest mode. The normalized wave amplitude for each mode can be defined as follows:

ck j jw ; NWAk ¼ Pm c k¼0 j wk j

ð13Þ

where

ck j ¼ jw

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffi ck w ^ k ¼ w ck : ck g2 þ If w ck g2 ¼ w ck w Rf w

ð14Þ

The Fourier coefficients of the new mode are initially set to zero and as the solution develops, the residual of this mode will begin to converge. Since it takes several iterations before the solution converges to three or four orders of magnitude, the frequency augmentation procedure is only enforced every several hundred iterations. Besides the capability of increasing the number of modes, the adaptation procedure allows for mode reduction. As the flow develops, some of the cells may initially require higher number of modes but the need seizes once the dominant features of the flow are established and the solution is within the desired level of accuracy. Having the ability to decrease and increase the number of modes, results in a unique mode distribution pattern in

144

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

Fig. 3. Entropy contours for laminar vortex shedding behind a 2D cylinder (Re = 100) for the non-adaptive NLFD method (a) m1, (b) m3, (c) m5, and (d) m7.

the physical domain. To eliminate the last mode, the normalized wave amplitude of the second last mode is compared to the reference normalized wave amplitude. If lower, then the mode is removed, ensuring that the last mode in each control volume is always less than the reference threshold. This approach has the unintended consequence of producing oscillations in the mode adaptations, where a mode augmentation is quickly followed by a mode reduction at the next adaptation cycle, thus potentially affecting the convergence rate of the solver. To eliminate the oscillations, the number of modes are held fixed after a specific number of oscillations occur at a particular control volume. The following pseudo-code, presents the mode augmentation and reduction approach employed in this work for the base adaptation procedure.

do Cell =1, Max m = Cell number of modes Compute SEm and SEm1 if (NWAm< RNWA) if (NWAm1< RNWA) m = m  1 (Decreasing the number of modes) end if else m = m + 1 (Increasing the number of modes) end if end do

145

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154 Table 3 Adaptive NLFD results at Re = 100 for various coefficients. Case name

Reference normalized wave amplitude

St

C ps

C pb

C Dt

C Dp

C Df

m1–m2 m1–m3 m1–m4 m1–m5 m1–m7 m1–m8 m1–m9 m1–m10 m1–m11 m1–m12 m1–m13 m1–m14

2

1.50  10 6.75  103 2.75  103 1.50  103 5.00  104 2.75  104 1.50  104 9.75  105 6.50  105 4.00  105 2.00  105 6.25  106

0.17585 0.16249 0.16626 0.16751 0.16653 0.16680 0.16699 0.16701 0.16701 0.16703 0.16704 0.16704

1.08666 1.08548 1.08643 1.08650 1.08646 1.08647 1.08649 1.08649 1.08649 1.08649 1.08649 1.08649

0.75329 0.70743 0.72126 0.74541 0.74821 1.08647 0.75126 0.75170 0.75217 0.75240 0.75272 0.75268

1.39895 1.35895 1.37133 1.39119 1.39355 1.39548 1.39684 1.39721 1.39762 1.39782 1.39809 1.39807

1.05418 1.02147 1.03196 1.04798 1.04996 1.05151 1.05260 1.05290 1.05323 1.05339 1.05361 1.05359

0.34478 0.33748 0.33937 0.34321 0.34359 0.34397 0.34424 0.34431 0.34439 0.34443 0.34448 0.34448

m1–m8⁄ m1–m8⁄⁄

6.25  106 1.00  107

0.16704 0.16704

1.08649 1.08649

0.75268 0.75268

1.39807 1.39807

1.05359 1.05359

0.34448 0.34448

To investigate the role of the higher modes on the accuracy of the solution, we propose a second alternate adaptation procedure. The strategy is identical to the base procedure described above with one exception and that is the alternate procedure has a limit on the maximum number of modes. Suitable values for the reference normalized wave amplitude as well as the limit are explored in the results section.

in either cell are below the reference normalized wave amplitude, we believe that the error is within the order of accuracy of the specified threshold. If a sufficiently small value for the reference normalized wave amplitude is employed, then the error incurred by the above flux adjustment is negligible.

2.4. Flux adjustment at temporal-mismatched cell faces

For the validation of the developed solver, the laminar vortex shedding behind a stationary 2D cylinder is investigated. The results are compared with the experimental results of Williamson [23–25], and the numerical results of: Park et al. [26] who used a second-order scheme for the spatial discretiztion; an eight-order scheme by Henderson [27]; and McMullen [14] who employed the NLFD method with second-order accurate fluxes. An O-type grid is used for the problem where 256 cells are used in the streamwise, i-direction, and 128 cells in the normal, j-direction, based on a grid study (Section 3.1). The ratio between the outer or farfield circle diameter to the cylinder diameter is 200. The problem is solved in non-dimensionalized form for M = 0.3 and a range of Reynolds numbers from Re = 60 up to Re = 180. In this range of Reynolds numbers, 2D laminar vortex shedding is observed. All presented cases, both non-adaptive and adaptive, were computed with a CFL of 3 for the five-stage Runge–Kutta solver. A gradient approach was developed by McMullen [14] for the class of problems where the time period of the unsteadiness is unknown a priori. In this method, the time period is modified in each iteration and converges to the exact value as part of the solution. The final formula for modifying the time period is as follows;

Apart from the frequency augmentation and reduction, another important issue is the evaluation of fluxes across faces where adjacent cells do not hold the same number of modes. Since the adaptive approach is a nonlinear technique, the fluxes and boundary conditions are computed and enforced in the time domain. Hence, at each temporal-mismatched cell face, a special treatment of the flux is required. For instance, if the left cell has three time steps or one mode, mL = 1, the solution is available at the t ¼ 0; T3, and 2T time instances, while the right cell has five time steps per period 3 or mR = 2, and therefore, the solutions are available at the t ¼ 0; T5 ; 2T ; 3T , and 4T instances. Apart from the solution at t = 0, 5 5 5 there is a mismatch when computing the flux across the face. To evaluate the right flux of the left cell, the highest mode in the right cell is removed resulting into a new number of modes, mR , where mR ¼ mL , and an inverse FFT is applied to recover the solution for a three time step case, as follows: 

wR ðtÞ

¼

mR X

2pk

^ k ei T t : w

k¼mR

The phases match and the corresponding right fluxes of the left   ^ L; w ^ R . The temporary re-sampled solucell are computed, F L ¼ F w tion for the right cell is discarded once the fluxes are evaluated. A similar strategy is applied for the left flux of the right cell, where an additional mode with the Fourier coefficients set to zero is added to the left cell resulting into a new number of modes, mL , such that mL ¼ mR , followed by an inverse FFT to recover the solution for a five time step case, as 

wL ðtÞ ¼

mL X

3. Results

T nþ1 ¼ T n  DT

@j c Rk j2 ; @T

ð15Þ

where

1 @j c Rk j2 2pV ^k  c Rk j; ¼ 2 kjw 2 @T T

ð16Þ

^ k is the state vector in the where c Rk is the unsteady residual and w frequency domain as defined previously.

2pk

^ k ei T t : w

k¼mL

The left fluxes  of the  right cell are then computed at all time in^ L; w ^ R , and the temporary re-sampled solution of stances, F R ¼ F w the left cell is discarded. This allows for a simple approach to compute the fluxes across cell faces with non-matching temporal resolutions; however, it renders the approach a non-conservative scheme. Since the normalized wave amplitudes of the last mode

3.1. Spatial grid study In order to determine the minimum grid size required to generate a grid independent solution, four grids of successive refinement were employed to simulate the flow field at Re = 100. The results are presented in Table 1 for six interested coefficients; the Strouhal number (St) and the time-averaged stagnation pressure coefficient ðCps Þ, base suction coefficient ðCpb Þ, and drag coefficients (total

146

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

Fig. 4. Entropy contours for laminar vortex shedding behind a 2D cylinder (Re = 100) for the adaptive NLFD method (a) m1–m5 (b) m1–m7, (c) m1–m9, and (d) m1–m11.

ðC Dt Þ, pressure ðC Dp Þ, and viscous or friction ðC Df Þ drag coefficients). To ensure the lowest temporal discretization error, 10 modes are considered in the simulations. In Section 3.2, a thorough temporal resolution study confirms that 10 modes is sufficiently small. The error percentage (I) corresponds to the error between the 256  128 grid and the finest grid, 384  192. This error is less than 1.5% for all of the coefficients. The error percentage (II) corresponds to the error between the finest grid and the selected experimental [28] and numerical [27] values, which are bolded. The obtained Strouhal number is within 0.13% of Williamson’s [28] experimental values, and 3.5% of McMullen’s numerical result [14]. The mean or time-averaged base suction coefficient, Cpb , and the viscous drag coefficient, C Df , are within 0.6% of Henderson’s results [27]. However a larger error is observed for the mean pressure drag coefficient and as a consequence mean total drag. The

difference is primarily due to the compressibility effects since the compressible Navier–Stokes equations are solved. It is observed that all of the coefficients approach the experimental and numerical values as the number of cells are increased. Based on the grid study, the 256  128 grid will be employed for all remaining simulations.

3.2. Temporal accuracy study In this subsection, a thorough analysis of the temporal accuracy of both the non-adaptive and adaptive NLFD methods are demonstrated. The temporal resolution study is performed for a flow over a stationary circular cylinder at Re = 100. The objective of this section is to validate the results and the adaptation strategy.

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

147

Fig. 5. Mode distribution contours for laminar vortex shedding behind a 2D cylinder (Re = 100) for the adaptive NLFD method (a) m1–m5 (b) m1–m7, (c) m1–m9, and (d) m1– m11.

3.2.1. The NLFD approach Table 2, presents the six interested coefficients by using the NLFD method and employing 1 up to 13 modes, corresponding to cases m1 through m13. As expected, when the number of modes are increased, the accuracy of the results improved. For instance, if the convergence criteria is set to an accuracy of up to five decimal places for all of the mentioned coefficients, eight modes are required. This confirms that the 10 modes employed to conduct the spatial grid study as reported in Section 3.1 was adequate to ensure that the error from the temporal discretization did not influence the spatial discretization error. Fig. 2 demonstrates the log of the error of the Strouhal number versus the log of the time step, where case m13 is considered to be the benchmark. The time step for each NLFD case, is evaluated by computing the time step in the temporal domain based on the ratio of the time period and the number of

samples (2m + 1), per period, where m is the number modes employed. A spectral-type convergence is observed, where the error decreases by four orders of magnitude with an approximate order reduction in the time step. An interesting observation in Table 2, is the required number of modes to achieve convergence for the six coefficients. Let us examine the values of pressure coefficients on two different points on the surface of the cylinder. At the stagnation point, four modes (m4) is sufficient for the mean stagnation pressure coefficient, while seven modes (m7) are required for the mean base pressure coefficient. The stagnation pressure is not as much influenced by an accurate resolution of the wake structure or the convection of the vortices downstream, while the mean base pressure and drag coefficients as well as the Strouhal number are strongly dependent on these flow features. Fig. 3 presents the entropy distribution con-

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

10

-3

10

-4

10

-5

10

-6

Case m8 of table 2 Case m1-m8 of table 3 Case m1-m14 of table 3 Case m1-m8* of table 3 ** Case m1-m8 of talbe 3

1

2

3

4

5

6

7

8

9

10

11

12

13

14

Mode Fig. 6. Average normalized wave amplitude versus the mode number for the nonadaptive case m8 and various adaptive NLFD cases; m1–m8, m1–m14, m1–m8⁄, and m1–m8⁄⁄.

Table 4 Accuracy of the Strouhal number for non-adaptive and adaptive NLFD cases. Case name

Strouhal number

Case Case Case Case Case

0.167039281189180 0.166804568900642 0.167037899424102 0.167037525077437 0.167039263598752

m8 of Table 2 m1–m8 of Table 3 m1–m14 of Table 3 m1–m8⁄ of Table 3 m1–m8⁄⁄ of Table 3

Table 5 Average, first, and second mode amplitudes of the total drag coefficient. m8

m1–m10

m1–m14

m1–m8⁄

C Dt Cd Dt1

1.39807

1.39548

1.39807

1.39807

0.00014

0.00018

0.00014

0.00014

Cd Dt2

0.00436

0.00425

0.00436

0.00436

tours with 1, 3, 5, and 7 modes. The figure shows that increasing the number of modes results in a better distribution of flow properties since the nonlinearities in the problem can be modelled more precisely. In Fig. 3(a), the resolution of the entropy in the region of the vorticies are less defined. As the number of modes are increased, the resolution greatly improves with an almost identical visual comparison between the simulations with 5 and 7 modes as demonstrated in Fig. 3(c and d). 3.2.2. The adaptive NLFD approach Next, we present the temporal accuracy study for the adaptive NLFD scheme in Table 3. The reference normalized wave amplitude was selected to produce an additional mode for each case as presented in Table 3. This was to ensure that a complete comparison between the non-adaptive and adaptive NLFD schemes could be conducted. The m1–mX signifies that the solution has a modal distribution containing one up to X modes. To achieve equivalent results for all six coefficients to the non-adaptive approach, case m1– m14 with a maximum of 14 modes are required at a reference normalized wave amplitude of 6.25  106. In this table, two addi-

tional cases, labeled as m1–m8⁄ and m1–m8⁄⁄, indicate an alternate adaptation approach and are explained at the end of this section. In Fig. 4, the contours of entropy distribution obtained from the following cases m1–m5, m1–m7, m1–m9, and m1–m11 are presented. As anticipated, decreasing the reference normalized wave amplitude improves the resolution of the wake structure, where the solutions from m1–m9 and m1–m11 are visually identical. A direct comparison between the two schemes, as demonstrated in Figs. 3(d) and 4(d), clearly illustrate identical entropy distributions, validating the adaptive approach. Fig. 5(a–d) illustrates the modal distribution for the corresponding adaptive cases. The figures demonstrate that as the reference normalized wave amplitude reduces, the region upstream of the cylinder is resolved with at most two modes due to the lower level of unsteadiness at the stagnation point; however, the wake region sees an increase in the higher modes with a concentration in the convection path of the vortices. The two lobes observed in the region of the wake are centred along the path of the shedding of the vortices with a concentration at the vortex cores. This confirms the observation in Table 2, as was seen for the NLFD simulations, that the mean stagnation pressure coefficient compared to the mean base pressure and viscous drag coefficients required a smaller number of modes to converge to an equivalent level of accuracy. As it was mentioned, Table 3 also contains two additional test cases, designated as cases m1–m8⁄, and m1–m8⁄⁄. For case m1– m8⁄, the reference normalized wave amplitude is set to a lower level of 6.25  106, equivalent to that employed for case m1–m14, but the achievable maximum number of modes is frozen at eight. The same is done for case m1–m8⁄⁄ with the difference that the reference normalized wave amplitude is set to 1.0  107. These cases represent the solutions with the alternate adaptive procedure, where a low reference normalized wave amplitude is employed, while freezing the highest attainable number of modes. This approach further reduces the computational cost of the adaptive approach. Here eight modes were chosen, since the non-adaptive case m8 of Table 2 resolved the accuracy of the six coefficients at the same level. As demonstrated in Table 3, cases m1–m8⁄ and m1–m8⁄⁄, reproduced identical results for the six coefficients. A thorough comparison of the various approaches are detailed in the next subsection.

Case m1-m8 of table 3 Case A * Case m1-m8 of table 3 Case B

15000

Number of Cells

Average normalized wave amplitude

148

10000

5000

0 1

2

3

4

5

6

7

8

Mode Fig. 7. Statistical distributions of modes among cells. Case A: post-processed result of case m8 of Table 2 at RNWA of case m1–m8 of Table 3 . Case B: post-processed result of case m8 of Table 2 at RNWA of case m1–m8⁄ of Table 3.

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

149

Fig. 8. The mode distribution contours for cases m1–m8 and m1–m8⁄ of Table 3, and their corresponding post-processed results of m8 of Table 2: (a) case m1–m8, (b) postprocessed result of case m8 with reference normalized wave amplitude of case m1–m8, (c) case m1–m8⁄, and (d) case m8 with reference normalized wave amplitude of case m1–m 8⁄.

3.2.3. Comparison of NLFD and adaptive NLFD approaches We compare the two approaches thoroughly to ensure that the adaptive approach recovers the non-adaptive solution with three different analysis. First, we compare the average normalized wave amplitudes for each mode between cases m8, m1–m8, and m1– m14 in Fig. 6. Average or mean normalized wave amplitude of a mode is obtained through an arithmetic averaging of the normalized wave amplitudes of the desired mode among all of the computational cells that have that mode in their Fourier series representation. The figure also contains the solutions from cases m1–m8⁄, and m1–m8⁄⁄ to demonstrate the advantages of the alternate adaptive procedure. The figure clearly illustrates that as the

reference normalized wave amplitude is reduced, the adaptive approach begins to recover the correct distribution of the average normalized wave amplitude at each mode of the benchmark case m8. The Strouhal numbers listed until 15 decimal places in Table 4, confirms that the accuracy of the Strouhal number improves as the reference normalized wave amplitude is reduced. This level of accuracy is certainly unnecessary; however, for the purpose to undeniably confirm the accuracy of the adaptive approach as well as understand the role of the reference normalized wave amplitude, it is necessary to provide a thorough validation. A comparison of cases m1–m14 and m1–m8⁄ confirms that by employing the

150

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

Table 6 NLFD results at Re = 160 for various coefficients.

Table 8 Comparison of the computed coefficients with previous experimental/numerical results for Re = 160.

Case name

St

C ps

C pb

C Dt

C Dp

C Df

m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14

0.20898 0.19441 0.18921 0.18920 0.19021 0.19050 0.19055 0.19058 0.19059 0.19060 0.19060 0.19061 0.19061 0.19061

1.06399 1.06425 1.06402 1.06398 1.06401 1.06403 1.06403 1.06403 1.06403 1.06403 1.06403 1.06403 1.06403 1.06403

0.93650 0.90941 0.91914 0.93426 0.93701 0.93659 0.93635 0.93637 0.93641 0.93643 0.93643 0.93643 0.93643 0.93643

1.39529 1.38231 1.38595 1.39708 1.39962 1.39946 1.39931 1.39933 1.39937 1.39939 1.39939 1.39939 1.39939 1.39939

1.11389 1.10294 1.10625 1.11568 1.11779 1.11765 1.11753 1.11754 1.11758 1.11759 1.11759 1.11759 1.11759 1.11759

0.28141 0.27937 0.27970 0.28141 0.28183 0.28181 0.28179 0.28179 0.28180 0.28180 0.28180 0.28180 0.28180 0.28180

St

C pb

C Dt

C Dp

C Df

Present study Previous studies

0.19061 0.18641a

0.93643 0.91063b

1.39939 1.33371b

1.11759 1.05256b

0.28180 0.28115b

Difference (%)

2.68

2.74

4.94

6.17

0.35

a b

Experimental results of Williamson [23]. Numerical results of Henderson [27].

-2

RNWA = 10-3 RNWA = 10 RNWA = 10-4 RNWA = 10-5 -6 RNWA = 10 Park et al. (1988) Williamson (1988)

same reference normalized wave amplitude but limiting the maximum attainable number of modes to eight, case m1–m8⁄ is able to recover a comparable distribution of the average normalized wave amplitude at each mode and yield an identical Strouhal number to case m1–m14 up to the desired five decimal places. The figure demonstrates that the presence of modes higher than eight do not contribute towards the accuracy of the Strouhal number and that what is of greater importance is the correct distribution of the average normalized wave amplitudes. Table 5 further validates case m1–m8⁄, by demonstrating that the mean, first, and second modes of the drag coefficient converge to the values obtained by cases m8 and m1–m14, while m1–m8, differs significantly. Decreasing the reference normalized wave amplitude further to 1.00  107, case m1–m8⁄⁄ provides two additional significant digits of accuracy. Second, we compare, as illustrated in Fig. 7, the distribution of the number of cells assigned to each mode for both approaches. In the case of the non-adaptive technique, at each control volume modes that contain a normalized wave amplitude lower than the reference normalized wave amplitude of the corresponding adaptive approach are removed. At a reference normalized wave amplitude of 2.75  104, the modal distribution of the post-processed case m8 (indicated as case A in Fig. 7) compares very well to case m1–m8. The same is true at the reference normalized wave amplitudes of 6.25  106 for the adaptive case m1–m8⁄ and its corresponding post-processed non-adaptive case m8 (indicated as case B in Fig. 7). As the reference normalized wave amplitude is reduced, the peaks gradually move towards higher modes. The jump

Strouhal Number

0.2

0.18

0.16

0.14 60

80

100

120

140

160

180

Re Fig. 9. Strouhal number versus Reynolds number for the adaptive NLFD method with various values of reference normalized wave amplitude.

observed at the last mode in the alternate adaption approach is due to the limit on the maximum number of modes. For case m1–m8⁄, around 4% of the cells reach the highest attainable number of modes. By reproducing comparable modal distributions throughout the computational domain, the adaptive method is further validated. Third, we compare contour plots of the modal distributions. Fig. 8(a) confirms that the modal distribution for case m1–m8 is al-

Table 7 Adaptive NLFD results at Re = 160 for various coefficients. Case name

Reference normalized wave amplitude

St

C ps

C pb

C Dt

C Dp

C Df

m1–m1 m1–m3 m1–m4 m1–m5 m1–m7 m1–m8 m1–m9 m1–m10 m1–m11 m1–m12 m1–m13 m1–m14 m1–m15 m1–m16 m1–m17 m1–m18 m1–m19 m1–m20

2

2.25  10 1.00  102 4.25  103 3.00  103 1.25  103 7.50  104 4.00  104 2.50  104 1.75  104 1.00  104 6.75  105 4.50  105 2.50  105 1.75  105 7.50  105 4.50  106 1.50  106 3.00  106

0.20898 0.17381 0.19016 0.19513 0.18943 0.18917 0.18954 0.19018 0.19068 0.19063 0.19056 0.19057 0.19056 0.19058 0.19058 0.19060 0.19060 0.19061

1.06399 1.06399 1.06440 1.06441 1.06408 1.06402 1.06400 1.06402 1.06403 1.06403 1.06403 1.06403 1.06403 1.06403 1.06403 1.06403 1.06403 1.06403

0.93650 0.84941 0.84308 0.89764 0.91781 0.92245 0.93385 0.93340 0.93671 0.93740 0.93700 0.93684 0.93629 0.93628 0.93635 0.93638 0.93641 0.93643

1.39529 1.33130 1.33108 1.37103 1.38234 1.38872 1.39686 1.39693 1.40000 1.39996 1.39964 1.39958 1.39924 1.39927 1.39932 1.39935 1.39937 1.39939

1.11389 1.05996 1.05947 1.09320 1.10316 1.10861 1.11547 1.11551 1.11772 1.11807 1.11780 1.11775 1.11747 1.11749 1.11753 1.11756 1.11758 1.11759

0.28141 0.27134 0.27161 0.27783 0.27918 0.28012 0.28139 0.28142 0.28200 0.28189 0.28180 0.28183 0.28178 0.28178 0.28179 0.28179 0.28180 0.28180

m1–m12⁄

3.00  106

0.19061

1.06403

0.93643

1.39939

1.11759

0.28180

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

-2

RNWA = 10-3 RNWA = 10-4 RNWA = 10 RNWA = 10-5 RNWA = 10-6 Henderson (1995) Williamson and Roshko (1990)

Base Pressure Coefficient (-Cpb )

0.9

0.8

0.7

0.6 60

80

100

120

140

160

180

Re Fig. 10. Base suction coefficient versus Reynolds number for the adaptive NLFD method with various values of reference normalized wave amplitude.

RNWA=10-2 -3 RNWA=10 RNWA=10-4 -5 RNWA=10 RNWA=10-6 Henderson (1995)

Viscous Drag Coefficient

0.4

0.35

151

of 3.00  106 (case m1–m20 of Table 7), yield identical results. At an equivalent reference normalized wave amplitude with a limit on the maximum number of modes to 12, case m1–m12⁄ of the alternate adaptive procedure, returns an identical set of coefficients. This again confirms that there is no effect of the 13th through the 20th modes of case m1–m20 on the final result and that at a lower reference normalized wave amplitude, case m1– m12⁄, is able to produce identical results. In this case 5.5% of the cells reach the 12th mode which is the maximum attainable number of modes. A comparison of the obtained numerical results with other experimental and numerical results is presented in Table 8. A larger discrepancy is observed when compared to the Re = 100 case due to the presence of three-dimensional effects as reported by Henderson [27]. A comparison of the solutions at Re = 100 and Re = 160, reveals that as the Reynolds number increases, the strength of the unsteadiness in the domain increases and consequently, a lower reference normalized wave amplitude as well as higher number of modes are required. In order to investigate the generality of the solver, the problem is solved for a range of Reynolds numbers from Re = 60 to Re = 180 with Reynolds number increments of 20. In Figs. 9–11 the results of the adaptive method with different values of reference normalized wave amplitude are compared with previous numerical results of Henderson [27], and experimental results of Williamson [23]. In Fig. 9, the Strouhal number as a function of Reynolds number is presented and the adaptive results converge towards Williamson’s proposed correlation (St = 3.3265/ Re + 0.1516 + 0.00016Re [28]). For comparison of the mean base suction and mean drag coefficients, the correlation functions of Henderson [27], f ðReÞ ¼ a0  a1  Rea2  expða3  ReÞ for C pb and f ðReÞ ¼ a0 =Rea1 for C Df , where the constants are listed in [27], are employed. Figs. 10 and 11, show that the time-averaged base pressure and viscous drag coefficients converge towards Henderson’s numerical correlation functions. 3.4. Computational cost and memory analysis

0.3

60

80

100

120

140

160

180

Re Fig. 11. Viscous drag coefficient versus Reynolds number for the adaptive NLFD method with various values of the reference normalized wave amplitude.

most identical to Fig. 8(b), which presents the modal distribution for the post-processed case m8. Similarly the modal distribution for case m1–m8⁄ is comparable to case m8 as illustrated in Fig. 8(c and d). Both adaptive procedures, the base and the alternate, produced equivalent coefficients, and spectral and modal distributions, while the later approach demonstrated the greater importance of the correct distribution of the normalized wave amplitude. 3.3. Reynolds number effect In this subsection, the effect of Reynolds number is investigated for further validation of the adaptive scheme. Tables 6 and 7, present the results at Reynolds number 160. At the same convergence criteria of five decimal place accuracy for all six coefficients, the non-adaptive method with 12 modes (case m12 of Table 6) and the adaptive method with a reference normalized wave amplitude

Based on the analysis conducted in subSection 3.2.3, of the difference between the non-adaptive and adaptive schemes, in this subsection we will provide a comprehensive comparison of the computational cost and memory usage. The convergence of the L2 norm of the residual corresponding to each mode, c Rk , is presented in Fig. 12(a) for case m8 and in Fig. 12(b) for case m1– m8⁄ at Re = 100. In the adaptive method, the adaptation is done at every 500 iterations, where a new mode is added to or removed from a control volume when the normalized wave amplitude of the highest modes is greater or smaller than the reference normalized wave amplitude. This is represented as region (I) in Fig. 12(b), where the jumps in the residuals are due to the addition of new modes, which are initially set to zero but quickly converges since the residual of all the other cells at the same mode level might be orders of magnitude lower. In region (II), neither adaptation nor redistribution of the modes are present. Its important to note that the adaptation and the reconstruction of fluxes across mismatched cell faces does not compromise the stability of the discretization of the flow solver. All cases presented in this work, achieved a machine zero convergence. Fig. 13 presents a quantitative comparison between the approaches. The error estimate is the accuracy of the mean base suction coefficient, ðCpb Þ, where m13 at Re = 100 and m20 at Re = 160 represent the benchmark solutions and the definition of DOF is provided in Section 2.3. Three important observations can be drawn from Fig. 13. First, at equivalent levels of accuracy, for instance three decimal places for the base suction coefficient at Re = 100 (the error should be 5  104), the adaptive approach offers a memory saving factor of 1.84 due to a reduction in the re-

152

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154 0

-2

10

10-4

100

Mean flow component st 1nd mode 2rd mode 3 mode 4th mode 5th mode th 6th mode 7th mode 8 mode

I T

-6

10

TI

-8

10

TI

10-10 TI

-12

T I

10

0

10000

20000

30000

T I

40000

T I

50000

L 2 norm of the average residual for each mode

L2 norm of the average residual for each mode

10

Mean flow component 1st mode 2nd mode rd 3th mode 4th mode 5 mode 6th mode 7th mode th 8 mode

-2

10

I

-4

10

I

T

T

-6

10

10-8

I T

-10

10

I T I

-12

10

T I

I

10-14

0

60000

II

10000

20000

30000

T

40000

Iteration Number

Iteration Number

(a)

(b)

50000

60000

NLFD (Re=160) adaptive NLFD (Re=160) NLFD (Re=100) adaptive NLFD (Re=100)

Error of Cp b

10-3

-4

10

10-5 0

200000

400000

600000

800000

1000000

DOF

L2 norm of the average residual for the mean flow component

Fig. 12. L2 norm convergence of the average residual for each mode versus iteration number for Re = 100: (a) Case m8, (b) Case m1–m8⁄.

10

0

Case m8 of table 2 Case m1-m14 of table 3 Case m1-m8* of table 3

-2

10

10-4 10-6

Acurate results up to 5 decimal places

-8

10

-10

10

-12

10

10-14

0

200000

400000

600000

800000

CPU Time

Fig. 13. Comparison of the error in the base suction coefficient (Cpb) versus degrees of freedom for the non-adaptive and adaptive NLFD methods.

Fig. 14. L2 norm convergence of the average residual for the mean flow component or zeroth mode versus CPU time at Re = 100 for cases m8, m1–m14, and m1–m8⁄.

quired DOF from 360,448 to 196,390. This translates to an effective computational speed-up of approximately 1.65, due to the additional cost of reconstructing the fluxes across mismatched cells. Second, the memory reduction and computational speed-up increases as the desired level of accuracy increases. If an accurate result up to five decimal places is required for the same case (error of ðCpb Þ should be 5  106) as observed in Fig. 13, the memory saving factor increases to 2.39 (DOF is dropped from 557,056 to 233,318) with a corresponding speed-up factor of 2. Here as the reference normalized wave amplitude is reduced, the ratio of the highest mode in the computational domain to the lowest mode increases and hence the computational advantage of the adaptive approach expands. Third, a slightly greater computational speed-up is observed at the higher Reynolds number of 160, due to the increased amount of unsteadiness observed in the flow as evident in Table 7, where a larger number of modes are required to resolve

the solution when compared to the Re = 100 case. From Fig. 13, it is observed that to acquire an accurate result up to five decimal places at Re = 160, the non-adaptive NLFD method requires 819,200 DOFs, while the adaptive method utilizes 305416. This results into a memory saving factor of 2.68 and a speed-up factor of 2.27 for the adaptive NLFD method. This substantiates the effectiveness of the adaptive approach for viscous periodic flows that have a need for high frequency content or large number of modes in localized regions of the flow. The computational edge of the adaptive approach grows when the amount or level of high frequency content as well as the ratio to the lowest mode increases. Fig. 14 shows the residual convergence of the mean flow component for cases m8, m1–m14, and m1–m8⁄ versus CPU time at Re = 100, demonstrating the approximate speed-up factor of two. The figure demonstrates better computational efficiency for the alternate adaptation procedure in comparison with the base adaption procedure.

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

normalized wave amplitude at the beginning of the simulation. Further investigations are required to validate the presented approach for a general class of cases.

Case m8 of table 2 Case m1-m14 of table 3 Case m1-m8* of table 3 Cases m1 to m7 of table 2

800000

4. Conclusion

CPU usage time

m7 600000

m6 m5 m4

400000

m3 m2 m1

200000

0

0

20000

153

40000

60000

Iteration Number Fig. 15. CPU time versus iteration number for cases m8, m1–m14, and m1–m8⁄.

In Fig. 15, the execution (CPU) times versus iteration numbers of the adaptive cases m1–m14 and m1–m8⁄ are compared to the non-adaptive solutions at Re = 100. The figure illustrates that cases m1–m14 and m1–m8⁄ provide the same accurate results as case m8 but at a reduced CPU time. It should be mentioned again that almost the same number of iterations are needed for all of the cases to reach the desired level of convergence or machine zero. Moreover, it is observed from Fig. 15 that the behavior of case m1–m14 is similar to that of the non-adaptive solution with a maximum of five modes, and case m1–m8⁄ converges asymptotically to the non-adaptive solution with a maximum of four modes. From these numbers, one may conclude that the adaptive solution with eight modes has the same memory and computational cost as the non-adaptive solution with four modes, m4, but with the accuracy of the solution with eight modes, m8. In conclusion, the results presented in this paper demonstrate the viability of the adaptive NLFD approach and the alternate adaptation procedure provides for a further reduction in the computational cost. However, the presented approach does require the user to provide apriori a value for the reference normalized wave amplitude as well as the limit on the maximum attainable number of modes. For the cases presented in this work, there is a linear correlation between the chosen reference normalized wave amplitude and the desired level of accuracy. Thus the chosen reference normalized wave amplitude should reflect this relationship and provide some initial guidance to the user. In the case of the limit on the number of achievable modes, the presented cases employed a pre-specified value since the objective was to present a comparison between the non-adaptive and adaptive approaches. However, in general, this limit could be determined as the solution develops. From Fig. 6 one may observe that the average normalized wave amplitude of the modes decreases as the mode number increases. Moreover, the number of cells that have higher mode contents decreases in the same manner (Fig. 7). Therefore, it may be considered that if the highest mode is represented by only a smaller than a predefined fraction of the control volumes, then a new upper limit is chosen, otherwise additional modes are allowed to be added to the control volumes. For instance, it was observed that at Re = 100 around 4% and at Re = 160 around 5.5% of the cells reach the attainable maximum number of modes. We thus recommend an attainable maximum number of modes of 6% for similar cases. In such an approach, governed by the desired level of accuracy of the final solution, the user is only required to provide a reference

This paper presents a stable adaptive nonlinear frequency domain technique for the Navier–Stokes equations. Laminar vortex shedding behind a stationary 2D cylinder is investigated for a range of Reynolds numbers from Re = 60 up to Re = 180. The results show good agreement with experimental data as well as other numerical results. Two adaptive approaches were explored, one specifying the reference normalized wave amplitude, while the other specifying both the reference normalized wave amplitude and a limit on the maximum attainable number of modes. The mode augmentation and reduction strategy allowed the adaptive method to recover at a comparable level the non-adaptive solution at equivalent reference normalized wave amplitudes. The presence of adapted control volumes does not compromise the stability of the solver. Numerical results show that the adaptive approach can provide a factor of two reduction in the memory and computational cost at equivalent convergence levels. Acknowledgements This research was supported by the McGill University Faculty of Engineering MEDA program and Fonds qubcois de la recherche sur la nature et les technologies, Projet de recherche en quipe. References [1] Yao J, Davis R, Alonso J, Jameson A. Unsteady flow investigations in an axial turbine using the massively parallel flow solver TFLO. Tech. rep., Technical Report 01-0529, AIAA 39th aerospace sciences meeting and exhibit, January 2001. [2] Adamczyk JJ. Model equation for simulating flows in multistage turbomachinery. Tech. rep., NASA technical memorandum 86869. (November 1984). [3] Hall KC, Thomas JP, Clark WS. Computation of unsteady nonlinear flows in cascades using a harmonic balance technique. Tech. rep., 9th international symposium on unsteady aerodynamics, aeroacoustics and aeroelasticity of turbomachines, Lyon, France, September 2000. [4] Chen T, Vasanthakumar P, He L. Analysis of unsteady blade row interaction using nonlinear harmonic approach. J Propul Power 2001;17(3):651–8. [5] Dowell E, Bland S, Williams MH. Linear/nonlinear behavior in unsteady transonic aerodynamics. Tech. rep., AIAA 22nd structures, structural dynamics and materials conference, April 1981. [6] Ning W, He L. Computation of unsteady flows around oscillating blades using linear and non-linear harmonic Euler methods. J Turbomach 1998;120(3):508–14. [7] He L. Modeling issues for computation of unsteady turbomachinery flows. Unsteady flows in turbomachines, von Karman Inst. Lecture Series 1996-05, von Karman Inst. for Fluid Dynamics, Brusseles, Belgium, 1996. [8] He L, Niang W. Efficient approach for analysis of unsteady viscous flows in turbomachines. AIAA J 1998;36(11):2005–12. [9] McMullen M, Jameson A, Alonso J. Acceleration of convergence to a periodic steady state in turbomachinery flows. AIAA paper 01-0152, 39th aerospace sciences meeting and exhibit, Reno (NV), January 2001. [10] McMullen M, Jameson A, Alonso J. Application of a non-linear frequency domain solver to the euler and Navier–Stokes equations. AIAA paper 02-0120, 40th aerospace sciences meeting and exhibit, Reno (NV), January 2002. [11] McMullen M, Jameson A. The computational efficiency of non-linear frequency domain methods. J Comput Phys 2006;212:637–61. [12] McMullen M, Jameson A, Alonso J. Demonstration of nonlinear frequency domain methods. AIAA J 2006;44(7):1428–35. [13] Hall K, Thomas J, Clark W. Computation of unsteady nonlinear flows in cascades using a harmonic balance technique. AIAA J 2002;40. [14] McMullen M. The application of non-linear frequency domain methods to the euler and Navier–Stokes equations. Ph.D. Dissertation, Department of Aeronautics and Astronautics, Stanford University, Stanford (CA), March 2003. [15] Jameson A. Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings, AIAA paper 91-1596, AIAA 10th computational fluid dynamics conference, Honolulu (HI), June 1991. [16] Nadarajah S. Convergence studies of the time accurate and non-linear frequency domain methods for optimum shape design. Int J Comput Fluid Dynam 2007;21(5–6):189–207.

154

A. Mosahebi, S. Nadarajah / Computers & Fluids 75 (2013) 140–154

[17] Nadarajah S, Jameson A. Optimum shape design for unsteady threedimensional viscous flows using a non-linear frequency domain method. AIAA J Aircraft 2007;44(5):1513–27. [18] Maple RC, King PI, Orkwis PD, Wolff JM. Adaptive harmonic balance method for nonlinear time-periodic flows. J Comput Phys 2004;193:620–41. [19] Maple RC, King PI, Oxley ME. Adaptive harmonic balance solutions to Euler’s equation. AIAA J 2003;41(9):1705–14. [20] Maple RC. Multigrid for adaptive harmonic balance CFD. Tech. rep., 16th AIAA computational fluid dynamics conference, June 2003. [21] Mohamadi AM, Nadarajah SK. An adaptive non-linear frequency domain method for viscous periodic steady state flows. Tech. rep., 48th AIAA computational fluid dynamics conference, January 2010. [22] Jameson A, Schmidt W, Turkel E. Numerical solutions of the Euler equations by finite volume methods with Runge–Kutta time stepping schemes. AIAA paper 81-1259, January 1981.

[23] Williamson. Vortex dynamics in the cylinder wake. Annu Rev Fluid Mech 1996;28:477–539. [24] Williamson C. Defining a universal and continuous Strouhal–Reynolds number relationship for the laminar vortex shedding of a circular cylinder. Phys Fluids 1988;31:2742.

[25] Williamson C, Brown G. A series in p1ffiffiffi to represent the Strouhal–Reynolds re number relationship of the cylinder wake. J Fluids Struct 1998;12:1073–85. [26] Park J, Kwon K, Choi H. Numerical solutions of flow past a circular cylinder at Reynolds numbers up to 160. KSME Int J 1998;12(6):1200–5. [27] Henderson RD. Details of drag curve near the onset of vortex shedding. Phys Fluids 1995;7:2102–4. [28] Williamson C. Oblique and parallel modes of vortex shedding in the wake of a circular cylinder at low Reynolds numbers. J Fluid Mech 1989;206:579–627.