A pseudo-transient pressure gradient method for solving the incompressible Navier–Stokes equations

Numerical solutions of incompressible viscous flow equations encounter a fundamental problem due to the lack of an evolutionary equation for the pressure field. That is, from a given pressure field, one can easily calculate the corresponding velocity field using the momentum equations; however, there is no explicit relation that can be utilized to update the pressure field from a known velocity field.¹ Numerous methods have been devised to resolve this issue. Certain methods modify the continuity equation to include a pressure related term, while others mathematically derive a pressure equation that complies to the divergence-free velocity field constraint. Some other methods include eliminating the pressure term from the governing equations. These methods can be categorized into two groups: (1) primitive variable formulation methods in which the governing equations are written in terms of the fundamental dependent variables (velocity components and pressure), and (2) non-primitive variable formulation methods in which the fundamental dependent variables in the governing equations are replaced with other non-primitive variables, such as vorticity, stream function, and potential function. Despite the acceptable results of these methods, they have some inherent limitations that can pose major difficulties in certain types of problems in computational fluid dynamics (CFD).

The non-primitive variables formulation relies on eliminating the pressure term from the governing equations by replacing the fundamental fluid variables, velocity components and pressure, with derived fluid variables such as the vorticity, stream function, and potential function.² In other words, it decouples the kinematic variables from the kinetic variables.³ One of the prominent examples of this formulation is the vorticity–stream function formulation.^4–6 In this approach, one takes the curl of the momentum equation to generate a second-order partial differential equation called the vorticity transport equation and utilizes the stream function in such a way that it satisfies the divergence-free velocity constraint for incompressible flows $∇·V→=0$⁠. This method is only applicable to two-dimensional problems as it is the only case where the stream function can analytically satisfy the continuity equation.² An extension of this approach to three-dimensional problems was introduced by Azziz and Hellums.⁷ Among the various variations of this approach, it is worth mentioning the dual-potential function formulation. This formulation exploits the Helmholtz decomposition theorem, in which a vector potential function consisting of a divergence-free part (rotational) and a curl-free part (irrotational) is used to generate the desired velocity field.^8,9

Gautam and Albdallah¹⁰ utilized the concept of the Helmholtz decomposition theorem to resolve the cell Reynolds number stability limitation encountered in the convection–diffusion Navier–Stokes equations. After substituting the proposed potential function into the transient term of the momentum equation, the potential component appears as a pressure-like term, hence, combines with the pressure term of the momentum equation. Also, the rotational component appears as a vortex-like term; therefore, it combines with the viscous term of the momentum equation. This formulation eliminates the convection term that invokes the cell Reynolds number limitation and transforms the unsteady momentum equation to an explicit steady-state equation.

Another formulation similar to the vector potential formulation is the vorticity–velocity formulation. This approach eliminates the pressure term from the equation and maintains the velocity components, which results in six partial differential equations to be solved simultaneously: three parabolic equations for vorticity and three Poisson-type equations for velocity.^11–14 It is worth mentioning that all non-primitive formulations solve for the velocity field; consequently, to get the pressure field, one needs to solve the pressure Poisson equation that is obtained by taking the divergence of the momentum equation.^15,16 Moreover, in these formulations, the vorticity boundary conditions rely on the velocity values within the flow domain; therefore, the vorticity values are always lagging by one iteration.

In primitive variable formulations, the governing equations are written in terms of the primary flow variables, velocity components and pressure; hence, these methods are characterized by their intuitive representation of the physical phenomena associated with the fluid flow. They can be classified under two distinct categories: (1) pseudo-compressibility methods introduced by Chorin¹⁷ in 1967 and (2) pressure Poisson equation methods introduced by Harlow and Welch¹⁸ in 1965. First, the core principle of the pseudo-compressibility methods lies in modifying the continuity equation in such a way that it includes a pressure term that gradually vanishes as solution converges. For instance, Chorin¹⁷ introduced a time-derivative pressure term to the continuity equation for steady flow problems; nonetheless, researchers have utilized it for many unsteady flow problems as well.^19,20 This technique, however, tends to exhibit spurious pressure oscillations when it is utilized on collocated grids.^21–23 Many approaches have been proposed to resolve this issue; for example, one-sided discretization, flux splitting,²⁴ regularized central differencing,²⁵ and consistent flux reconstruction.²⁶ 

