A Simple Strategy for Capturing the Unstable Steady Solution of an Unsteady Flow

A simple strategy for capturing unstable steady solution is proposed in this paper. By enlarging the real time step, the unstable high frequency modes of the unsteady flow can be effectively suppressed because of the damping action on the real time level. The steady component will not be changed during real time marching. When the unsteady flow solution converges to a steady state, the partial derivatives of real time in unsteady governing equations will disappear automatically. The steady solution is the exacted unstable steady solution. The method is validated by the laminar flow past a circular cylinder and a transonic buffet of NACA0012 airfoil. This approach can acquire high-precision unstable steady solutions quickly and effectively without reducing the spatial discretization accuracy. This work provides the basis for unstable flow modeling and analysis.


Introduction
Computational fluid dynamics (CFD) has been developed for over 30 years and is able to intensively simulate the complex flow and unstable process of temporal evolution.However, in many engineering studies, particularly for multidisciplinary analysis such as fluid-structure interaction and flow control issues, the direct CFD method is time consuming and difficult to implement for qualitative analysis and parameterization design.Thus, in the last 10 years, reduced order methods [1]- [3] based on the unsteady flow field have been developed rapidly and have been widely used in various engineering fields [4]- [8].
Most reduced order methods are executed on the basis of external excitation for the response to a steady state.These methods are still linear models, which are convenient for stability analysis and control law designing.
Modelling and reducing orders for the unsteady flow-field perturbations based on a steady-state solution are the essence of the linear method.Therefore, obtaining the steady-state numerical solution is the fundamental job of this method.Obtaining the steady-state solution is simple for stable flows.By contrast, acquiring the steadystate solution is usually difficult for unstable flows.However, most problems involved in engineering are closely related to unstable flows, such as vortex-induced vibration problems [9] [10] and buffet problems in aeronautics [11]- [15].Analyzing the flow physical mechanisms and flow-field stability control is significant to model unstable flows.
Unstable flows can be modelled in two ways.One approach involves linearizing the control equations on the basis of steady-state solutions and obtaining the linear equations of small disturbance unsteady components for further reduction [14]- [16].The other approach uses closed-loop control to transform the unstable flow to stable flow, creates a random perturbation for the stable flow, and institutes dynamic models by the system identification method [17] [18].The key point of the second approach is that the designed actuating pattern and control law can turn the unstable flow to a stable flow.Designing the control laws of a complex, nonlinear, and unsteady flow is also difficult.The foundation of the first modelling approach is obtaining the steady-state solution of unstable flows, which usually involves two methods.The first method imposes forced symmetry boundary conditions [19].However, this method is only suitable for symmetric problems, such as flows past a circular cylinder.This method is invalid for asymmetric problems, such as dislocation tandem cylinders, airfoil at high angles of attack with dynamic separation, and transonic buffet problems.The second method needs a feedback control process to obtain a steady solution [17] [20].Designing an appropriate actuating pattern and control law to keep the flow stable is also difficult.A few literature works proposed the use of the time mean flow solution instead of the unstable steady solution.However, the topology of both solutions is significantly different from each other.Moreover, the characteristics of disturbance solutions are observed [19] [21]- [23].Other literature works have adopted the subcritical flow before instability to substitute for the unstable steady solution [24] [25].This method is simple to implement but only takes an expectant effect in the vicinity of an instability critical point; away from the critical point, the method will produce significant errors.
This paper proposes a simple strategy for capturing an unstable steady solution.For the implicit dual timestepping scheme, the high frequency unstable modes will be filtered by enlarging the real time step substantially such that only the steady component remains.Thus, we can acquire a high-precision steady solution that satisfies the control equations and boundary conditions.The effectiveness of the method is verified by a typical flow-field solution including the flow past a circular cylinder and the transonic buffet problem of NACA0012 airfoil.Furthermore, we also analyze the mechanism of the impact on the real time step for the unstable flows.

