TIME-PERIODIC INCOMPRESSIBLE FLOW SIMULATION USING NON-LINEAR FREQUENCY DOMAIN APPROACH

Original scientific paper The capability of nonlinear reduced frequency domain (NLFD) approach in predicting the flow field of time-periodic incompressible flow is investigated. The efficient NLFD method uses Fourier series representation in time considering the assumption of solution periodicity and the resulting nonlinear equations are solved using pseudo-spectral approach. In this manner, the desired periodic solution can directly be obtained without the need to solve the initial transient part. Thus, without loss of generality a considerable computational cost reduction can be achieved in comparison with the common timeaccurate methods. In the present algorithm (INLFD) the NLFD approach is used to solve the flow field in an incompressible formulation and the validation is performed using some periodic test cases with analytical solutions. The results show that capturing only a limited number of temporal modes can provide an accurate estimation of the flow field. Finally, the computational costs of the INLFD and a time-accurate method are also compared which demonstrate the efficiency of the INLFD method in representing nonlinear flow field physics by savings in computational time up to 50 percent.


Introduction
Despite improvements in steady flow simulations, efficient modelling of unsteady flow has been yet a challenge in computational fluid dynamics.In numerical simulation of unsteady flow, the solution time history must be resolved accurately.Thus, the added dimension (time) increases the computational cost effectively.In time-accurate solvers the equations are discretized in space and the time term is left.This would construct a coupled ordinary differential equations system.In order to resolve the solution time history, the method of advancing in time (method of lines) must be used.
By decaying the initial transient state, some unsteady flows reach a periodic steady state solution.In many of these flows, the periodic solution is of more interest than the initial transient one.Since the time which is needed to reach the periodic solution is usually much larger than the steady state period, most of the computational time is spent on resolving the initial transient solution [1].
In order to reduce the long computational time of solving the periodic unsteady flow field, periodic methods can be used.In this approach the periodic solution is calculated directly while the transient decay is not a component of the solution.Considering the solution periodicity, these methods use trigonometric interpolants such as Fourier series instead of the periodic solution which results in a spectral temporal accuracy of the solution versus the number of harmonics used [2][3][4].Periodic methods can be subdivided into three main categories: linearized, semi-nonlinear and non-linear frequency domain methods.
The least expensive approach in calculation of these types of flows is to linearize the flow field about the mean flow steady solution with the assumption that the unsteady (perturbed) component magnitude is much smaller than that of the steady one.In linearized frequency domain methods the mean flow components of the solution are solved first.Any temporal frequency of the solution can then be solved from mean flow components [5].This method has been also extended to potential, Euler and full Navier-Stokes equations for any regime of flow (subsonic, transonic and supersonic) [6][7][8][9][10].
However, the linearized frequency domain methods cannot adequately account for the strong nonlinearities which are typically encountered in flows.To overcome the error associated with the linearized form, seminonlinear and nonlinear forms of frequency domain approach have been developed.In [11][12][13][14][15] semi-nonlinear frequency domain analysis of periodic unsteady flow has been extended.Adamczyk [11] formed deterministic stress terms, using several different linearization and averaging operators and then modelled these terms.These terms effort to quantify the unsteady field effects on the time averaged solution.Ning et al. [12], He et al. [13] and Chen et al. [14] also used averaging operators to form the deterministic stress terms and calculated these terms with a modified version of a linearized frequency domain solver.In deterministic-stress methods, no direct coupling between the harmonics is evident mathematically.So, this method will poorly represent a physical phenomenon wherein significant energy is transferred between different harmonics.
An efficient periodic solution of fully non-linear system of equations was first proposed by Hall et al. [16] using the harmonic balance technique on 2-D turbomachinery cascades.McMullen et al. [17][18][19] also proposed a nonlinear frequency domain method (NLFD) which represents a form of residual in the frequency domain.He focused on 2-D turbomachinery flows and demonstrated the efficiency of this technique to represent complex non-linear flow solutions using a minimum number of harmonic modes.He [4] also showed that flow simulation in frequency domain is an order of magnitude faster than in time domain for the same level of temporal accuracy.
In spectral methods the approximations are defined in terms of a truncated series expansion, such as some quantity (error or residual) which is forced to be zero in an approximate sense [20][21][22].In NLFD method, considering the assumption of periodicity in time, the state vector and fluxes can be represented by truncated Fourier series expansions.The fluxes are evaluated through spatial discretization methods and then transformed to frequency domain by a fast Fourier transform (FFT) technique.The state vector is updated in the frequency domain and transformed back to the time domain.
Murman [23] used the procedure followed by Hall et al. [15] and extended the application to 3-D Euler equations to calculate dynamic stability derivatives of some configurations in an inviscid flow.Also, the results show that in some cases for viscous flows, the solution with just a single mode would generally be irrelevant and higher (but a limited) number of modes should be considered.Mavriplis et al. [24] has also shown that only capturing the dominant modes can provide accurate estimations of fluid properties.
Kharati et al. [25,26] have also used perturbation technique in nonlinear frequency domain method to estimate the solution at high harmonics in transonic flows.The density and velocity field at high harmonics are perturbed about those of low harmonics to form semilinear governing equations.The proposed method has reduced the computational time in comparison with common nonlinear frequency domain approaches.
All previous researches have used non-linear reduced frequency approach in compressible form which solves energy equation alongside with momentum and continuity (or pressure) equations.However, by using incompressible form the energy equation does not need to be solved which results in numerical cost reduction.The innovation of this research is to develop a CFD code which simulates unsteady periodic incompressible flows using non-linear reduced frequency approach (INLFD code).

