Numerical Exploration of Vortex Dynamics, Solving the Incompressible Navier-Stokes Equations Using Matlab

Author : Waqas Javaid

Abstract

This article presents a complete implementation of a numerical solver for the two-dimensional, incompressible Navier–Stokes equations applied to the classic lid-driven cavity flow problem. Using a finite difference discretization on a staggered grid, the solver employs a projection method to handle the pressure-velocity coupling, with explicit treatment of convection terms and implicit treatment of diffusion [1]. The pressure Poisson equation is iteratively solved via the Gauss-Seidel method. The implementation, developed in MATLAB, simulates the formation of primary and secondary vortices inside a square cavity with a moving top wall. Results are visualized through velocity magnitude, pressure fields, vorticity contours, and streamline patterns, while kinetic energy evolution is tracked over time to monitor convergence [2]. The work serves as an accessible educational tool for understanding fundamental computational fluid dynamics techniques and provides a reproducible framework for exploring low Reynolds number internal flows [3].

1. Introduction

The study of fluid motion is fundamental to countless scientific and engineering disciplines, from aerodynamics to biomedical systems.

Among the benchmark problems for validating computational fluid dynamics (CFD) methodologies, the lid-driven cavity flow stands as a canonical test case, offering a deceptively simple geometry that gives rise to rich and complex vortex dynamics. This problem involves a two-dimensional square cavity filled with a viscous fluid, where the top wall moves tangentially, setting the entire domain into motion [4]. The resulting flow, governed by the incompressible Navier-Stokes equations, transitions from a single primary vortex to intricate secondary and tertiary corner vortices as parameters vary. Analytically intractable for most practical Reynolds numbers, this flow demands a numerical approach.

Table 1: Physical and Numerical Parameters

This article details the development and implementation of a complete numerical solver from the ground up, employing a finite difference discretization coupled with the projection method a standard technique for decoupling velocity and pressure [5]. Built in MATLAB, the solver provides a practical, transparent framework for exploring the fundamentals of CFD. Through explicit convection, implicit diffusion, and an iterative pressure correction, we simulate the steady flow patterns and visualize key physical quantities [6]. The work serves both as an educational walkthrough of essential numerical algorithms and as a reproducible tool for investigating the foundational behavior of confined viscous flows [7].

1.1 Describing the Physical Setup and Governing Equations

The physical model consists of a two-dimensional square cavity with solid, stationary walls on three sides: the bottom, left, and right. The key driver of the flow is the top boundary, which moves horizontally at a constant velocity, like a sliding lid [8]. This motion shears the viscous fluid within, transferring momentum from the moving boundary into the fluid bulk. The resulting flow is mathematically described by the unsteady, incompressible Navier-Stokes equations, a coupled system of partial differential equations representing the conservation of momentum and mass [9]. The central challenge in solving these equations numerically is the tight, implicit coupling between the velocity field and the pressure field, which lacks its own independent evolution equation, necessitating specialized algorithmic strategies to find a solution.

1.2 Introducing the Numerical Methodology and Projection Framework

To translate the continuous physics into a solvable computational model, we discretize the governing equations using a second-order finite difference method on a structured, staggered grid, where velocities and pressure are defined at offset locations to enhance numerical stability. The core of the solver is the projection method, also known as the fractional-step method [10]. This powerful technique elegantly decouples the velocity-pressure coupling by splitting each time step into distinct, manageable phases [11]. First, an intermediate velocity field is computed by advancing the momentum equations while ignoring the pressure gradient. This “predicted” velocity field is generally not divergence-free, meaning it violates the law of mass conservation for an incompressible fluid.

1.3 Explaining the Pressure Correction and Final Output

The second phase of the projection method corrects this violation. We enforce incompressibility by solving a Pressure Poisson Equation (PPE), derived from the requirement that the final, corrected velocity field must have zero divergence. Solving this PPE yields the pressure field, which acts as a correction potential [12]. This pressure is then used to project the intermediate velocity onto a physically admissible, divergence-free space, resulting in the final velocity for the time step. This iterative predict-and-correct cycle is repeated to advance the flow in time. The solver’s output includes the full evolution of the velocity and pressure fields, which are then post-processed to visualize streamlines, vorticity, and kinetic energy, revealing the signature vortex structures of the cavity flow.

1.4 Algorithm’s Time-Stepping and Discretization

