Pseudo-characteristic upwind scheme for incompressible flows

For enhancing the performance of the characteristic-based method, the multidimensional characteristics of incompressible flows are applied for developing a novel efficient bicharacteristic-based method. The multidimensionality of information propagation is embedded in the convective flux computation by the directional derivatives along the bicharacteristics in conjunction with the cross derivative. The well-known benchmark problems (i.e. the cross flow around a circular cylinder, the lid-driven cavity flow, and the flow around a sphere) are solved to validate the proposed scheme in the presence of strong multidimensional interaction of separated vortical flow. The results represent significant consensus with the benchmark solutions. Furthermore, the proposed scheme results in significantly faster convergence and higher resolution in comparison with the characteristic-based method.


Introduction
The artificial compressibility method (ACM) 1 introduces pseudo-acoustic waves with finite propagation speeds into incompressible flow simulations, which results in a time lag between flow disturbances and their influences on the flow field. 2 The resultant interrupted interaction of pressure field and viscous boundary layer evince inaccuracy and delayed convergence, particularly in separated and viscous internal flows. 2 Therefore, enhancing the flux computation and convergence acceleration are of major interest in pseudo compressible flow simulations.
The main objective of this article is to improve the efficiency of the characteristic-based method (CBM) developed by Drikakis et al. 3 The CBM implements the grid-aligned compatibility equations for upwinding the convective fluxes along the pseudo characteristics. This method was further extended to three-dimensional flows, 4 heat transfer simulations, 5 multigrid techniques, 6 unstructured grids, 5,7 variable density and multispecies flows 8,9 and turbulent flows. 10,11 Despite many successful implementations of CBM, slow convergence was observed in some studies. [12][13][14] Neofytou 15 reported that the source of this inefficiency lies in the mathematical bias of the scheme and proposed a revision for the CBM. The comparative studies of the original and revised CBM revealed that both the methods demonstrate the same performances in terms of accuracy and convergence rate. 16,17 Razavi et al. 13 suggested that the slow convergence of the CBM originates from the assumption of local one-dimensionality of wave propagation parallel to the grid normals, in which the effects of disturbances traveling along other directions are excluded. Their hypothesis was based on the similar problem observed previously in the gridaligned compressible upwind schemes. 18 To incorporate the multidimensional physics of propagation of information, multidimensional characteristic-based method (MCBM) was developed for Cartesian, 13 non-Cartesian, 14 and unstructured grids. 19 Results showed that in addition to the faster convergence, the MCBM takes advantage of an elevated Courant-Friedrichs-Lewy (CFL) number and higher resolution as compared to the CBM. In the MCBM, the wave propagation directions are set to the grid normal and parallel directions, along which four pseudo-acoustic projected Riemann invariants are applied for the estimation of the interface values.
In this article, a novel multidimensional bicharacteristic-based method (BCBM) is developed which differs inherently from the MCBM in both the mathematical bias and the numerical model. From mathematical point of view, for derivation of the multidimensional characteristics of incompressible flow, we followed the same procedure applied for unsteady compressible flows, 20 which results in compatibility relations with two directional derivatives along two independent directions. Since coefficient matrices along coordinate directions cannot be diagonalized simultaneously, 21 the influence of information propagating normal and parallel to the grid surface is taken into account by the directional discretization along the bicharacteristics and cross derivatives (i.e. directional derivatives normal to the wave front), respectively. Furthermore, both the pseudo-acoustic and the shear waves are incorporated to attain a comprehensive wave model that completely covers the physical behavior of the pseudo wave propagation. The proposed method is applied in a cell-centered finite-volume method for the convective flux treatment.

Governing equations
The governing equations of two-dimensional incompressible flow, modified by the artificial compressibility, 1 in the finite-volume form can be delineated as equation (1) dU dt where A ij and k are the cell area and the cell interface number, respectively. Here, the multidimensional bicharacteristic-based upwind scheme is developed for the computation of the convective fluxes.