Methodology 2.1 Governing equations
The primitive variable formulation of the incompressible Navier-Stokes equations in the absence of body forces are given by Eq. (1).
where u, p, ρ and ν are velocity field, pressure, density and kinematic viscosity, respectively.
The conservative two-dimensional form of Eqs.(1a) and (1b) can be written in Cartesian coordinates as where u and v are respectively x and y-components of the velocity vector u.Integrating Eqs. ( 2b) and (2c), over a control volume V and using finite volume approach in which the continuous surface integrals are represented by a discrete summation of fluxes across a finite number of control surfaces, the semi-discrete Eqs.(3.a) and (3.b) will be obtained.
where S x and S y are projected surfaces of control volume in  and  directions, respectively.

Let W=[u,v]
T denote the answer vector.Also the spatial operator R is introduced as a function of space and time including both the convective and dissipative fluxes (Eq.( 4)).
So the momentum equations are simplified in a semidiscrete form:

Transformation of equations into frequency domain
Considering the assumption of periodicity of the solution vector W and spatial operator R in time, both can be represented by truncated Fourier series expansions: where i 2 =−1.
Since the trial (or basis) functions exp(ikωt) are known, determination of the expansion coefficients, k Ŵ and k R , yields acquiring W and R. Therefore, each of the Eqs.(6a) and (6b) would contain N complex coefficients which have to be determined.Since the functions W and R are assumed to be real, two Fourier coefficients with an opposite value of k, are complex conjugate, i.e., Ŵ −k =Ŵ k and The coefficients Ŵ 0 and 0 R being obviously real, the expansions (6) would contain N real unknowns.Eq. ( 7) is also obtained by substituting these discrete Fourier series into Eq.( 5).
Moving time derivative inside the series summation and using Fourier terms orthogonality property, Eq. ( 8) is obtained for each wave number k.In the absence of time derivative term, this equation can be solved just like a steady equation.= 0 Since R is a nonlinear function of W(t), due to convective terms, each coefficient k R depends on all answer transform coefficients, Ŵ k .Thus, Eq. (8) represents a nonlinear set of equations which must be solved iteratively.
Here is a problem of expressing each k R in terms of all coefficients, Ŵ k Hall proposed two approaches to calculate these coefficients [16].In the first approach k R is directly calculated from all Ŵ k using a complex series of convolution sums, which has a massive complexity and computational cost.The second procedure is to use a pseudo-spectral method which relies on the computational efficiency of FFT.This approach is described in the following section.
The other problem to calculate R in incompressible formulation arises from the unknown pressure field, p.There is no explicit equation to obtain pressure, however, the pressure field can be calculated indirectly.Substituting the correct pressure field into momentum equations, the resulting velocity field would satisfy the continuity equation.So, an additional equation which is referred to as pressure Poisson equation (PPE) must be solved.This equation is derived from taking the momentum Eq. (1b) divergence along with continuity Eq. (1a): ( ) So the correct velocity and pressure field can be obtained by solving momentum and pressure Poisson equations iteratively.
The resultant PPE formulation in conjugate with the momentum equation is equivalent to a set of continuity and momentum equations if and only if proper boundary conditions for pressure are specified such that the velocity divergence becomes zero ( 0 = ∇u ) at the boundaries (theorem 1 in [27]).Projecting the momentum equation on the normal vector n at the boundaries, the Neumann boundary condition would be obtained [27][28]:

