Divergence-Free Pressure Poisson Equation for Viscous Incompressible Flows with Solid Wall Boundaries

2.1. PPE Reformulations

In conventional PPE solution procedures, the pressure field is uniquely determined by the velocity field of the same moment, independent of the flow’s time history [17]. Some researchers have pointed out that this characteristic constitutes an advantage facilitating the realization of high-order solutions in time [18]. For incompressible flows, the pressure is also considered a Lagrange multiplier to satisfy the continuity constraint, rendering the pressure merely an auxiliary function to correct the velocity field. While these observations may be correct in mathematics, they do not contribute effectively to the proper solution of the pressure field. Orszag et al. [20] have criticized the “parabolic” PPE scheme, saying that it forgets the effect of initial conditions immediately and only reflects the effect of boundary conditions at the domain edge. In his review paper [5], Rempfer has proposed that “the value of the pressure at any point in space and time does not only depend on the value of the velocity field everywhere in the domain; it also depends on its history in time.” As discussed in the last chapter, the boundary-condition-free requirement for PPE necessitates the pressure field being corrected from the one in the last iteration.

Inspired by these thoughts, we aim to reformulate the PPE and the Neumann boundary conditions in this study to incorporate the pressure field from the preceding iteration. The momentum equation for viscous incompressible flows in conservation form is

$${\displaystyle \frac{\partial \mathit{u}}{\partial t}}+\nabla \cdot (\mathit{u}\mathit{u})+\nabla P/\rho =\nu {\nabla }^2\mathit{u}+\mathit{f}\hspace{1em}\hspace{1em}\hspace{1em}t \ge 0, \mathit{x}\in \Omega$$

(3)

in which the velocity u = (u(x, y, z), v(x, y, z), w(x, y, z)), ν being the kinematic viscosity and f the body force per unit mass acting on the fluid. Taking the divergence of (3) with the velocity dilatation rate

$$\mathsf{\Phi }=\nabla \cdot \mathit{u}$$

(4)

we can obtain (5) without assuming Ф = 0

$${\nabla }^{2}P/\rho =-(\nabla \mathit{u}):{(\nabla \mathit{u})}^T+\nabla \cdot \mathit{f}+\nu {\nabla }^2\mathsf{\Phi }-{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial t}}-\mathit{u}\cdot \nabla \mathsf{\Phi }-{\mathsf{\Phi }}^2\hspace{1em}\hspace{1em}t \ge 0, \mathit{x}\in \Omega$$

(5)

With nonzero Ф, (5) can be considered a PPE of the velocity field, which is not converged and solenoidal. When the divergence-free condition is satisfied, we should have Ф = 0 and ∂Ф/∂t = 0 throughout the domain, and (5) becomes

$${\nabla }^{2}P/\rho =-(\nabla \mathit{u}):{(\nabla \mathit{u})}^T+\nabla \cdot \mathit{f}$$

(6)

(6) is the so-called standard PPE documented in the literature and can be considered the PPE for the solenoidal velocity field. By comparing (5) and (6), we notice that the velocity source term –(▽u):(▽u)^Tand the body force term ▽·f are the common components in both equations. The two PPE only differ by the Ф terms. By eliminating the common terms on their RHS, (5) and (6) suggest a revised PPE in an iterative scheme, with k denoting the iteration index and n+1 referring to the time level at t+∆t.

$${\nabla }^2{P}^{\mathit{n}+1,k+1}/\rho ={\nabla }^2{P}^{\mathit{n}+1,k}/\rho -\nu {\nabla }^2{\mathsf{\Phi }}^{\mathit{n}+1,k}+{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial t}}+{\mathit{u}}^{\mathit{n}+1,k}\cdot \nabla {\mathsf{\Phi }}^{\mathit{n}+1,k}+{({\mathsf{\Phi }}^{\mathit{n}+1,k})}^2\hspace{1em}\hspace{1em}t \ge 0, \mathit{x}\in \Omega$$

(7a)

For steady-state computations using a time-marching approach with unity density, (7a) is simplified as