Pseudo-characteristics of incompressible flows
Taking the linear combination of the reduced governing equations (i.e. without the viscous terms) yields 3,20 The disturbance vectors W j (j = 1, 2, 3), columns of W in equation (3), are defined with their components being the coefficients of the partial derivatives in equation (2) Equation (2) denotes the compatibility relation if the disturbance vectors lie on a characteristic surface. 20 Hence, the parameters m 1 , m 2 , and m 3 are determined to have n 3 W = 0, in which n = ½n x , n y , n t denotes the outward normal to the characteristic surface. After multiplication and introducing x = n t + un x + vn y , one obtains x À un x À vn y n x n y b 2 n x x 0 b 2 n y 0 x To have a nontrivial solution, the determinant of the coefficient matrix must vanish, which yields where the following roots are found Each of x 0 , x + , and x À provides only one equation to determine the three unknown components of n, resulting in three envelopes of normals perpendicular to the characteristic surfaces, which are represented by the linear and quadratic factors of equation (5). Due to the homogeneity of these functions with respect to n, the corresponding characteristic velocities are found as 20 ð7 À bÞ Equation (7-a) shows the stream surface, which contains the pseudo velocity vector defined by V = (u, v, 1). The characteristic velocities l 6 x and l 6 y are the generators of two elliptic Mach conoids defined by equation (8). The characteristics and compatibility equations with the superscripts ( + ) and (2) correspond to the domains of dependence and influence, respectively 20 dx dt The compatibility equation on the stream surface is derived by finding m 1 , m 2 , and m 3 from equation (4) and substituting them in equation (2). Replacing x 0 = 0 in equation (4) reveals that the rank of the coefficient matrix is two. Hence, choosing m 2 = n y arbitrarily and solving equation (4) give m 1 = 0 and m 3 = À n x . By substituting m 1 , m 2 , and m 3 in equation (2) and introducing u 0 = tan À1 (n y =n x ) as the wave front angle, the compatibility equation on the stream surfaces is derived as equation (9) sinu 0 du dt In a similar way, substituting x 6 = (1=2)(a 6 c) in equation (4) indicates that one parameter should be specified arbitrary. Setting m 1 = À (a 6 c)=2 and solving equation (4) result in m 2 = b 2 cos u 6 and m 3 = b 2 sinu 6 . By substituting m 1 , m 2 , and m 3 in equation (2), the compatibility equation on the Mach conoid takes the following form b 2 cos u 6 du dt 6 + sin u 6 dv dt 6 À a 6 c 2 dp dt The terms Z 0 and Z 6 denote the cross derivatives, which include directional derivatives normal to the wave fronts. Equations (9) and (10) represent the twodimensional compatibility equations of the shear and pseudo-acoustic waves, respectively.

Characteristic discretization
Since a unique wave pattern does not exist in multidimensional flows, it is necessary to specify the wave model. The number of independent waves cannot exceed the number of governing equations. Here, one shear and two pseudo-acoustic waves are utilized. The wave fronts, u 0 and u + , are set to the grid outward normal direction, (u n ). The neighboring cells, which are not along bicharacteristic directions and specified by (j) in Figure 1, affect the solution point through the cross derivatives. Their influences on the interface values are excluded in the grid-aligned methods, while in the BCBM, they are taken into account, retaining the multidimensionality of the scheme.
The stencil of BCBM scheme is presented in Figure 1. The solution point, the vertex of the Mach cone, placed at the grid surface, where the interface values, p Ã , u Ã , and v Ã , should be determined. The pseudo path-line and two generators of the Mach cone are extended rearward from the solution point to intersect the initial value surface at the characteristic points 0, 1, and 2, respectively.
Setting u 0 0 = u n and discretizing equation (9) along the pseudo path-line result in equation (11-a). In a similar manner, equations (11-b) and (11-c) are derived by setting u + 1 = u n , u + 2 = u n + p, and discretizing equation (10) along the corresponding bicharacteristics Solving equation (11) for the interface values results in For a generic parameter u, standing for the flow variables and coefficients of equation (12), one obtains ð13-aÞ First À order accuracy : u L = u i, j , u R = u i + 1, j ð13-bÞ Second À order accuracy : in which k denotes the wave number. For the first-and second-order accuracy, the flow parameters, p k , u k , and v k , are computed by equations (13-b) and (13-c), respectively. A k and the velocity components of B k are estimated by the first-order accuracy, equation (13-b). The cross derivatives are computed at the upwind cells by the interface values obtained in the previous time level using the finite-volume method (i.e. equation (1)).