Solution procedure
The goal is to achieve a solution, W, for which the system of Eqs.(8) converges to zero for all wave numbers.However, at any iteration during the solution process, the unsteady residual, k Î which is defined in Eq.
In order to calculate k R in terms of Ŵ k (for all k), the pseudo-spectral method is used.The first assumption in this method is that Ŵ k is known for all wave numbers.Using inverse FFT, the coefficients Ŵ k transform back to physical space resulting in a state vector, W(t) which is sampled at evenly distributed intervals over the time period.The spatial operator R at each time instance can then be computed from the known W(t) and the pressure p(t), which is obtained by solving the PPE.FFT is also used to transform the spatial operator into the frequency domain.So k R would be obtained for all wave numbers.
In order to compute the unsteady residual k Î at each wave number, k R is added to the spectral representation of the temporal derivative ikVωŴ k .A diagram of the transformations used by the pseudo spectral approach for one iteration is shown in Fig. 1.
Figure 1 The pseudo-spectral approach diagram for incompressible flows The spatial operators computational cost in pseudospectral approach is the product of one steady-state spatial operator cost and the number of time instances used to represent the solution, N. Also, the cost of FFT is proportional to Nln(N) [16].So, for more practical values of N (under 10) the pseudo-spectral approach cost is mainly determined by the cost associated with the calculation of spatial operators.One of the advantages of the pseudo-spectral approach is its flexibility in adopting different forms of non-linear operators.In this research the finite volume formulation is used, however the application of the pseudo-spectral approach is also well suited for finite difference or other types of spatial operators.In addition, different turbulence models which would be difficult to be explicitly expressed in the frequency domain are more easily handled in the pseudo-spectral approach.

Pseudo-time derivatives
Instead of solving Eq. (11) directly a pseudo-time derivative can be added to the unsteady residual.The application of the pseudo-time derivative is consistent with established convergence acceleration techniques used to solve steady-state problems [16].Adding this term, a time-stepping approach can be employed to numerically integrate the resulting Eq. (12).= 0 In Eq. ( 12) the pseudo-time derivative acts as a gradient to drive the absolute value of the unsteady residual to zero for each wave number.

Staggered grid
In order to overcome the problem of oscillating pressure field, a staggered grid arrangement is used in which the scalar variables (pressure data) are located at the cell centres of the control volume, whereas the velocity and momentum variables (u and v velocities) are placed on the cell faces.This is different from a collocated grid arrangement, where all variables are stored in the same position and a finite volume-based discretization employs special interpolation schemes to determine the flux across the cell edges.A staggered storage is mainly used on structured grids.Using a staggered grid is a simple way to avoid odd-even decoupling between the pressure and velocity.Odd-even decoupling is a discretization error which can occur on collocated grids and can then lead to a decoupled (checkerboard) pressure field and oscillations in solutions [29][30].

Test cases description, results and discussion
In order to validate the INLFD code, in this section the incompressible oscillating flow is simulated for some test cases.The emphasis of this research is to investigate the capability of the developed INLFD code in predicting the details of unsteady periodic incompressible flow field and its strength in cost reduction.Although these cases are simple, they have the advantage of involving exact analytical solutions.