Governing Equations and Finite-Volume Formulation
The 2D Navier-Stokes flow equations can be written in integral form as follows: ( ) ( ) where Ω is the control volume, ∂Ω is the boundary of the control volume, and In the cell-centered finite-volume method, the computational domain is divided into non-overlapping control volumes that completely cover the domain.The interface variables are derived from the average values of the grid cells to calculate the flux of control volumes.The equations of the integral form are translated to linear ordinary differential equations through spatial discretization, and the flow variables are obtained by the timemarching method.
For the second-order finite-volume method, the semi-discrete finite-volume formulation of the flow equations is expressed as follows: where i Q denotes the cell-centered value of the control volume i,

( )
N i is the sum of cell faces and , i m Γ is the interface area.
The numerical flux can be evaluated by the upwind scheme.According to the Godunov-type method, the interface normal flux is calculated by the Riemann flux: In the current paper, we mainly adopt the Roe scheme [26] to evaluate the numerical flux.For the secondorder accuracy method, the middle point of the control volume interface is used to calculate the numerical integral.

Dual Time-Stepping for Implicit Scheme
For the semi-discrete Equation ( 2), we adopt the implicit discretization for flux terms and the three-point backward-difference approximation of the time derivative: ( ) where The symbol t ∆ denotes the real time step, and the superscript 1 n + is the value of the new time level.Equation ( 4) is solved by the implicit dual time-stepping scheme, which is a particular technique often employed for unsteady flows proposed by Jameson [27].The main idea of the dual time-stepping approach is adding the pseudo-time level to the left of the governing equations.Thus, we obtain the following: ( ) where τ denotes the pseudo-time variable.The implicit time-marching scheme can be employed for the solution of the system of Equation ( 6) on the pseudo-time level.When the scheme is fulfilled, i.e., , we obtain a steady state in pseudo time.The real solution of the new time level is also obtained.The significant advantage of this approach is that the physical time step is unrestricted as in explicit methods.

Mechanism Analysis of the Impact on Real Time Step for the Response of the Unsteady Flow
For the unsteady flow, the unsteady effect in numerical method is mainly reflected by the real time derivative in the governing equations.However, the explicit time-marching scheme must adopt a small time step to ensure the numerical stability because of the restriction of the CFL condition.Thus, the calculation efficiency is low.Therefore, the implicit time-marching scheme is generally used in the unsteady flow-field solution.The process for solving numerical governing equations includes decomposing the computational domain with small nonoverlapping control volumes and transforming the continuous partial differential equations into discrete linear algebraic equations to solve.Considering the truncation error caused by the equation approximation and the dissipation introduced by the numerical scheme, the numerical solution cannot completely and exactly describe the real flow field.Hence, we acquired an approximate simulation.With regard to the physical unsteady flow, such as periodic flow, vortex shedding, and high angle of attack separation flow, the strength of its unsteady characteristics is mainly determined by the instability of flow itself.The necessary condition for the numerical method to accurately simulate the unsteady flow is that the dissipation introduced by the method must be less than the flow instability characteristics.A higher numerical method accuracy corresponds to a smaller solution dissipation.
In general, the steady solution of the unstable flow cannot be obtained by solving the steady governing equations.However, given the dissipation action of the numerical scheme, which embodies the damping effect in flow-field solving, solving the steady governing equations directly can still acquire the steady-state solution if the numerical damping is greater than the flow instability for a few weak unsteady flows.However, for a strong unsteady flow, numerical damping can insufficiently suppress the flow instability, particularly the high-preci-sion method with low dissipation.Furthermore, the steady-state solution cannot be obtained naturally.Thus, many studies adopt coercive boundary conditions to achieve this solution.Expanding the scale of computational grids or using low-order accuracy spatial discrete schemes is the most straightforward way to increase the numerical dissipation.Although this manner may also acquire the steady-state solution, such an approach will greatly reduce the numerical solution accuracy and result in incorrect computation.
We can quickly, efficiently, and accurately obtain the unstable steady solution by enlarging the real time step for unsteady flow-field equations.For the unsteady flows, the flow field contains different magnitude of frequency modes.We need to adopt very small real time step to exactly capture the high frequency modes.While, when Enlarged the real time step, the high frequency modes of unsteady flow are filtered out, whereas the low frequency stable mode will not change.Therefore, enlarging the real time step is equivalent to increasing the numerical dissipation of the real time level.With the proceeding of time marching, when all of the unstable frequency modes decay and only the stable frequency modes remain, the flow field tends to a steady state.In Equation (10), when iterations in the pseudo time-level converge, the pseudo-time derivative tends to be zero and the real solution in the new time level is obtained.When all of the unstable modes are filtered out, the real time derivative ( ) also tends to be zero, and the unstable steady solution is obtained.Furthermore, the steady solution satisfies the steady flow control equations and boundary conditions and does not reduce the numerical accuracy.Thus, we obtain the steady solution for the unstable flow field.