Extension to three-dimensional flows
With a similar approach implemented for the twodimensional compatibility relations, the threedimensional characteristics read: n y du À n x dv + Z 1 sh = 0, along l = V z 1 sh = dt n y p x À n x p y À Á ð14-aÞ n z du À n x dw + Z 2 sh = 0, along l = V z 2 sh = dt n z p x À n x p z ð Þ ð14-bÞ n x du + n y dv + n z dw À a + c 2b 2 dp + Z ac = 0 in which V = ½u, v, w and n = ½n x , n y , n z denote velocity vector and cell-face unit outward normal directions, repectively, and a = n Á V and c = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 2 + 4b 2 p . Equations (14-a) and (14-b) are two shear waves, corresponding to two independent tangential velocities in three dimensions. Equation (14-c) represents pseudoacoustic compatibility equation. The three-dimensional wave model includes two pseudo-acoustic waves with opposite propagation directions and two shear waves, which altogether provide four equations for computation of the four unknowns. The discretization is the same as the two-dimensional version.

Diffusive fluxes, time integration, and boundary conditions
The diffusive fluxes are computed by averaging on the secondary cells. 14 To update the flow field, equation (1) is discretized using the fifth-order explicit Runge-Kutta method. 3 The maximum time step is obtained by the CFL stability criterion based on the maximum characteristic velocity. At solid boundaries, the no-slip and the normal momentum conditions are imposed for the velocity and the pressure, respectively. At inlet far-field boundaries, the velocity is set to the free-stream value, and the pressure is extrapolated from the domain. At outlet, the pressure is fixed, and the velocity is extrapolated.

Results and discussion
The BCBM scheme is applied for solution of three test cases: the cross flow over a circular cylinder, the liddriven cavity flow, and the flow around sphere. In our numerical experiments, the fastest convergence is achieved with b 2 = 1. The BCBM and CBM are applied in the second-order form, unless otherwise stated. Convergence criteria are defined as u t j j\10 À9 . An Intel 2.2-GHz Dual-Core central processing unit (CPU) with 8-GB of RAM is used for the simulations.

Cross flow over a circular cylinder
A clustered 60 3 60 O-grid with the outer diameter of 20D (D denotes the cylinder diameter) is used for the solution of the steady cross flow around a circular cylinder. Nondimensionalization is carried out using the cylinder diameter and the free-stream values.
The pressure coefficients and the vorticity distribution along the cylinder surface are validated by the results of Fornberg, 22 as shown in Figures 2 and 3 respectively. The results show remarkable concurrence with the benchmark solutions. The streamlines at Re = 20 are plotted in Figure 4. In Table 1, the drag  coefficients and the length of vortices are compared to the benchmark results. Figure 5 shows the convergence histories at Re = 40. The permissible CFL numbers observed for the firstand second-order BCBM, CBM, and averaging method were 1.1, 1.2, 0.6, and 0.8, respectively. The second-order BCBM method converged approximately three times faster than the CBM. Furthermore, the convergence of the second-order BCBM is obtained with 20% less computational time as compared to the averaging method.

