Pseudo-characteristic upwind scheme for incompressible flows

Abstract

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.

Keywords

Incompressible flow multidimensional upwinding bicharacteristic-based upwinding artificial compressibility method Navier–Stokes equations

Introduction

The artificial compressibility method (ACM)¹ 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.² The resultant interrupted interaction of pressure field and viscous boundary layer evince inaccuracy and delayed convergence, particularly in separated and viscous internal flows.² 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.³ 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,⁴ heat transfer simulations,⁵ multigrid techniques,⁶ 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–14 Neofytou¹⁵ 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.¹³ 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 grid-aligned compressible upwind schemes.¹⁸ To incorporate the multidimensional physics of propagation of information, multidimensional characteristic-based method (MCBM) was developed for Cartesian,¹³ non-Cartesian,¹⁴ and unstructured grids.¹⁹ 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,²⁰ which results in compatibility relations with two directional derivatives along two independent directions. Since coefficient matrices along coordinate directions cannot be diagonalized simultaneously,²¹ 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,¹ in the finite-volume form can be delineated as equation (1)

\begin{matrix} \frac{d U}{dt} + \frac{1}{A_{ij}} \sum_{k} (F_{k} Δ y_{k} - G_{k} Δ x_{k}) - (R_{k} Δ y_{k} - S_{k} Δ x_{k}) = 0 \\ U = [\begin{matrix} p \\ u \\ v \end{matrix}], F = [\begin{matrix} β^{2} u \\ u^{2} + p \\ uv \end{matrix}], G = [\begin{matrix} β^{2} v \\ uv \\ v^{2} + p \end{matrix}], \\ R = \frac{1}{Re} [\begin{matrix} 0 \\ \frac{\partial u}{\partial x} \\ \frac{\partial v}{\partial x} \end{matrix}], S = \frac{1}{Re} [\begin{matrix} 0 \\ \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial y} \end{matrix}] \end{matrix}

(1)

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

\begin{matrix} m_{1} p_{t} + m_{2} p_{x} + m_{3} p_{y} + m_{2} u_{t} + (m_{1} β^{2} + m_{2} u) u_{x} \\ + m_{2} v u_{y} + m_{3} v_{t} + (m_{3} u) v_{x} + (m_{1} β^{2} + m_{3} v) v_{y} = 0 \end{matrix}

(2)

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)

W = [\begin{matrix} m_{2} β^{2} m_{1} + m_{2} u m_{3} u \\ m_{3} m_{2} v β^{2} m_{1} + m_{3} v \\ m_{1} m_{2} m_{3} \end{matrix}]

(3)

Equation (2) denotes the compatibility relation if the disturbance vectors lie on a characteristic surface.²⁰ Hence, the parameters $m_{1}, m_{2}, and m_{3}$ are determined to have $n \times W = 0$ , in which $n = [n_{x}, n_{y}, n_{t}]$ denotes the outward normal to the characteristic surface. After multiplication and introducing $χ = n_{t} + u n_{x} + v n_{y}$ , one obtains

[\begin{matrix} χ - u n_{x} - v n_{y} n_{x} n_{y} \\ β^{2} n_{x} χ 0 \\ β^{2} n_{y} 0 χ \end{matrix}] [\begin{matrix} m_{1} \\ m_{2} \\ m_{3} \end{matrix}] = 0

(4)

To have a nontrivial solution, the determinant of the coefficient matrix must vanish, which yields

χ (χ^{2} - α χ - β^{2} (n_{x}^{2} + n_{y}^{2})) = 0, α = u n_{x} + v n_{y}

(5)

where the following roots are found

χ^{0} = 0, χ^{\pm} = \frac{1}{2} (α \pm c), c = \sqrt{α^{2} + 4 β^{2} (n_{x}^{2} + n_{y}^{2})}

(6)

Each of $χ^{0}$ , $χ^{+}$ , and $χ^{-}$ 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²⁰

λ^{0} = ({\frac{dx}{dt}}^{0}, {\frac{dy}{dt}}^{0}) = (u, v)

(7-a)

\begin{matrix} λ^{\pm} = ({\frac{dx}{dt}}^{\pm}, {\frac{dy}{dt}}^{\pm}) \\ = (\frac{u (- α \pm c) - 4 β^{2} \cos φ^{\pm}}{\pm 2 c}, \frac{v (- α \pm c) - 4 β^{2} \sin φ^{\pm}}{\pm 2 c}) \end{matrix}