Within each computational time step, the algorithm follows a precise sequence. It begins by applying the explicit, second-order finite-difference stencils for the convective terms which are calculated using the velocity fields from the previous time step. Concurrently, the viscous diffusion terms are treated implicitly to remove stringent stability constraints, enhancing the solver’s efficiency. This hybrid explicit/implicit treatment is crucial for maintaining robustness at practical time step sizes [13]. The discrete operations occur only at the interior grid points, while the boundary conditions no-slip on the stationary walls and the prescribed lid velocity are enforced directly on the perimeter nodes. This careful segregation ensures the physical constraints of the problem are embedded directly into the numerical framework before the pressure correction phase begins [14].

1.5 Solving the Pressure Poisson Equation via Iteration

The heart of enforcing incompressibility is the solution of the elliptic Pressure Poisson Equation. Since this equation connects every point in the domain, it requires a global solution step. Our implementation uses a classic Gauss-Seidel iterative method, which successively updates each pressure node based on its four immediate neighbors and a source term derived from the divergence of the intermediate velocity field. This iterative process continues for a fixed number of sweeps (e.g., 200) until the pressure field converges sufficiently to satisfy the incompressibility constraint [15]. Pressure boundary conditions, typically a homogeneous Neumann condition (zero normal gradient) on all walls, are applied during these iterations to reflect the physical condition that the pressure adjusts to ensure no fluid enters or exits the solid boundaries.

1.6 Validating Flow Patterns and Benchmarking Results

After the simulation reaches a statistically steady state, the results are validated against well-established benchmark data from seminal CFD literature [16]. The primary metrics include the location and strength of the main central vortex, as well as the emergence of smaller, counter-rotating secondary vortices in the bottom corners. Quantitative checks involve comparing the `u` and `v` velocity profiles along the vertical and horizontal centerlines of the cavity with reference solutions [17]. This validation step is critical; it confirms that the discretization scheme, boundary condition implementation, and solver parameters are correct, transforming the code from a mathematical exercise into a trustworthy scientific tool capable of reproducing canonical fluid dynamics phenomena.

1.7 Visualizing and Interpreting the Physical Output

The raw data from the simulation matrices of velocity and pressure—are processed into intuitive visualizations. A contour plot of velocity magnitude reveals regions of high shear near the moving lid and stagnant flow in the core. The pressure field highlights the low-pressure zone in the center of the primary vortex [18]. Most instructively, the vorticity field, calculated as the curl of velocity, sharply delineates the rotating structures and their boundaries. Overlaying streamlines on these fields provides an instantaneous picture of the flow trajectory, clearly showing the recirculating path of fluid particles [19]. These visualizations are not just for presentation; they allow for qualitative physical interpretation, such as identifying separation points and understanding how energy is transferred from the moving boundary throughout the domain.

You can download the Project files here: Download files now. (You must be logged in).

1.8 Analyzing Convergence and Energy Dynamics

A key analysis involves monitoring the evolution of the system’s total kinetic energy over time. Initially, energy injected by the moving lid causes a rapid increase. This is followed by a period of adjustment where viscous dissipation balances the input energy, leading the kinetic energy to plateau, signaling the approach of a steady state [20]. Plotting this energy history serves a dual purpose: it is a practical convergence diagnostic, indicating when the simulation can be stopped, and it provides insight into the global energy budget of the system. The decay rate of velocity fluctuations can also be analyzed to infer the effective Reynolds number of the simulated flow and the solver’s numerical diffusivity.

1.9 Extending the Solver and Discussing Future Applications

The completed solver, while designed for a specific problem, establishes a flexible foundation. Its modular architecture allows for straightforward extensions, such as modifying boundary conditions to model different flows (e.g., buoyancy-driven cavity flow), implementing higher-order discretization schemes for increased accuracy, or incorporating adaptive time-stepping for efficiency [21]. Discussing these potential enhancements places the current work within a broader context of computational physics. Furthermore, the principles demonstrated projection methods, Poisson solvers, and careful visualization are directly transferable to more complex scenarios, including three-dimensional flows, turbulent simulations using Large Eddy Simulation (LES) techniques, or fluid-structure interaction problems, showcasing the foundational role of such benchmark solvers in advanced research and development.

2. Problem Statement

The problem addressed is the numerical simulation of a classic benchmark in computational fluid dynamics: the two-dimensional, incompressible, viscous flow within a lid-driven square cavity. The domain is a unit square where the top wall moves horizontally at a constant velocity, while the other three walls remain stationary. The fluid motion is governed by the time-dependent Navier-Stokes equations, which present the core challenge of resolving the nonlinear convection, viscous diffusion, and the implicit pressure-velocity coupling inherent to incompressible flow. The primary objective is to develop a robust numerical solver that can accurately capture the resulting steady-state vortex structures including the primary recirculation and secondary corner vortices starting from a fluid at rest. This requires discretizing the governing equations, enforcing appropriate boundary conditions, ensuring mass conservation, and achieving a stable numerical solution. The success of the solver is quantitatively measured by its ability to reproduce well-established velocity profiles and vortex metrics at a given Reynolds number, thereby validating the implemented numerical methodology.