Lid-driven cavity flow
The governing equations are nondimensionalized using the velocity of the moving wall and the width of the cavity.
The numerical experiments at Re = 400, 1000, and 3200 are carried out on a uniform grid with 121 3 121 cells, while for Re = 5000, a finer grid with 201 3 201 cells is used.
The velocity profiles along the vertical and horizontal centerlines of the cavity are shown in Figures 6 and  7, respectively. The solutions obtained by the BCBM are in a good agreement with the benchmark results. 26,27 For example, at Re = 1000, the total deviations of u-velocity profiles from the benchmark results 26 are 0.63%, 2.31%, and 9.29% for the BCBM, averaging method, and second-order CBM, respectively. At Re = 5000, discrimination is more significance with the deviations of 1.87%, 10.28%, and 19.63% for the BCBM, averaging method, and secondorder CBM, respectively. The accuracy of the BCBM is compared to the first-, second-, and third-order CBM in Figure 8, which indicates the better performance of the BCBM.       To investigate the resolution of the BCBM, the streamlines at Re = 1000 are plotted in Figure 9. The vortex patterns obtained by the BCBM are similar to the previous reports. 27 The center locations and the sizes of the primary and secondary vortices are presented in Table 2. The maximum deviation from the results of Ghia et al. 26 is 4.55%.
The convergence history and the computational time of the first-and second-order BCBM at Re = 1000 are compared to the CBM and averaging method in Figure 10. The first-and second-order BCBM and averaging method are implemented with the CFL number equal to 1.6. The CBM shows the lower stability limit of the CFL number equal to 0.8. The second-order BCBM converges about seven times faster than the CBM.
The convergence rate of the BCBM is compared to the MCBM in Figure 11 for different grid sizes at Re = 1000. The comparison is made by the same convergence criteria (i.e. the normalized value of ∂p=∂t to be less than 10 À4 ) implemented in the work of Razavi et al. 13 The convergence of the BCBM is achieved after 1298 and 2595 iterations for 41 3 41 and 61 3 61 grids, respectively; while, for the MCBM, the solution converged after 5667 and 9582 iterations. 13 The faster convergence of the BCBM is mainly due to its more comprehensive wave model, which includes shear waves and pseudo-acoustic waves.
The effect of the artificial compressibility parameter on the convergence rate and accuracy of the scheme at Re = 1000 is represented in Figure 12. As expected, the accuracy of the solution is not affected by the artificial compressibility parameter. However, relatively faster convergence is achieved with b 2 = 1.

Flow around a sphere
To validate the three-dimensional scheme, the steady flow around a sphere is solved using a clustered O-shaped grid, with outer diameter of 20D (D denotes cylinder diameter), including 60 3 60 3 60 cells.
The vorticity and pressure distributions along the sphere surface at the symmetric plane, z = 0, are represented in Figure 13. All three applied methods show good agreement with the results of Rimon and Cheng. 28 However, the BCBM behaves more favorably. The isopressure surfaces and the contours of the pressure coefficient at Re = 100 are plotted in Figure  14, revealing significant concurrence with the results of Johnson and Patel. 29 The permissible CFL numbers observed for the BCBM, CBM, and averaging method were 1.4, 0.8, and 0.9, respectively. Figure 15 shows the convergence histories at Re = 100, revealing two times faster convergence of the BCBM as compared to the CBM.

Conclusion
In this study, the BCBM is proposed by implementation of the multidimensional characteristics of pseudo compressible flows for enhancing the convective fluxes. These equations contain directional derivatives along the bicharacteristics and normal to the wave fronts. The independence of these directions permits the incorporation of multidimensionality of information propagation in the flux computations. This approach is more consistent with the governing equation, as compared to MCBM, due to minimizing the influence of external controlling parameters such as data interpolations along parallel direction to the grid surface. In addition, implementation of the shear waves in the BCBM, which are excluded in the MCBM, provides more physical and more inclusive wave model, resulting in faster convergence of the BCBM. The proposed method     exploits higher resolution and faster convergence (up to seven times in some experiments) as compared to gridaligned CBM. These advantages will be very valuable in unsteady, turbulent, and three-dimensional flow simulations, where the reduction of the computational time is of major interst. 2

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.