Case 1: Stokes second problem
In this case (Fig. 2) an infinite flat plate at the bottom of an infinitely deep sea of fluid with a linear harmonic motion parallel to itself is considered (Stokes' second problem).Since the objective is only to achieve the periodic steady solution, initial conditions are not needed to be satisfied.The velocity distribution is governed by the diffusion Eq. ( 13).
The analytical solution of the Stokes' second problem is given by Eq. ( 14) [31].
The "cc" and "ss" notations are defined as follows: Figure 5 Geometry for unsteady flow between two oscillating plates Figure 6 Velocity profile of the periodic flow between two oscillating plates using different harmonics Fig. 6 shows the exact solution at t = T and the velocity profiles for 3, 5 and 7 harmonics.According to the analytical solution, it is expected that similar to case 1, the numerical solution with 3 harmonics is adequate to provide a good approximation of the flow field.Fig. 7 shows both the exact and INLFD solutions in different instances for the flow between two oscillating plates.Here, two parallel static infinite plates at y = ±h, with an incompressible viscous fluid between them are considered (Identical to test case 2 but with static plates).
The continuity equation is automatically satisfied and the momentum equations are reduced to: where u and p are only functions of x and y, respectively.To produce a pulsating pressure gradient, an oscillating pressure can be considered: The analytical solution is then given by Eq. ( 20) [33].Like previous cases the numerical solution using only 3 harmonics can sufficiently provide a good approximation of the flow field (Fig. 8).Fig. 9 illustrates the INLFD and exact solutions in different instances for the flow produced by a pulsating pressure gradient.As it is seen the INLFD solution is exactly in agreement with the analytical one.

Case 4: Flow between two oscillating co-axial cylinders
Two infinite length concentric circular cylinders of radius r 1 and r 2 are considered (Fig. 10) and the gap between them is filled with an incompressible Newtonian fluid with kinematic viscosity, ν.The outer cup oscillates with velocity U•cos(ωt) while the inner one is kept stationary.Using polar coordinates, the continuity and momentum equations are: Since the flow is axi-symmetric, partial derivatives with respect to θ are zero and the continuity equation results in v r = 0.So the momentum equations are reduced to: Applying the separation of variables method, Eq. (22b) breaks into two independent ordinary differential equations: a t-dependent differential equation and a rdependent one which forms a Bessel differential equation.
Considering the boundary conditions, the analytical solution is: where ( )  Figs. 13 and 14 compare the velocity and pressure profiles for different harmonics at t = T which shows using higher number of harmonics, the numerical solution becomes closer to the exact solution.It is also seen that using 7 harmonics, the exact solution is achieved with a good approximation.Thus using higher number of harmonics is unnecessary.

Conclusion
In this research a CFD code on the basis of non-linear reduced frequency domain approach was developed for simulating periodic incompressible flow.The capability of the developed INLFD code in predicting the details of the unsteady periodic incompressible flow field was then investigated.For this purpose the INLFD code validation was performed for some test cases with exact analytical solutions.For all these test cases, there was a good agreement between the INLFD code results using a limited number of time varying modes and the analytical solutions, which confirm INLFD capability to accurately resolve the flow field.The results show that capturing only the dominant harmonics can provide accurate estimations of the velocity and the pressure field.Also, it is demonstrated that the solution can be approximated with just 3 harmonics for simple cases whereas for more complex cases further harmonics (e.g. 7 harmonics in case 4) are to be used to achieve the desired level of accuracy.As an expected result, increasing the harmonics number brings about more accuracy of the INLFD solution.Afterwards the computational time reduction is also examined for one of the test cases, by comparing the INLFD cost with that of the general time dependent methods.The results show that the INLFD is approximately 2-6 times faster (based on the harmonics number used) than the time-accurate method which proves the claim of the effective computational cost reduction.

Figure 2 Figure 3 Figure 4
Figure 2 Coordinate system of the Stokes' second problem

Figure 7
Figure 7 Velocity profile of the periodic flow between two oscillating plates for different instances 3.3 Case 3: Channel flow with a pulsatile pressure gradient

Figure 8
Figure 8 Velocity profile of the pulsatile pressure gradient channel flow using different harmonics

Figure 9
Figure 9 Velocity profile of the pulsatile pressure gradient channel flow for different instances

Figure 10
Figure 10 Geometry of the flow between two oscillating co-axial cylinders 24b) In Figs.11 and 12 the INLFD and exact solutions are plotted in different instances for the flow between two concentric cylinders.As it is seen the INLFD solution completely follows the exact solution.

Figure 11 Figure 12
Figure 11 Velocity profile of the flow between two oscillating co-axial cylinders for different instances

Figure 13 Figure 14 Figure 15 Figure 16 Figure 17
Figure13 Velocity profile of the flow between two oscillating co-axial cylinders using different harmonics