Laminar Flow past a Circular Cylinder
Laminar flow past a stationary circular cylinder is one of the most classical subjects of fluid mechanics and is the fundamental model of investigating the complex flow.Modeling and analyzing the flow past a cylinder can also deepen the understanding of flow physical mechanisms.This section mainly analyzes two states of the Reynolds number including Re = 60 and Re = 100.
In this case, the laminar flow past a stationary cylinder has a free stream of 0.1 and the Reynolds number is 60 and 100.The computational domain is divided by 40324 control volumes and 200 points on the cylinder surface.Furthermore, a hybrid mesh is adopted wherein the first height of the boundary layer is 5 × 10 −3 , the layer sum is 11, and the total height is 0.1.The meshes of the computational domain and the grids near the cylinder are shown in   From the above results, when the real time step is small, the unsteady flow field can be accurately simulated by the numerical method.The shedding vortex behind the cylinder is significantly captured and the lift coefficient results are calculated accurately.When increasing the real time step to the vortex shedding period, the unstable flow field will transform to the steady state and do not change with time.The lift coefficient curve also decreases rapidly to nearly zero.The vorticity distribution of both the time-mean flow and unstable steady flow are all very close to the referenced article (Figure 4(b)).Therefore, the steady flow field obtained by this method does not reduce the accuracy and has high precision.
We perform a calculation by using different sizes of real time step and depict the change of lift coefficient amplitude in Figure 5.In this figure, the y-axis denotes the amplitude of the lift coefficient and the x-axis denotes the multiple of the shedding vortex period T. For the weak unstable flow, Re = 60, the flow field will tend to show a steady state when the real time step increases to 4%T.For the strong unstable flow, i.e., Re = 100, the flow field will tend to show a steady state when the real time step increases to T/2.Thus, the critical real time step of acquiring the steady flow field is related to the instability strength of the flow field.A stronger instability characteristic corresponds to a larger critical real step size.
After obtaining the unstable steady solution, we can recalculate the real unsteady flow field on the basis of the solution as the initial flow field.The lift coefficient response of the cylinder with Re = 100 and dtreal = 0.01T is   shows that a slow and detailed process can be caught from the stable state to the unsteady flow field.By using the unstable steady solution, we can make a better analysis of period flow.