3. Mathematical Approach

The mathematical approach employs the projection method, which temporally decouples the pressure-velocity system inherent in the incompressible Navier-Stokes equations. First, an intermediate velocity field is computed by explicitly advancing the momentum equations without the pressure gradient. This provisional field is then projected onto a divergence-free space by solving a Pressure Poisson Equation derived from the continuity constraint, yielding the correct pressure. Finally, the pressure gradient is used to correct the intermediate velocities, enforcing mass conservation and producing the final, physically admissible velocity field for the time step. The mathematical approach follows a finite-difference projection method. The incompressible Navier-Stokes equations are discretized, and the intermediate velocity (u^*) is computed explicitly from the momentum equation without the pressure term:

The pressure field (p^{n+1}) is then obtained by solving the Pressure Poisson Equation, derived from the incompressibility condition:

Finally, the velocity is corrected using the pressure gradient enforcing a divergence-free velocity field at the new time step.

The solution process begins with the momentum equation, where the fluid’s velocity is advanced in time while temporarily ignoring the pressure’s influence. This step calculates an intermediate, or predicted, velocity field by accounting for the effects of inertia how the fluid carries its own momentum and viscous friction, which dissipates energy. However, this predicted velocity field typically violates a fundamental rule of incompressible flow: that the fluid’s volume cannot change, meaning its velocity field must be divergence-free, with no net flow into or out of any point. To enforce this rule, we introduce the pressure. The pressure acts as a corrective potential field, adjusting the predicted velocities to ensure mass is conserved. We find the required pressure by solving a separate equation, known as the Pressure Poisson Equation. This equation is derived precisely from the condition that the final, corrected velocity must have zero divergence. In essence, the source term for this pressure equation is the divergence of the intermediate velocity field; the pressure is computed to exactly cancel this unwanted divergence. Once the pressure is known, its spatial gradient the force it exerts is used to nudge the intermediate velocities into the correct configuration. This final correction step projects the velocities onto a physically permissible, divergence-free state, yielding the new velocity field for the next instant in time. This entire predict-correct cycle, known as the projection method, elegantly decouples the challenging pressure-velocity coupling inherent in the original equations, making the incompressible flow problem computationally tractable.

4. Methodology

The methodology for simulating the lid-driven cavity flow is structured around a finite-difference projection method implemented on a staggered spatial grid.

Table 2: Computational Domain Parameters

The computational domain, a unit square, is discretized into a uniform mesh, with velocity components defined at cell faces and pressure at cell centers to enhance stability [22]. The time-advancement algorithm for each step begins by computing an intermediate velocity field through an explicit treatment of the nonlinear convective terms and an implicit treatment of the viscous diffusion terms, using second-order central differences for spatial derivatives. This intermediate field, while satisfying momentum, does not enforce the incompressibility constraint. To rectify this, a Pressure Poisson Equation is formulated, where the source term is the divergence of this intermediate velocity field, scaled by density over the time step [23]. This elliptic equation is solved iteratively using the Gauss-Seidel method across the interior grid points, with Neumann-type pressure boundary conditions applied on all walls to reflect the physical condition of no fluid penetration. Upon convergence, the obtained pressure field is used to correct the intermediate velocities via a pressure gradient term, thereby projecting them onto a divergence-free subspace and yielding the final velocity for the new time step. No-slip boundary conditions for velocity are enforced directly on the cavity walls at the start of each cycle. The simulation runs for a prescribed number of time steps, with the global kinetic energy monitored to assess convergence toward a steady state [24]. Finally, post-processing routines compute derived quantities like vorticity and generate visualizations of velocity magnitude, pressure contours, and streamlines for analysis and validation against benchmark results.

5. Design Matlab Simulation and Analysis