$${\nabla }^2{P}^{\mathit{n}+1}={\nabla }^2{P}^{\mathit{n}}-\nu {\nabla }^2{\mathsf{\Phi }}^{\mathit{n}}+{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial t}}+{\mathit{u}}^{\mathit{n}}\cdot \nabla {\mathsf{\Phi }}^{\mathit{n}}+{({\mathsf{\Phi }}^{\mathit{n}})}^2\hspace{1em}\hspace{1em}t \ge 0, \mathit{x}\in \Omega$$

(7b)

The viscous term in (7b) conserves the velocity diffusion. Though it produces negligible difference in the solutions, it can boost numerical stability and permit larger pseudo time steps in the computation. The Ф convective term in (7b) is necessary for convection-dominated flows when PPE is directly coupled with the momentum equations. The temporal Ф term in (7b) is indispensable and can be approximated by

$${\displaystyle \frac{\partial \mathsf{\Phi }}{\partial t}}={\displaystyle \frac{{\mathsf{\Phi }}^{\mathit{n}+1}-{\mathsf{\Phi }}^{\mathit{n}}}{\mathsf{\Delta }t}}$$

(8)

In standard PPE practice since its introduction by Harlow and Welch [19], Фⁿ⁺¹ is set to zero to enforce continuity, while Фⁿ is retained as nonzero to maintain numerical stability. When additional Ф terms are involved in (7b), our experience shows that three approaches are viable for handling the temporal Ф terms in (8): retaining both terms, removing both, or retaining only Фⁿ⁺¹. Notably, removing both terms leads to slow convergence with occasional instability, whereas retaining only Фⁿ⁺¹ yields the fastest convergence. And the diffusion term can be evaluated at either time level n or n+1, depending on the computational needs.

Compared with the conventional PPEs of (5) or (6), the new PPE (7) includes only Ф terms on its RHS, eliminating the need to calculate the velocity source term and the divergence of the body force f. The influence of f on pressure is instead accounted for in the momentum equation calculations. The iterative nature of (7) establishes a mechanism that drives the Ф field towards zero during the pressure evolution. Once Ф converges to a trivial value, (7) becomes equivalent to the continuity equation. The use of the variable Ф avoids the direct coupling between the pressure and the velocity, and moderates the pressure updating through its spatial distribution, transport, and diffusion. A similar technique was employed in the pressure gradient method study [21] to achieve divergence-free solutions without the need for specifying pressure gradient boundary conditions. Since |Ф| << 1 typically holds in the domain interior, the nonlinear terms and viscous terms in (7b) are generally small in magnitude. Defining P’ = Pⁿ⁺¹– Pⁿ, (7b) retaining only the temporal term, is simplified as