Transonic Buffet Flow of NACA0012 Airfoil
Transonic buffet is a phenomenon of forced vibration because of the interaction between shock waves and boundary layer.This phenomenon is a technical challenge that is often encountered in aircraft design.Modeling and analyzing the buffet flow are particularly important to investigate the mechanism and mitigate or control the buffet vibration.Obtaining the steady-state solution of the unstable flow is the fundamental work of this method.This section mainly demonstrates the transonic buffet flow as the case of NACA0012 airfoil.
In this case, the flow past a NACA0012 airfoil has an angle of attack of α = 5.5˚ and the free stream M ∞ = 0.7 and Reynolds number is 3.0 × 10 6 .The computational domain is divided by 8568 control volumes and 180 points on the airfoil surface.A hybrid mesh is adopted wherein the first height of the boundary layer is 1 × 10 −5 and the total height is 0.02.The meshes of the computational domain and grids near the cylinder are shown in Figure 7.We adopt the SA turbulence model, and other computational method settings are similar to those in Section 3.1.
Given that the reduced frequency of this state is k = 0.195, the corresponding flow period is T = 23, which is uniformed by the non-dimensional parameter of the chord length of the airfoil and the free-stream sound speed.From the above results, when the real time step is equal to 0.01T, the flow field quickly develops to the unstable state and a buffet phenomenon significantly occurs.When the real time step increases to dtreal = T, the flow field no longer presents an unsteady state.The vortex is stationary on the upper airfoil, and the lift coefficient converges and does not change with time.By comparing pressure contours and pressure coefficients, significant     the multiple of the shedding vortex period T, where Cl_min denotes the minimum amplitude of the lift coefficient and Cl_max denotes the maximum amplitude of the lift coefficient.The critical real time step is dtreal = 5%T.When the real time step is larger than the critical step, the flow field will tend to remain stable and have the same size of the lift coefficient.

Conclusion
This paper proposes a fast, effective, and accurate way of capturing the unstable steady solution.By enlarging the real time step size, the unstable frequency modes of the unsteady flow field are filtered out because of the damping in the real time level.When the numerical iteration converges, the time derivatives of unsteady term dismiss automatically and the stable solution suited for steady governing equations are obtained.The method is validated by two unsteady flow cases including the laminar flow past a circular cylinder and a transonic buffet of NACA0012 airfoil.The critical real time step of acquiring the steady flow field is related to the instability strength of the flow field.A stronger instability characteristic corresponds to the need for a larger a critical real step size.The unstable steady solution can be obtained by this method.Thus, we can recalculate the unsteady flow field.Moreover, the detailed process from the linear state to the unsteady nonlinear saturated state can be captured, thus providing the foundation for the modeling and analysis of unstable flows.
the control volume boundary.The vector Q are conservative variables, inviscid fluxes and viscous fluxes respectively.

Figure 2 .
Figure 2. Results of the circular cylinder for the real time step dtreal = T. (a) Vorticity of flow field behind the cylinder; (b) Lift coefficient response of the cylinder.

Figure 3 .
Figure 3. Pressure distribution and streamlines of the time mean flow with Re = 100.

Figure 4 .
Figure 4. Results of the circular cylinder for the real time step dtreal = T. (a) Pressure distribution and streamlines of the flow field; (b) Vorticity comparisons of the unstable steady flow and time-averaged flow: the upper results are obtained from the referenced paper [23], and the lower results are obtained by the present method; (c) Lift coefficients response of the cylinder.

Figure 5 .
Figure 5. Change of lift coefficient amplitude with different sizes of real time step.
First, we set the real time step as dtreal = 0.01, T = 0.23, i.e., 100 points are computed in a single period.The results are shown in Figure 8.The pressure contours and streamlines of four different times in a period are shown in Figure 8(a), and the lift coefficient response of the airfoil is shown in Figure 8(b).The pressure distribution and streamlines of the time mean flow solution at this state are shown in Figure 9.However, when we enlarge the real time step to dtreal = T, a stable state flow field is obtained and the lift coefficient converges.The results are shown in Figure 10.The pressure contours and streamlines are shown in Figure 10(a), and the lift coefficient response of the airfoil is shown in Figure 10(b).The comparison of pressure coefficients with time-mean flow and unstable steady flow is shown in Figure 10(c).

Figure 6 .
Figure 6.Lift coefficient response of the cylinder based on the steady solution for Re = 100.

Figure 8 .
Figure 8. Results of the transonic buffet of NACA0012 airfoil for the real time dtreal = 0.01T.(a) Pressure contours and streamlines of four different times in a period; (b) Lift coefficient response of the airfoil.

Figure 9 .Figure 10 .
Figure 9. Pressure distribution and streamlines of the time mean flow for NACA0012.

Figure 11 .
Figure 11.Minimum and maximum amplitudes of the lift coefficient corresponding to different sizes of the real time step.