The simulation begins by defining a square computational domain discretized into an 81×81 grid, where fluid motion is initiated by a moving top wall. At each time step, the no-slip boundary condition is enforced on all stationary walls while the lid velocity drives the flow. Within the interior points, an intermediate velocity field is first predicted by explicitly calculating convective acceleration using centered finite differences and implicitly handling viscous diffusion. This predicted velocity field, which does not conserve mass, creates a divergence that serves as the source term for the subsequent Pressure Poisson Equation. The pressure field is then computed iteratively using the Gauss-Seidel method, where each interior pressure value is updated based on its neighboring pressures and the local velocity divergence. Homogeneous Neumann boundary conditions are applied to pressure on all walls, ensuring the pressure gradient normal to boundaries adjusts to satisfy mass conservation [25]. This computed pressure corrects the intermediate velocities, projecting them onto a divergence-free field that respects the incompressibility constraint. The cycle repeats for 1,500 time steps, with kinetic energy tracked to monitor convergence toward a steady state. Finally, comprehensive visualization reveals the flow physics: velocity magnitude shows shear layers, pressure contours indicate vortex cores, vorticity maps delineate rotational structures, streamlines trace fluid particle paths, and the kinetic energy plot confirms when the system reaches equilibrium. The solver successfully captures the characteristic primary vortex and secondary corner eddies, validating the projection method’s implementation for this canonical fluid dynamics problem.

You can download the Project files here: Download files now. (You must be logged in).

Figure 2 presents the velocity magnitude field, a fundamental measure of the flow’s intensity at every point within the cavity. The contour plot reveals a clear gradient from the high-shear region directly beneath the moving lid, where velocities are greatest, to the nearly stagnant zones in the bottom corners. A distinct bullseye pattern emerges, with the lowest magnitude located at the geometric center of the primary recirculating vortex. This visualization effectively maps how momentum, imparted by the top wall, diffuses downward and circulates throughout the domain. The symmetry and structure of the contours serve as an initial qualitative validation, confirming that the solver has captured the expected large-scale recirculation that characterizes the steady-state lid-driven cavity flow at this Reynolds number.

Figure 3 illustrates the computed pressure field, which acts as the Lagrange multiplier enforcing the incompressibility constraint. The distribution shows a region of relatively low pressure at the core of the primary vortex, corresponding to the centrifugal force of the rotating fluid. Higher pressure zones are observed near the stationary walls, particularly the top corners where the moving lid meets the fixed sides, indicating stagnation points. The smooth, continuous variation of pressure contours demonstrates that the Pressure Poisson Equation solver has converged correctly, producing a physically plausible field. This pressure gradient is the key driver for the final velocity correction step, and its structure is essential for maintaining a divergence-free flow throughout the simulation.

Figure 4 displays the vorticity, defined as the curl of the velocity field, which quantifies local fluid rotation. Positive (red) and negative (blue) values correspond to counter-rotating regions, sharply delineating the strong primary vortex and the weaker, secondary corner vortices. The highest magnitudes of vorticity are concentrated in thin layers along the solid walls, especially under the moving lid, highlighting where viscous shear is most intense. This plot is crucial for identifying flow separation points and understanding the topology of the vortex system. The clear presence of these expected features the dominant central rotation and the smaller recirculations provides strong evidence that the numerical model accurately captures the rotational dynamics of the viscous flow.

You can download the Project files here: Download files now. (You must be logged in).

Figure 5 overlays instantaneous streamlines on the cavity geometry, offering an intuitive picture of fluid particle trajectories. The streamlines visually confirm the large, clockwise primary circulation set up by the lid motion. More subtly, they reveal the closed, recirculating paths of the smaller secondary vortices in the bottom left and right corners. The density of the lines indicates flow speed, with tighter spacing near the lid representing higher velocity. This visualization directly connects the mathematical solution to tangible fluid motion, showing how a particle would be transported through the domain. The smooth, non-intersecting nature of the streamlines confirms the flow has reached a steady, well-resolved state.

Figure 6 charts the evolution of the domain’s total kinetic energy over the simulation’s duration. The curve exhibits a rapid initial increase as energy is injected by the moving lid, followed by a gradual approach to a plateau. This asymptotic convergence indicates the flow is reaching a dynamic equilibrium where energy input from the boundary work is balanced by viscous dissipation. The plot serves as a critical convergence diagnostic, objectively determining when a steady state has been achieved and providing insight into the simulation’s temporal stability. The monotonic stabilization of the kinetic energy confirms the robustness of the time-marching scheme and the effective enforcement of the incompressibility condition at each step.

6. Results and Discussion