$${\nabla }^{2}P\text{'}={\displaystyle \frac{\partial \mathsf{\Phi }}{\partial t}}\approx {\displaystyle \frac{\mathsf{\Phi }}{\mathsf{\Delta }t}}\hspace{1em}\hspace{1em}t \ge 0, \mathit{x}\in \Omega$$

(7c)

This approximates the corrective PPE used in the well-known SIMPLE method, in which zero velocity divergence is assumed at the next iteration and Ф represents the velocity dilatation of the intermediate velocity field. Therefore, the theoretical soundness of (7b) is partially validated by the proven success of the SIMPLE algorithm. And the viscous terms in (7b) can be extended to the SIMPLE series method and the projection methods to improve numerical stability. However, it should be noted that (7c) can only be solved separately from the momentum equation as in the SIMPLE algorithm, due to its lack of the Ф convection term required for mass transport.

2.2. Corrective Pressure Boundary Condition

Boundary condition specification for the pressure is required for solving (7), particularly the Neumann boundary condition for wall-bounded flows. As discussed in Section I, directly applying the conventional boundary condition of (2) will not yield a divergence-free solution, though it is theoretically consistent with the Navier-Stokes equations. We rederive the pressure gradient expression from the momentum equation (3) with all Ф terms retained as

$${\left.{\displaystyle \frac{\partial P}{\partial \mathit{n}}}\right|}_{\partial \mathsf{\Omega }}=\rho \mathit{n}\cdot [-\nabla \cdot (\mathit{u}\mathit{u})+\nu {\nabla }^2\mathit{u}]{|}_{\partial \mathsf{\Omega }}=\rho \mathit{n}\cdot [-\mathit{u}\mathsf{\Phi }-(\mathit{u}\cdot \nabla )\mathit{u}+\nu (-\nabla \times \nabla \times \mathit{u}+\nabla \mathsf{\Phi })]{|}_{\partial \Omega }$$

(9a)

or

$${\left.{\displaystyle \frac{\partial P}{\partial \mathit{n}}}\right|}_{\partial \Omega }={\left.\rho [-{\mathit{u}}_{\mathit{n}}\mathsf{\Phi }-({\mathit{u}}_{\mathit{n}}\cdot \nabla ){\mathit{u}}_{\mathit{n}}+\nu (-\nabla \times \omega +{\displaystyle \frac{\partial \mathsf{\Phi }}{\partial \mathit{n}}})]\right|}_{\partial \Omega }$$

(9b)

in which ω=▽◊u is the curl. (9b) can be interpreted as the Neumann pressure boundary condition when Ф≠0 and the flow field is non-solenoidal prior to convergence. When Ф=0, the above equation corresponds to the converged and solenoidal flow field, taking the form:

$${\left.{\displaystyle \frac{\partial P}{\partial \mathit{n}}}\right|}_{\partial \Omega }={\left.\rho [-({\mathit{u}}_{\mathit{n}}\cdot \nabla ){\mathit{u}}_{\mathit{n}}+\nu (-\nabla \times \omega )]\right|}_{\partial \Omega }$$

(9c)

Following the same methodology used in deriving (7), we combine (9b) and (9c) and obtain an evolving Neumann boundary condition in an iterative scheme with time marching.

$${\left.{\displaystyle \frac{\partial P^{\mathit{n}+1}}{\partial \mathit{n}}}\right|}_{\partial \Omega }={\left.{\displaystyle \frac{\partial P^{\mathit{n}}}{\partial \mathit{n}}}\right|}_{\partial \Omega }+{\left.\rho ({\mathit{u}}_{\mathit{n}}{\mathsf{\Phi }}^{\mathit{n}}-\nu {\displaystyle \frac{\partial {\mathsf{\Phi }}^{\mathit{n}}}{\partial \mathit{n}}})\right|}_{\partial \Omega }$$

(10)

Compared with boundary conditions (1) and (2), (10) involves only 1st-order derivatives, incorporates the tangential velocity effects through the Ф term, and imposes no larger complexity in implementation. It is evident that this boundary constraint will be weakened as Ф gets close to zero globally, making the system well-posed. This is the same principle as the artificial compressibility method and (7): the boundary pressure update will diminish when Ф is approaching zero. (10) also accounts for the influence of the wall normal velocity via the u_nФ term. With respect to the numerical compatibility, we have found that this boundary condition (10) alone will not produce a divergence-free solution without coupling with the new PPE (7), as will be shown in Section IV. For implementations, the normal derivatives of P and Ф can utilize one-sided spatial discretization schemes of any high order. This new pressure boundary condition can also be promoted to the SIMPLE family to achieve divergence-free solutions and to the artificial compressibility method. Orszag et al. [20] have proposed an iterative Neumann pressure boundary condition in a similar form to (10) for high-order temporal accuracy, without providing any derivations and numerical validation. But their formulation features an undetermined parameter and misses the u_n term for wall normal speed.

Since the calculation of (7) and (10) needs the Ф values at the boundaries, the boundary condition for Ф is equally critical. The boundary Ф specification must be numerically consistent with the interior Ф calculation to make the Ф field continuous and smooth in the computation. For wall-bounded flows, setting Ф directly to zero at wall boundaries skips the non-solenoidal evolution process and may fail to yield a divergence-free solution. A Rubin boundary condition for Ф has been proposed in [22] using the energy method. It is equivalent to setting ∂Ф/∂n = 0 at the stationary solid wall boundary. It is not compatible with (10), because ∂Ф/∂n = 0 implies a constant ∂P/∂n in (10) from its initial value at the boundaries. Therefore, we adopt the Dirichlet boundary condition for Ф on solid walls and compute it with (4) and one-sided difference scheme.

On a staggered mesh, the computation of the Ф field should be simple and straightforward. However, for collocated mesh, the pressure and Ф field solutions have been observed in numerical experiments with some spurious oscillations near certain solid walls. They stem from the numerical discontinuities between the interior and the boundary Ф values, introduced by the artificial viscosity term used to supplement the central scheme for interior Ф computations. When the magnitude of Ф diminishes, any non-decaying numerical errors become significant, contaminating the domain and causing solution wiggles. As we know, an upwind or one-sided scheme can be decomposed into one symmetrical part and one antisymmetric part. Conversely, the antisymmetric central scheme augmented with a symmetric artificial viscosity term is equivalent to an upwind or one-sided scheme. For illustration, consider a finite difference central scheme (CS) with uniform grid spacing ∆x, defined as

$$CS={\displaystyle \frac{u_{i-2}-8u_{i-1}+8u_{i+1}-u_{i+2}}{12\mathsf{\Delta }x}}+O{(\mathsf{\Delta }x)}^4$$

(11a)

and a typical 4th-order artificial viscosity term (AVT) as

$$AVT={\displaystyle \frac{u_{i-2}-4u_{i-1}+6u_i-4u_{i+1}+u_{i+2}}{12\mathsf{\Delta }x}}+O{(\mathsf{\Delta }x)}^4$$

(11b)

As shown in (12), a 3rd-order forward/backward one-sided scheme (used for the left/right boundary) can be expressed as a linear combination of the above CS and AVT terms. In practice, however, these combinations typically involve terms of different orders and include a scaling coefficient.

$${\left.{\displaystyle \frac{\partial u}{\partial x}}\right|}_i\approx {\displaystyle \frac{-2u_{i-1}-3u_i+6u_{i+1}-u_{i+2}}{6\mathsf{\Delta }x}}=CS-AVT$$

(12a)

$${\left.{\displaystyle \frac{\partial u}{\partial x}}\right|}_i\approx {\displaystyle \frac{u_{i-2}-6u_{i-1}+3u_i+2u_{i+1}}{6\mathsf{\Delta }x}}=CS+AVT$$

(12b)

The backward one-sided scheme (12b) for boundaries with ending coordinate indices equals CS + AVT. This is the same combination used in the domain interior to compute ∂u/∂x for Ф. In contrast, the forward one-sided scheme (12a) for boundaries with starting coordinate indices conflicts with CS + AVT, being responsible for the oscillatory Ф field near these boundaries. Though the numerical schemes (including compact schemes) applied in implementations may be different, the core inconsistency between boundary and interior discretizations persists. The artificial dissipation term can be gradually removed with an exponentially decaying coefficient as performed by Rogers et al. [23], but the decaying parameter requires some case-by-case tuning, making this technique not very feasible.

We have found a simple technique to resolve this inconsistency. As pointed out by Mazumder [24], in a physical problem, when a boundary derivative (e.g., ∂u_n/∂n) is not a prescribed value but computed with one-sided schemes, additional numerical error other than the truncation error will arise. And the solution of this derivative will be discontinuous from the interior to the boundaries, because the governing equation is not satisfied at the boundary. Mazumder has provided a corrective treatment: enforce the governing equation in the one-sided scheme for 1st-order derivatives by incorporating a 2nd-order derivative, which may be already known from the physical problem. For one-dimensional implementation with uniform grid spacing ∆x, a 2^nd-order derivative can be embedded into the one-sided difference scheme for the 1st-order u derivative by cancelling the 3^rd-order derivatives in the following two Taylor series expansions:

$$u_{i+1}=u_i+(\mathsf{\Delta }x){\left.{\displaystyle \frac{du}{dx}}\right|}_i+{\displaystyle \frac{{(\mathsf{\Delta }x)}^2}{2!}}{\left.{\displaystyle \frac{d^{2}u}{d{x}^2}}\right|}_i+{\displaystyle \frac{{(\mathsf{\Delta }x)}^3}{3!}}{\left.{\displaystyle \frac{d^{3}u}{d{x}^3}}\right|}_i+{\displaystyle \frac{{(\mathsf{\Delta }x)}^4}{4!}}{\left.{\displaystyle \frac{du}{dx}}\right|}_i+\mathrm{...}$$

(13a)

$$u_{i+2}=u_i+(2\mathsf{\Delta }x){\left.{\displaystyle \frac{du}{dx}}\right|}_i+{\displaystyle \frac{{(2\mathsf{\Delta }x)}^2}{2!}}{\left.{\displaystyle \frac{d^{2}u}{d{x}^2}}\right|}_i+{\displaystyle \frac{{(2\mathsf{\Delta }x)}^3}{3!}}{\left.{\displaystyle \frac{d^{3}u}{d{x}^3}}\right|}_i+{\displaystyle \frac{{(2\mathsf{\Delta }x)}^4}{4!}}{\left.{\displaystyle \frac{du}{dx}}\right|}_i+\mathrm{...}$$

(13b)

so that

$${\left.{\displaystyle \frac{du}{dx}}\right|}_i={\displaystyle \frac{-7u_i+8u_{i+1}-u_{i+2}-2{(\mathsf{\Delta }x)}^2{\left.(d^{2}u/d{x}^2)\right|}_i}{6\mathsf{\Delta }x}}+O{(\mathsf{\Delta }x)}^4$$

(14)

In computing ∂u_n/∂n for (4), the 2nd-order term ∂²u_n/∂n² is accessible in the momentum equation in the form of (2) and can be replaced with (∂P/∂n)/μ for a solid wall boundary. Therefore, Ф at the wall boundary with the starting coordinate index is computed as

$${\left.\mathsf{\Phi }{|}_i={\displaystyle \frac{\partial {\mathit{u}}_{\mathit{n}}}{\partial \mathit{n}}}\right|}_i={\displaystyle \frac{-7{\mathit{u}}_{\mathit{n},i}+8{\mathit{u}}_{\mathit{n},i+1}-{\mathit{u}}_{\mathit{n},i+2}-2{(\mathsf{\Delta }x)}^2{\left.(\partial P/\partial \mathit{n})\right|}_i/\mu }{6\mathsf{\Delta }x}}$$

(15a)

By replacing the d²u/dx² term with the pressure gradient, (15a) establishes a new local coupling that improves numerical accuracy and convergence rate. Though (2) should not be used as a proper boundary condition to solve PPE, it is still a useful relation in this treatment. For finite volume implementations, in which the boundary conditions are applied on the cell faces rather than the cell centers, ∂u_n/∂n is computed differently. On a one-dimensional grid of equal cell spacing ∆x, Ф at the left (W) wall boundaries (starting coordinate index) is calculated as

$${\left.{\mathsf{\Phi }}_W={\displaystyle \frac{\partial {\mathit{u}}_{\mathit{n}}}{\partial \mathit{n}}}\right|}_W={\displaystyle \frac{104{\mathit{u}}_{\mathit{n},W}-108{\mathit{u}}_{\mathit{n},O}+4{\mathit{u}}_{\mathit{n},E}-9{(\mathsf{\Delta }x)}^2{\left.(\partial P/\partial \mathit{n})\right|}_W/\mu }{48\mathsf{\Delta }x}}$$

(15b)

with the centers of the boundary cell and the cell next to it denoted as O and E.

Based on the above discussions, we choose to compute Ф with (15) at the boundaries of starting coordinate indices, while employing uncorrected one-sided schemes to calculate Ф at the solid wall boundaries of ending indices. The uncorrected one-sided schemes can be of any order relative to the central scheme for computing the interior Ф. With this implementation on collocated grids, the pressure and Ф fields become very smooth without any wiggles, and we can minimize the Ф magnitude to nearly machine accuracy. Although (15) will turn out to be a 3rd-order one-sided scheme if the d²u/dx² term is substituted with a one-sided discretization, this 3rd-order scheme will not work as effectively as (15). And if we apply the correction procedure to the opposite one-sided scheme for wall boundaries of ending indices, the Ф solution will be a uniform nontrivial value, since the correction disrupts the prior CS + AVT consistency.

For non-solid wall boundaries, such as an inflow with a static pressure value and a domain symmetry with zero pressure gradient, the theoretical requirement of avoiding explicit pressure boundary conditions is violated. In these cases, the solution of the PPE (7) should be restricted to the domain interior and solid wall boundaries. For inflow or outflow boundaries, Ф can be set to zero, assuming the far-upstream or far-downstream flows satisfy the continuity as in fully developed or undisturbed regions. This treatment maintains consistency with the iterative evolution of Φ toward zero in the domain interior.

4.1. Lid-Driven Cavity Flow

The first validation case is the classic top-lid-driven cavity flow. This wall-bounded flow is simple in geometry but rich in flow features. The square domain consists of boundary conditions of all Neumann type for the pressure. It also features numerical singularities at the two upper corners, which often pose difficulties for conventional PPE solutions. This problem is first tested using our proposed formulations on a 129 × 129 mesh. The obtained centerline velocity profile results for different Reynolds numbers are plotted in Figure 1 and compared with the benchmark data of Ghia et al. [31], which are based on the vorticity-stream function method and strictly divergence-free. The two data sets show excellent agreement. Figure 2(a) presents the convergence history for the case with Re = 1000. Figrures 2(b) and 2(c) demonstrate the domain solutions of the velocity dilatation, confirming the capability of our new formulations for producing a genuinely divergence-free velocity field. The resolved magnitude level of Ф is independent of the grid spacing. And it is not affected by the corner singularity from the domain geometry. We have also noticed that before Ф converges to the level of nearly machine zero, its contours bear a rough resemblance to the pressure contours. This Ф-P quantitative relationship can be described by the penalty method. Although the artificial viscosity involved in computing (4) may be nontrivial in certain regions, it does not compromise the fact that this formulation is able to drive the given Ф expression to near machine zero.

We then conducted a comparative study in numerics with this cavity flow problem on a 101× 101 mesh. Table 1 elaborates the case setup employing the PPE and pressure boundary conditions from the proposed method, previous methods, and their mixtures. We also examine the curl-curl and the homogeneous pressure boundary conditions in this study. All cases utilize identical formulations except for the PPE and pressure boundary condition. For Cases 2 and 3, the discretization of (2a) includes ghost points computed by applying the continuity equation at the boundary. For Case 4, the boundary vorticity is calculated for the curl-curl boundary condition (2b). Table 1 also lists the maximum stable time step sizes in serial computing for comparison, with equal underrelaxations in variable updating for each case. The well-converged velocity divergence solutions for Case 1~6 are presented in Figrures 3(a)~(f). It is evident that our proposed PPE and pressure boundary condition outperform the others and represent the only combination capable of achieving divergence-free solutions (Case 1). The rest of the cases exhibit high magnitudes of velocity dilatation up to the level of ~10² near the two upper corners, which exceed typical numerical errors and indicate poor solution quality under the influence of boundary singularities. To better visualize the solution differences, the contour range for Case 2 to 6 is restricted to -0.1 ~ 0.1. The conventional PPE with the homogeneous pressure boundary condition yields the poorest result, as it additionally produces numerical boundary layers adhering to the two side walls. We also find that our proposed boundary condition (10) alone will not generate divergence-free solutions with the conventional PPE (Case 6), nor will our proposed PPE perform well with the conventional boundary condition (Case 2). Because the pressure correction in (7b) and (10) should be consistent between the domain interior and boundaries.

From the time step comparison in Table 1, our proposed PPE and boundary condition (Case 1) demonstrate a significant advantage in stability over the conventional setups (Case 3 and Case 4), albeit Case 4 is also theoretically proven unconditionally stable [6]. For this problem, our new PPE functions adequately with only the linear terms in (7b), and its stability is considerably enhanced when more nonlinear terms are included. Without external ghost points from the continuity, Cases 2 and 3 will have even smaller allowable time steps. The homogeneous pressure boundary condition (Case 5) permits the largest time steps. This is attributed to its numerical simplicity, which enhances the diagonal dominance of the system matrix. From these comparisons, the pressure solutions for Case 2~5 are uniquely determined but flawed. Based on the PDE theory, we conclude that the pressure boundary conditions ((1) and (2)) are enforced excessively, making the system ill-posed instead of being underdetermined as previously hypothesized. The homogeneous pressure boundary condition (Case 5) imposes the most rigid over-constraint and consequently renders the largest errors.

Figure 3. Velocity divergence solutions with the problem setup prescribed in Table 1 for Re = 100 on a 101 x 101 uniform mesh.

Figure 3. Velocity divergence solutions with the problem setup prescribed in Table 1 for Re = 100 on a 101 x 101 uniform mesh.

4.2. Flow Passing a Circular Cylinder

The second validation case focuses on the uniform flow passing a circular cylinder. It is a fundamental external flow problem characterized by flow separation and wake formation. The accurate prediction of these flow features heavily relies on the velocity solution quality and the proper specification of the pressure boundary condition on the cylinder surface. This problem is investigated under steady-state condition at a Reynolds number of 40. As illustrated in Figure 4, the cylinder diameter is set to unity, and the circular computational domain is designed with a radius 25 times the cylinder diameter to mitigate the adverse influence of the finite domain effects on flow predictions. The 541 × 301 structural mesh is clustered around the cylinder and refined in the wake region to better capture the flow details. The near wake region has a radial resolution of 0.01. And all intersecting mesh lines are configured to be orthogonal to one another to reduce discretization errors. The external boundary condition configurations for each variable are also specified in Figure 4. The left semicircle of the domain serves as the inflow boundary with given velocity, static pressure, and Ф. The right semicircle functions as the outflow boundary with the velocity and pressure values obtained by streamwise extrapolation from the domain interior.

As outlined in Table 2, the simulations are carried out using our proposed PPE/pressure boundary condition and the conventional formulations for comparative study. All cases are run with identical time-marching durations to ensure sufficient convergence. The experimental measurements by Coutanceau and Bouard [32] are available as baseline data on the flow separation angle and reattachment distance. Solution flow fields from the three simulations are visualized in Figure 5 with streamlines and velocity divergence contours. Our proposed PPE and boundary condition produce a high-level divergence-free solution (Figure 5(a)). If more iterations are executed, the velocity divergence magnitude can be further reduced. For conventional PPE and boundary condition (Figure 5(b)), the velocity divergence reaches a level of 10^-2. And for conventional PPE with a homogeneous boundary condition (Figure 5(c)), this magnitude increases to the level of 10^-1, with high values concentrated round the cylinder surface.

Table 2 compares the separation angle (θ) and the reattachment distance measured from x = 0 as defined in Figure 5(a) for the three cases. The flow separation point on the cylinder surface and the reattachment point are identified by streamline bifurcations. Our proposed formulation predicts the reattachment distance closest to the experimental measurement, whereas the conventional PPE with homogeneous pressure boundary condition results in the largest deviation. For the separation angle, results of all three cases are very close. They differ from the experimental data by approximately 1 degree, partly attributed to the mesh resolution. This numerical example demonstrates that the correct pressure boundary condition and divergence-free velocity solution are critical in the flow separation prediction, which in turn influences wake length and drag estimates. Our proposed formulations thus provide an important improvement in this predictive capability.

Figure 6 compares the pressure and normal pressure gradient on the upper half-cylinder surface along the central angle, which starts from the rear stagnation point. In Figure 6(a), Case 1 and Case 2 display nearly identical pressure curves, while the curve of Case 3 is slightly smaller in magnitude. Overall, their predictions are in close agreement. That also indicates very similar tangential pressure gradients along the cylinder surface. However, differences between the three cases become more pronounced for the normal pressure gradients, which are post-processed from the pressure solutions using a one-sided scheme. Case 1 resolves a peak normal pressure gradient approximately 10% larger than Case 2, whereas the pressure gradient in Case 3 is trivial due to constraints imposed by its simulation setup. A comparison of Figrures 6(a) and (b) reveals that nearly identical boundary pressure values can correspond to distinct normal gradients, which subsequently impact both pressure and velocity solutions. This observation casts doubt on the validity of the practice adopted by some researchers [10,16], namely, replacing the Neumann pressure boundary condition on solid walls with a Dirichlet boundary condition, even if the Dirichlet-specified pressure values are accurate. Based on the curl-free relationship for the pressure gradients, ∂Pn/∂τ = ∂P_τ/∂n, this distinction in normal pressure gradients also accounts for the difference in tangential pressure gradients next to the wall surface, which subsequently affects the flow separation predictions.