In the pressure Poisson methods, the continuity equation is replaced by a second-order elliptic non-homogeneous partial differential equation derived by taking the divergence of the momentum equation. In this method, the pressure Poisson equation is used to compute the pressure field whereas the velocity field is computed from the momentum equations. Typically, the boundary conditions associated with this pressure equation are of the Neumann type and obtained by rewriting the momentum equation at the boundary of the flow domain; thus, the solution of this equation exists only when the boundary conditions and the source term are correlated by a compatibility condition. This condition affirms that the numerical integration of the source term over the flow domain is equal to the numerical integration of the derivative boundary conditions (Green's theorem). S. Abdallah¹⁶ introduced a novel method to accommodate this condition on a collocated grid. Alternatively, to eliminate the need for this compatibility constraint, one could use the so-called pressure gradient method in which the Neumann pressure boundary conditions are replaced by Dirichlet pressure gradient boundary conditions. The notion of using the pressure gradient as a dependent variable was first introduced by Shih²⁷ and later independently by Said.²⁸ It was found that at high Reynolds numbers, oscillatory velocity and pressure fields are produced. Yoshihiro²⁹ modified the pressure gradient method by utilizing a weighted average pressure gradient term in order to damp these oscillations. Later, Marshaal and Shaaban³⁰ extended the modified version to the generalized curvilinear coordinates system.

In summary, methods for solving incompressible flow equations can be categorized into two groups: the non-primitive variables group, which often faces challenges in deriving the corresponding boundary conditions and is typically limited to two-dimensional problems, and the primitive variables group, which is applicable to both two-dimensional and three-dimensional problems but sometimes fails to properly satisfy the continuity constraint. This study introduces a new primitive variable approach that ensures the continuity constraint is met to machine precision at every grid point within the flow domain. The principle of this method lies in treating the pressure gradient as a stand-alone variable where a pseudo-transient pressure gradient (PTPG) term is added to the gradient of the continuity equation. Hence, we call this the PTPG method. Upon convergence, the continuity equation is satisfied to a constant value enforced at some point within the flow domain. In Sec. II, we show the mathematical formulation of the new method and the finite difference discretization of the new system of equations. In Sec. III, we validate the presented method by solving three steady-state benchmark problems and comparing the results with data available in the literature.

Incompressible flow in a square cavity is a canonical flow problem that has been extensively used for validating numerous incompressible flow solvers. This problem consists of a unit cavity with no-slip and no-flux boundary conditions on all walls, except for the top wall, which represents a moving lid with a tangential unit velocity. The prominence of this problem is attributed to its simple geometry and its complex flow structure, which consists of a large primary vortex near the center of the cavity with two small rotating vortices at the bottom corners. Furthermore, numerous reliable numerical solutions are available in the literature, which provide an extensive source of validation data. Numerical solutions for Re = 100, 400, and 1000 are obtained in this work. A uniform grid with 201 × 201 elements is used for Re = 100 and 400, and 301 × 301 elements for Re = 1000, Table II shows the grid sensitivity analysis for each Reynolds number.

Figure 4 presents the convergence history of dilation (divergence of the velocity field), velocity components, and pressure gradients for the case of Re = 1000, which represents the flow physics at moderate Re. The $L∞$ error, computed as the difference between the current value and the previous value of a variable, is used for residual representation. It is evident that the presented dependent variables successfully attain different levels of convergence, on the order of $10−16$ for velocity components, $10−14$ for dilation, and $10−13$ for pressure gradients. Figures 5(a) and 5(b) compare the velocity divergence contours obtained using the presented method and the SIMPLE scheme in ANSYS Fluent, respectively. It has been noted that the divergence of the velocity field reaches very small levels over the entire domain, less than 10⁻¹³, when the new method is used. On the other hand, when the SIMPLE scheme is employed, the divergence of the velocity field is not satisfied in the cells near the top corners, whereas in the rest of the flow domain it ranges from $10−3$ to $10−9$⁠. Thus, the new numerical method can provide well satisfied continuity constraint over the entire domain.

The velocity profiles at the horizontal and vertical centerlines are obtained using the proposed method and compared with those obtained using the stream function–vorticity method by Ghia et al.³¹ in Fig. 6. The results demonstrate excellent agreement between the two methods. For instance, at Re = 1000, the inflection point of the horizontal velocity profile obtained by the stream function–vorticity approach is at $y=0.171$⁠, with a corresponding velocity $u=−0.382$⁠. The new method depicts a very close result, where the inflection point is at $y=0.172$⁠, with a corresponding velocity $u=−0.386$⁠. Similar behavior is observed in the vertical velocity profile where the two peaks are attained at the same location by the two methods. Moreover, the two profiles match the physical interpretation of the problem. For example, the horizontal velocity component shown in Fig. 6(a) comprises two parts, a positive part and a negative part. Integrating this profile along the $y$-axis yields a value of $3.5×10−5$⁠, illustrating the principle of conservation of mass. Similarly, integrating the vertical velocity profile, shown in Fig. 6(b), along the $x$-axis yields a value of − $3.82×10−6$ depicting the same conservation principle. Finally, the horizontal and vertical velocity contours align perfectly with the velocity contours found in the literature,¹⁰  Fig. 7.

To evaluate the applicability of the new method for three-dimensional fluid flow problems, the three-dimensional driven cavity flow is analyzed at a Reynolds number (Re) of 100. The centerline horizontal velocity within the central plane exhibited a trend similar to that obtained in the case of two-dimensional driven cavity flow, demonstrating the method's consistency and reliability across different dimensions, Fig. 8.

For many CFD problems, the use of non-orthogonal curvilinear grid is inevitable. To examine the utility of the presented method on such grids, the flow inside a 2D constricted channel with Reynolds number equal to 100 is solved. The employed computational domain along with the corresponding boundary conditions is illustrated in Fig. 9. The length of the channel is 6 times its height, and the bump height is 0.1 on each side. Mesh independence is attained with $201×1201$ grid points. The convergence history of all dependent variables and the contour of velocity divergence are shown in Fig. 10. Evidently, all dependent variables converge well, and the velocity divergence inside the computational domain is less than $10−13$⁠. A comparison case has been configured using ANSYS Fluent for the same Reynolds number. The continuity contour obtained using the SIMPLE scheme shows very high values, on the order of $10−2$⁠, compared to that obtained using the new method, Fig. 11. The obtained velocity profiles exhibit excellent agreement with those calculated using the new method, Fig. 12. Figure 13 shows the velocity contour inside the domain. This contour aligns perfectly with the intuitive physical understanding of such problem. This is to say for a constant flow rate, the fluid velocity increases when it passes through a constricted area, conservation of mass principle.

A new method called the pseudo-transient pressure gradient (PTPG) method is developed based on the spirit of the artificial compressibility method and presented as a new approach to solve the incompressible Navier–Stokes equations. The principle of the described method lies in its consideration of the pressure gradient as a pseudo-transient term added to the gradient of the continuity equation, which creates an additional set of differential equations that must be solved with the momentum equations. Upon convergence, the pseudo-transient term vanishes, and the continuity constraint is satisfied to machine zero. This numerical method eliminates the compatibility and boundary condition problems associated with primitive variable approaches and improves the convergence of all dependent variables. The finite difference discretization scheme is utilized on collocated grids to solve three steady-state validation problems: lid-driven cavity flow for Re = 100, 400, and 1000, fully developed channel flow for Re = 100, and flow in a constricted channel for Re = 100. The obtained solutions are free of oscillations and compare well with the benchmark results available in the literature. Furthermore, the divergence of the velocity field approaches machine zero in all examined problems. Therefore, the PTPG method is a robust primitive variables method that can be used to solve both the steady and the unsteady incompressible Navier–Stokes equations, and it is simple to apply in Cartesian and curvilinear coordinates systems, as well as in 2D and 3D problems.

The authors have no conflicts to disclose.

Ahmad M. Alsaghir: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Shashank Mishra: Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal). Shaaban Abdallah: Conceptualization (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). Je-Hyeong Bahk: Conceptualization (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.