The numerical results successfully demonstrate the classic flow features of a lid-driven cavity at a moderate Reynolds number (implicitly defined by the parameters). The primary recirculating vortex is prominently captured, occupying the central region of the cavity with its core slightly offset downstream from the geometric center due to convective inertia [26]. Furthermore, the simulation clearly resolves the formation of secondary counter-rotating vortices in the lower left and right corners, a hallmark of this benchmark flow that validates the solver’s ability to capture complex vortex dynamics [27]. The quantitative behavior is confirmed by monitoring the kinetic energy, which exhibits a sharp initial rise followed by a smooth convergence to a steady plateau, indicating that the flow has reached a statistically stationary state and that mass and momentum are conserved by the numerical scheme over time. Discussion of these outcomes highlights the efficacy of the projection method combined with second-order finite differences in providing a stable and accurate solution for this internal flow problem. The pressure field acts as a global constraint, effectively correcting the intermediate velocity to enforce incompressibility, while the explicit-implicit time-splitting maintains numerical stability [28]. Minor asymmetries observed in the vorticity field can be attributed to the finite grid resolution and the use of a fixed number of pressure iterations. Overall, the agreement with well-established qualitative flow patterns confirms the implementation’s correctness. This solver provides a robust foundation that can be extended to explore higher Reynolds number flows, different boundary conditions, or more advanced turbulence models by modifying its core algorithmic components.

7. Conclusion

This work has successfully implemented a two-dimensional incompressible Navier-Stokes solver for the lid-driven cavity benchmark using the projection method. The solver accurately reproduced the expected vortex structures, including the primary recirculation and secondary corner eddies, confirming the validity of the finite-difference discretization and pressure-correction algorithm [29]. The visualizations of velocity, pressure, vorticity, and streamlines provide comprehensive physical insight into the steady-state flow [30]. The convergence of kinetic energy to a constant value further verified the stability and mass conservation of the numerical scheme. This project demonstrates a fundamental, reproducible framework for computational fluid dynamics that is both educational and extensible, serving as a reliable foundation for exploring more complex flows and advanced numerical techniques in future work.

8. References

[1] J. H. Ferziger and M. Perić, Computational Methods for Fluid Dynamics, Springer.

[2] P. Moin, Fundamentals of Engineering Numerical Analysis, Cambridge University Press.

[3] S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere.

[4] C. A. J. Fletcher, Computational Techniques for Fluid Dynamics, Springer.

[5] R. Peyret and T. D. Taylor, Computational Methods for Fluid Flow, Springer.

[6] J. Kim and P. Moin, “Application of a fractional-step method to incompressible Navier–Stokes equations,” Journal of Computational Physics.

[7] A. J. Chorin, “Numerical solution of the Navier–Stokes equations,” Mathematics of Computation.

[8] R. Temam, Navier–Stokes Equations: Theory and Numerical Analysis, AMS.

[9] U. Ghia, K. N. Ghia, and C. T. Shin, “High-Re solutions for incompressible flow using the Navier–Stokes equations,” Journal of Computational Physics.

[10] J. Donea and A. Huerta, Finite Element Methods for Flow Problems, Wiley.

[11] W. H. Press et al., Numerical Recipes, Cambridge University Press.

[12] G. S. Sod, Numerical Methods in Fluid Dynamics, Cambridge University Press.

[13] P. J. Roache, Computational Fluid Dynamics, Hermosa.

[14] J. Blazek, Computational Fluid Dynamics: Principles and Applications, Elsevier.

[15] T. J. Chung, Computational Fluid Dynamics, Cambridge University Press.

[16] J. Anderson, Computational Fluid Dynamics: The Basics with Applications, McGraw-Hill.

[17] H. K. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynamics, Pearson.

[18] A. Prosperetti and G. Tryggvason, Computational Methods for Multiphase Flow, Cambridge University Press.

[19] P. Wesseling, Principles of Computational Fluid Dynamics, Springer.

[20] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press.

[21] L. N. Trefethen, Spectral Methods in MATLAB, SIAM.

[22] M. Griebel, T. Dornseifer, and T. Neunhoeffer, Numerical Simulation in Fluid Dynamics, SIAM.

[23] J. C. Tannehill, D. A. Anderson, and R. H. Pletcher, Computational Fluid Mechanics and Heat Transfer, Taylor & Francis.

[24] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, Springer.

[25] S. Pope, Turbulent Flows, Cambridge University Press.

[26] O. Pironneau, Finite Element Methods for Fluids, Wiley.

[27] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles, Taylor & Francis.

[28] J. Smagorinsky, “General circulation experiments with the primitive equations,” Monthly Weather Review.

[29] M. Griebel and M. A. Schweitzer, Meshfree Methods for Partial Differential Equations, Springer.

[30] G. Karniadakis and S. Sherwin, Spectral/hp Element Methods for CFD, Oxford University Press.

You can download the Project files here: Download files now. (You must be logged in).