(7-b)

φ^{\pm} = \tan^{- 1} (\frac{n_{y}}{n_{x}}), α = u \cos φ^{\pm} + v \sin φ^{\pm}, a n d c = \sqrt{α^{2} + 4 β^{2}}

Equation (7-a) shows the stream surface, which contains the pseudo velocity vector defined by $V = (u, v, 1)$ . The characteristic velocities $λ_{x}^{\pm} and λ_{y}^{\pm}$ are the generators of two elliptic Mach conoids defined by equation (8). The characteristics and compatibility equations with the superscripts (+) and (−) correspond to the domains of dependence and influence, respectively²⁰

{({\frac{dx}{dt}}^{\pm} - \frac{u (- α \pm c)}{\pm 2 c})}^{2} + {({\frac{dy}{dt}}^{\pm} - \frac{v (- α \pm c)}{\pm 2 c})}^{2} = \frac{4 β^{4}}{\pm c^{2}}

(8)

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 $χ^{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 $φ^{0} = ta n^{- 1} (n_{y} / n_{x})$ as the wave front angle, the compatibility equation on the stream surfaces is derived as equation (9)

\begin{matrix} \sin φ^{0} {(\frac{du}{dt})}^{0} - \cos φ^{0} {(\frac{dv}{dt})}^{0} + Z^{0} = 0 \\ {(\frac{d}{dt})}^{0} = \frac{\partial}{\partial t} + λ_{x}^{0} \frac{\partial}{\partial x} + λ_{y}^{0} \frac{\partial}{\partial y}, Z^{0} = \sin φ^{0} p_{x} - \cos φ^{0} p_{y} \end{matrix}

(9)

In a similar way, substituting $χ^{\pm} = (1 / 2) (α \pm c)$ in equation (4) indicates that one parameter should be specified arbitrary. Setting $m_{1} = - (α \pm c) / 2$ and solving equation (4) result in $m_{2} = β^{2} \cos φ^{\pm} and m_{3} = β^{2} \sin φ^{\pm}$ . By substituting $m_{1}, m_{2}, and m_{3}$ in equation (2), the compatibility equation on the Mach conoid takes the following form

\begin{matrix} β^{2} (\cos φ^{\pm} {(\frac{du}{dt})}^{\pm} + \sin φ^{\pm} {(\frac{dv}{dt})}^{\pm}) \\ - \frac{α \pm c}{2} {(\frac{dp}{dt})}^{\pm} + Z^{\pm} = 0 \\ {(\frac{d}{dt})}^{\pm} = \frac{\partial}{\partial t} + λ_{x}^{\pm} \frac{\partial}{\partial x} + λ_{y}^{\pm} \frac{\partial}{\partial y} \\ Z^{\pm} = \frac{β^{2}}{\pm c} (u \sin φ^{\pm} - v \cos φ^{\pm}) \\ (\sin φ^{\pm} p_{x} - \cos φ^{\pm} p_{y}) + \frac{β^{2} (α \pm c)}{\pm 2 c} \\ [\cos φ^{\pm} (u \sin φ^{\pm} - v \cos φ^{\pm}) \mp c \sin φ^{\pm}] \\ (\sin φ^{\pm} u_{x} - \cos φ^{\pm} u_{y}) + \frac{β^{2} (α \pm c)}{\pm 2 c} \\ [\sin φ^{\pm} (u \sin φ^{\pm} - v \cos φ^{\pm}) \pm c \cos φ^{\pm}] \\ (\sin φ^{\pm} v_{x} - \cos φ^{\pm} v_{y}) \end{matrix}

(10)

The terms $Z^{0}$ and $Z^{\pm}$ denote the cross derivatives, which include directional derivatives normal to the wave fronts. Equations (9) and (10) represent the two-dimensional 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, $φ^{0}$ and $φ^{+}$ , are set to the grid outward normal direction, ( $φ_{n}$ ). The neighboring cells, which are not along bicharacteristic directions and specified by (▪) 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.

Figure 1.

Two-dimensional characteristics of incompressible flows and the stencil of the bicharacteristic-based method.

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 $φ_{0}^{0} = φ_{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 $φ_{1}^{+} = φ_{n}$ , $φ_{2}^{+} = φ_{n} + π$ , and discretizing equation (10) along the corresponding bicharacteristics

\sin φ_{n} (u^{*} - u_{0}) - \cos φ_{n} (v^{*} - v_{0}) + B_{0} = 0

(11-a)

\begin{matrix} \cos φ_{n} (u^{*} - u_{1}) + \sin φ_{n} (v^{*} - v_{1}) \\ - A_{1} (p^{*} - p_{1}) + B_{1} = 0 \end{matrix}

(11-b)

\begin{matrix} \cos φ_{n} (u^{*} - u_{2}) + \sin φ_{n} (v^{*} - v_{2}) \\ + A_{2} (p^{*} - p_{2}) - B_{2} = 0 \end{matrix}

(11-c)

B_{0} = Δ t {(Z^{0})}_{0}, A_{j} = \frac{α_{j} + c_{j}}{2 β^{2}}, B_{j} = \frac{Δ t}{β^{2}} {(Z^{+})}_{j}, j = 1, 2

Solving equation (11) for the interface values results in

p^{*} = \frac{1}{A_{1} + A_{2}} (A_{1} p_{1} + A_{2} p_{2} - \cos φ_{n} (u_{1} - u_{2}) - \sin φ_{n} (v_{1} - v_{2}) + B_{1} + B_{2})

(12-a)

\begin{matrix} u^{*} = \frac{\cos φ_{n}}{A_{1} + A_{2}} [A_{1} A_{2} (p_{2} - p_{1}) + \cos φ_{n} (A_{2} u_{1} + A_{1} u_{2}) \\ + \sin φ_{n} (A_{2} v_{1} + A_{1} v_{2}) - (A_{2} B_{1} - A_{1} B_{2}) \end{matrix}] + \sin φ_{n} (\sin φ_{n} u_{0} - \cos φ_{n} v_{0} - B_{0})

(12-b)

\begin{matrix} v^{*} = \frac{\sin φ_{n}}{A_{1} + A_{2}} [A_{1} A_{2} (p_{2} - p_{1}) + \cos φ_{n} (A_{2} u_{1} + A_{1} u_{2}) \\ + \sin φ_{n} (A_{2} v_{1} + A_{1} v_{2}) - (A_{2} B_{1} - A_{1} B_{2}) \end{matrix}] + \cos φ_{n} (\cos φ_{n} v_{0} - \sin φ_{n} u_{0} + B_{0})

(12-c)

For a generic parameter $θ$ , standing for the flow variables and coefficients of equation (12), one obtains

\begin{matrix} {(θ_{k})}_{i + \frac{1}{2}, j} = \frac{1}{2} [(1 + sign (γ_{k})) θ_{L} + (1 - sign (γ_{k})) θ_{R}] \\ k = 0, 1, 2 \end{matrix}

(13-a)

First - order accuracy : θ_{L} = θ_{i, j}, θ_{R} = θ_{i + 1, j}

(13-b)

\begin{matrix} Second - order accuracy : \\ θ_{L} = \frac{3}{2} θ_{i, j} - \frac{1}{2} θ_{i - 1, j}, θ_{R} = \frac{3}{2} θ_{i + 1, j} - \frac{1}{2} θ_{i + 2, j} \end{matrix}

(13-c)

γ_{k} = {(λ_{x})}_{k} \cos φ_{n} + {(λ_{y})}_{k} \sin φ_{n}

(13-d)

λ_{i + \frac{1}{2}, j} = \frac{1}{2} (λ_{i, j} + λ_{i + 1, j})

(13-e)

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 two-dimensional compatibility relations, the three-dimensional characteristics read:

\begin{matrix} n_{y} du - n_{x} dv + Z_{sh}^{1} = 0, along λ = V \\ z_{sh}^{1} = dt (n_{y} p_{x} - n_{x} p_{y}) \end{matrix}

(14-a)

\begin{matrix} n_{z} du - n_{x} dw + Z_{sh}^{2} = 0, along λ = V \\ z_{sh}^{2} = dt (n_{z} p_{x} - n_{x} p_{z}) \end{matrix}

(14-b)

\begin{array}{l} n_{x} d u + n_{y} d v + n_{z} d w - \frac{α + c}{2 β^{2}} d p + Z_{a c} = 0 a l o n g λ = \frac{(- α + c) V - 4 β^{2} n}{2 c} \\ Z_{a c} = \frac{(α + c) d t}{2 c} [n_{x} (V - (α - c) n) \cdot \nabla u + n_{y} (V - (α - c) n) \cdot \nabla v + n_{z} (V - (α - c) n) \cdot \nabla w - c \nabla \cdot V] + \frac{d t}{c} [(V - α n) \cdot \nabla p] \end{array}

(14-c)

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 $α = n \cdot V$ and $c = \sqrt{α^{2} + 4 β^{2}}$ .

Equations (14-a) and (14-b) are two shear waves, corresponding to two independent tangential velocities in three dimensions. Equation (14-c) represents pseudo-acoustic 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.¹⁴ To update the flow field, equation (1) is discretized using the fifth-order explicit Runge–Kutta method.³ 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 lid-driven cavity flow, and the flow around sphere. In our numerical experiments, the fastest convergence is achieved with $β^{2} = 1$ . The BCBM and CBM are applied in the second-order form, unless otherwise stated. Convergence criteria are defined as $| u_{t} | < 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 \times 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,²² 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 2.

Comparison of the pressure coefficients along the cylinder surface for the bicharacteristic-based, characteristic-based, and averaging methods and Fornberg:²² (a) Re = 20 and (b) Re = 40.

Figure 3.

Comparison of the span-wise vorticity distributions along the cylinder surface for the bicharacteristic-based, characteristic-based, and averaging methods and Fornberg:²² (a) Re = 20 and (b) Re = 40.

Figure 4.

Streamlines of the flow over a circular cylinder at Re = 20: (a) bicharacteristic-based, (b) averaging, and (c) characteristic-based methods.

Table 1.

Drag coefficients (C_d ) and length of vortices (L_w ) for the flow around a circular cylinder.

	Re = 20		Re = 40
	$C_{d}$	$\frac{L_{w}}{D}$	$C_{d}$	$\frac{L_{w}}{D}$
BCBM	2.05	0.93	1.53	2.25
CBM	2.10	0.98	1.59	2.41
Averaging	2.06	0.93	1.54	2.27
Fornberg²²	2.00	0.91	1.50	2.24
Dennis and Chang²³	2.05	0.94	1.52	2.35
Calhoun²⁴	2.19	0.91	1.62	2.18
Russell and Jane Wang²⁵	2.13	0.94	1.60	2.29

BCBM: bicharacteristic-based method; CBM: characteristic-based method.

Figure 5 shows the convergence histories at Re = 40. The permissible CFL numbers observed for the first- and 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.

Figure 5.

Convergence histories of the bicharacteristic-based, characteristic-based, and averaging methods for the flow around a circular cylinder at Re = 40.

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 \times 121$ cells, while for $Re = 5000$ , a finer grid with $201 \times 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²⁶ 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 second-order 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.

Figure 6.

u-velocity profiles along the vertical centerline of the cavity for the bicharacteristic-based, characteristic-based, and averaging methods and benchmark results:^26,27 (a) Re = 400, (b) Re = 1000, (c) Re = 3200, and (d) Re = 5000.

Figure 7.

v-velocity profiles along the horizontal centerline of the cavity for the bicharacteristic-based, characteristic-based, and averaging methods and benchmark results:^26,27 (a) Re = 400, (b) Re = 1000, (c) Re = 3200, and (d) Re = 5000.

Figure 8.

Comparison of the (a) u-velocity and (b) v-velocity profiles for the bicharacteristic-based method, characteristic-based methods, and Ghia et al.²⁶ at Re = 1000.

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.²⁷ 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.²⁶ is 4.55%.

Figure 9.

Streamline plots at Re = 1000: (a) bicharacteristic-based, (b) averaging, (c) Erturk et al.,²⁷ (d) first-order characteristic-based, and (e) second-order characteristic-based methods.

Table 2.

Center location and size of the primary and secondary vortices for the cavity flow.

		Primary	Bottom-right		Bottom-left		Top-left
		x; y	x; y	h; v	x; y	h; v	x; y	h; v
Re = 400	BCBM	0.5557; 0.6083	0.8924; 0.1269	0.2602; 0.3244	0.0491; 0.0458	0.1291; 0.1101	–	–
	Ghia et al.²⁶	0.5547; 0.6055	0.8906; 0.1250	0.2617; 0.3203	0.0508; 0.0469	0.1273; 0.1081	–	–
	Deviation	0.18%; 0.46%	0.20%; 1.52%	0.57%; 1.28%	3.35%; 2.34%	1.41%; 1.66%	–	–
Re = 1000	BCBM	0.5319; 0.5674	0.8685; 0.1131	0.2982; 0.3620	0.0822; 0.0759	0.2223; 0.1682	–	–
	Ghia et al.²⁶	0.5313; 0.5625	0.8594; 0.1094	0.3034; 0.3536	0.0859; 0.078	0.2188; 0.1680	–	–
	Deviation	0.09%; 0.87%	0.74%; 3.38%	1.70%; 2.37%	4.30%; 2.81%	1.59%; 0.11%	–	–
Re = 3200	BCBM	0.5185; 0.5421	0.8237; 0.0898	0.3425; 0.4117	0.0841; 0.1131	0.2948; 0.2395	0.0525; 0.8933	0.0825; 0.2089
	Ghia et al.²⁶	0.5165; 0.5469	0.8125; 0.0859	0.3406; 0.4102	0.0859; 0.1094	0.2844; 0.2305	0.0547; 0.8984	0.0859; 0.2057
	Deviation	0.38%; 0.87%	1.37%; 4.55%	0.54%; 0.36%	2.00%; 3.38%	3.60%; 3.90%	4.02%; 0.56%	3.95%; 1.55%
Re = 5000	BCBM	0.5159; 0.5362	0.8134; 0.0758	0.3505; 0.4302	0.0717; 0.1397	0.3308; 0.2643	0.0633; 0.9074	0.1189; 0.2660
	Ghia et al.²⁶	0.5117; 0.5352	0.8086; 0.0742	0.3565; 0.4180	0.0703; 0.1367	0.3184; 0.2643	0.0625; 0.9102	0.1211; 0.2693
	Deviation	0.82%; 0.18%	0.59%; 2.15%	1.68%; 2.91%	1.99%; 2.19%	3.89%; 0.00%	1.28%; 0.30%	1.81%; 1.22%

BCBM: bicharacteristic-based method.

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.

Figure 10.

Convergence behavior of the bicharacteristic-based, averaging, and characteristic-based methods in the cavity flow at Re = 1000.

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 $\partial p / \partial t$ to be less than $10^{- 4}$ ) implemented in the work of Razavi et al.¹³ The convergence of the BCBM is achieved after 1298 and 2595 iterations for $41 \times 41$ and $61 \times 61$ grids, respectively; while, for the MCBM, the solution converged after 5667 and 9582 iterations.¹³ The faster convergence of the BCBM is mainly due to its more comprehensive wave model, which includes shear waves and pseudo-acoustic waves.

Figure 11.

Comparison of the convergence rates of the bicharacteristic-based method for different grid sizes.

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 $β^{2} = 1$ .

Figure 12.

Effect of the artificial compressibility parameter on the convergence and accuracy of the bicharacteristic-based method.

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 \times 60 \times 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.²⁸ 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.²⁹

Figure 13.

Comparison of the pressure coefficient (left) and vorticity (right) profiles along the sphere surface at z = 0.5 with the benchmark results:²⁹ (a) Re = 40 and (b) Re = 100.

Figure 14.

Pressure field of flow around a sphere at Re = 100 at z = 0: (a) pressure coefficient contours of bicharacteristic-based method, (b) pressure coefficient contours of Johnson and Patel,²⁹ and (c) the isopressure surfaces of bicharacteristic-based method.

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.

Figure 15.

Comparison of the convergence behavior of the bicharacteristic-based, averaging, and characteristic-based methods at Re = 100 for the flow around a sphere.

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 grid-aligned 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.²

Footnotes

Academic Editor: Roslinda Nazar

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.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Chorin

. A numerical method for solving incompressible viscous flow problems. J Comput Phys 1967; 2: 12–26.

Kwak

Kiris

Kim

. Computational challenges of viscous incompressible flows. Comput Fluids 2005; 34: 283–299.

Drikakis

Govatsos

Papantonis

. A characteristic-based method for incompressible flows. Int J Numer Meth Fl 1994; 19: 667–685.

Drikakis

. A parallel multiblock characteristic-based method for three-dimensional incompressible flows. Adv Eng Softw 1996; 26: 111–119.

Zhao

Zhang

. A high-order characteristics upwind FV method for incompressible flow and heat transfer simulation on unstructured grids. Comput Method Appl M 2000; 190: 733–756.

Drikakis

Iliev

Vassileva

. A nonlinear multigrid method for the three-dimensional incompressible Navier–Stokes equations. J Comput Phys 1998; 146: 301–321.

Tai

Zhao

. Parallel unsteady incompressible viscous flow computations using an unstructured multigrid method. J Comput Phys 2003; 192: 277–311.

Shapiro

Drikakis

. Artificial compressibility, characteristics-based schemes for variable density, incompressible, multi-species flows. Part I. Derivation of different formulations and constant density limit. J Comput Phys 2005; 210: 584–607.

Shapiro

Drikakis

. Artificial compressibility, characteristics-based schemes for variable-density, incompressible, multispecies flows: part II. Multigrid implementation and numerical tests. J Comput Phys 2005; 210: 608–631.

10.

Nithiarasu

Liu

C-B

. An artificial compressibility based characteristic based split (CBS) scheme for steady and unsteady turbulent incompressible flows. Comput Method Appl M 2006; 195: 2961–2982.

11.

Atashbar Orang

. Genuinely characteristic-based scheme for the incompressible turbulent flows. Comput Fluids 2014; 103: 175–185.

12.

Neofytou

Drikakis

. Non-Newtonian flow instability in a channel with a sudden expansion. J Non-Newton Fluid 2003; 111: 127–150.

13.

Razavi

Zamzamian

Farzadi

. Genuinely multidimensional characteristic-based scheme for incompressible flows. Int J Numer Meth Fl 2008; 57: 929–949.

14.

Zamzamian

Razavi

. Multidimensional upwinding for incompressible flows based on characteristics. J Comput Phys 2008; 227: 8699–8713.

15.

Neofytou

. Revision of the characteristics-based scheme for incompressible flows. J Comput Phys 2007; 222: 475–484.

16.

Zhao

Huang

. On the characteristics-based ACM for incompressible flows. J Comput Phys 2007; 227: 1–11.

17.

Razavi

Orang

. A comparative investigation of hydrofoil at angles of attack. Int J Numer Meth Fl 2012; 68: 1087–1101.

18.

Levy

Powell

van Leer

. Use of a rotated Riemann solver for the two-dimensional Euler equations. J Comput Phys 1993; 106: 201–214.

19.

Hashemi

Zamzamian

. A multidimensional characteristic-based method for making incompressible flow calculations on unstructured grids. J Comput Appl Math 2014; 259(Part B): 752–759.

20.

Chapman

. High speed flow. 1st ed. New York: Cambridge University Press, 2000.

21.

Deconinck

Hirsch

Peuteman

. Characteristic decomposition methods for multidimensional Euler equations. In: Zhuang

Zhu

(eds) Tenth international conference on numerical methods in fluid dynamics, vol. 264 (Lecture notes in physics). Berlin: Springer, 1986, pp.216–221.

22.

Fornberg

. A numerical study of steady viscous flow past a circular cylinder. J Fluid Mech 1980; 98: 819–855.

23.

Dennis

Chang

G-Z

. Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100. J Fluid Mech 1970; 42: 471–489.

24.

Calhoun

. A Cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions. J Comput Phys 2002; 176: 231–275.

25.

Russell

Jane Wang

. A Cartesian grid method for modeling multiple moving objects in 2D incompressible viscous flow. J Comput Phys 2003; 191: 177–205.

26.

Ghia

Shin

. High-resolutions for incompressible flow using the Navier-Stokes equations and a multigrid method. J Comput Phys 1982; 48: 387–411.

27.

Erturk

Corke

Gökçöl

. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers. Int J Numer Meth Fl 2005; 48: 747–774.

28.

Rimon

Cheng

. Numerical solution of a uniform flow over a sphere at intermediate Reynolds numbers. Phys Fluids 1969; 12: 949–959.

29.

Johnson

Patel

. Flow past a sphere up to a Reynolds number of 300. J Fluid Mech 1999; 378: 19–70.