High-Resolution CFD Analysis of Incompressible Flow in a Lid-Driven Cavity Using Explicit RK2 and SOR Techniques

Author : Waqas Javaid

Abstract:

This study presents a detailed computational investigation of two-dimensional incompressible flow inside a square lid-driven cavity using the streamfunction–vorticity formulation of the Navier–Stokes equations. The numerical model is developed entirely in MATLAB to provide a flexible and efficient research framework for high-resolution fluid dynamic simulations. A second-order explicit Runge–Kutta (RK2) method is employed for time integration, while the Poisson equation for the streamfunction is solved iteratively using the Successive Over-Relaxation (SOR) technique. Computational Fluid Dynamics (CFD) provides a numerical framework to analyze fluid flow behavior under various physical conditions [1]. The top wall moves at a constant velocity, generating vortical structures that evolve into steady-state recirculating flow patterns. The solver accurately captures the formation of the primary vortex at the cavity center and secondary eddies near the corners. Spatial discretization employs a uniform structured grid with central differencing for diffusion and convection terms, ensuring numerical stability and accuracy. The governing principles of viscous, incompressible flow are derived from the Navier–Stokes equations, which describe momentum conservation in fluids [2]. Simulations are performed for Reynolds numbers up to 1000, highlighting the transition from laminar to more complex flow behavior. Separate output visualizations, including streamfunction contours, vorticity fields, velocity vectors, and 3D surface plots, demonstrate the physical validity of the computed results. The proposed model achieves convergence with minimal numerical oscillations and demonstrates robustness for long-term time integration. Comparison with benchmark data confirms the reliability of the present approach for incompressible flow problems. This research provides a reusable MATLAB platform for advanced computational fluid dynamics studies and forms a strong basis for extending to turbulent and three-dimensional flow analysis. The finite volume and finite difference approaches have been widely used to discretize and solve flow equations efficiently [3].

1.    Introduction:

 The study of incompressible fluid motion within confined geometries has long served as a benchmark problem in Computational Fluid Dynamics (CFD), providing valuable insights into flow structure, vorticity transport, and numerical stability. Among such problems, the two-dimensional lid-driven cavity flow remains one of the most extensively investigated cases due to its simple geometry yet complex internal flow behavior. It exhibits a rich combination of primary and secondary vortex formations, shear-layer interactions, and recirculating flow zones that are ideal for validating numerical solvers. Kim and Moin introduced the fractional-step method for time advancement, which significantly improved solution accuracy for incompressible flows [4]. Patankar’s work on numerical heat transfer established a foundation for pressure–velocity coupling in CFD solvers [5]. The present research focuses on the numerical simulation of the lid-driven cavity flow using the streamfunction–vorticity formulation of the Navier–Stokes equations implemented in MATLAB.

This approach eliminates pressure as a primary variable, simplifying the solution procedure while maintaining high accuracy for incompressible flow conditions. The governing equations are discretized using a finite-difference scheme on a uniform structured grid, ensuring second-order spatial accuracy. Peyret and Taylor provided detailed formulations for the numerical treatment of viscous flows and benchmark cavity problems [6]. Time integration is performed using the explicit second-order Runge–Kutta (RK2) method, while the Poisson equation for the streamfunction is efficiently solved using the Successive Over-Relaxation (SOR) technique. The numerical code is designed to be modular, computationally stable, and easily extendable for higher Reynolds numbers and transient flow studies. The cavity top wall is driven at a constant velocity, generating circulation patterns that evolve toward a steady-state solution. The simulation captures the essential physics of viscous shear-driven motion and demonstrates excellent convergence characteristics. Detailed post-processing, including contour plots, velocity profiles, and three-dimensional visualizations, provides deeper physical interpretation of the flow field

This research aims to establish a reliable MATLAB-based CFD framework that can serve both as an educational tool and as a foundation for advanced doctoral-level investigations in fluid mechanics. . Hirsch emphasized that the accuracy of CFD simulations depends strongly on spatial resolution and boundary condition implementation [7].

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

1.1    Background and Motivation:

The lid-driven cavity problem is one of the most fundamental and extensively studied benchmarks in Computational Fluid Dynamics (CFD). It provides valuable insights into the behavior of incompressible viscous flows within a simple yet physically rich domain. Ghia et al  [8] produced benchmark results for the lid-driven cavity flow, which are widely used for validating numerical solvers. The configuration consists of a square cavity where the top lid moves tangentially at a constant velocity, while the remaining walls remain stationary. Despite its geometric simplicity, this problem exhibits complex flow patterns such as recirculating vortices, shear layers, and corner eddies, making it an ideal case for validating numerical schemes and studying internal flow dynamics.

1.2    Physical and Mathematical Description:

The flow inside the cavity is governed by the two-dimensional incompressible Navier–Stokes equations, which represent the conservation of mass and momentum. To simplify the problem and eliminate the pressure term, the equations are reformulated in the streamfunction–vorticity form. This transformation allows the flow field to be expressed using only two scalar variables: streamfunction (psi) and vorticity (omega). The lid-driven cavity remains an ideal test case for assessing numerical methods in laminar and transitional regimes [9]. The approach efficiently handles the boundary conditions and provides stable numerical convergence for incompressible flows.

1.3    Numerical Methodology:

In this study, a finite-difference-based numerical scheme is implemented using MATLAB. The computational domain is discretized with a uniform structured grid to achieve second-order spatial accuracy. Time integration of the vorticity transport equation is performed using the second-order Runge–Kutta (RK2) method, which offers a balance between accuracy and computational cost. The streamfunction equation, which is a Poisson equation, is solved iteratively using the Successive Over-Relaxation (SOR) technique to ensure faster convergence to steady-state. Harlow and Welch developed one of the earliest finite-difference models for viscous incompressible flow using a staggered grid arrangement [10].

1.4    Boundary Conditions and Flow Parameters:

The boundary conditions applied in this simulation correspond to the physical behavior of the lid-driven cavity. The top boundary (lid) moves with a constant velocity, generating the main circulation, while the bottom and side walls remain stationary and enforce no-slip conditions. The simulation is performed for a Reynolds number of ( Re = 1000 ), representing a steady laminar regime with strong recirculations and distinct secondary vortices.

 1.5    Research Significance and Objectives:

The main objective of this study is to develop a robust and extensible MATLAB-based CFD solver capable of simulating the lid-driven cavity flow using the streamfunction–vorticity formulation. This solver not only serves as a reliable research tool but also as an educational platform for understanding the fundamentals of CFD. The Runge–Kutta method used in the present study ensures second-order temporal accuracy and numerical stability [11]. The SOR (Successive Over-Relaxation) algorithm is applied to solve the Poisson equation efficiently for streamfunction computation [12]. The presented results, including six detailed visual outputs streamfunction, vorticity, velocity vectors, centerline velocity profiles, and 3D surface plots provide a comprehensive understanding of the cavity flow characteristics. The study ultimately demonstrates MATLAB’s capability as a powerful environment for implementing, testing, and visualizing advanced CFD algorithms.

2.    Problem Statement:

The lid-driven cavity problem represents a classical benchmark in Computational Fluid Dynamics (CFD) for analyzing two-dimensional incompressible viscous flow in a confined domain. Despite its geometric simplicity, the flow exhibits complex patterns, including primary and secondary vortices, which make it a challenging numerical test case. The main difficulty lies in accurately resolving the velocity field, pressure gradients, and vorticity transport within the cavity, especially at higher Reynolds numbers. Traditional finite-volume solvers often require significant computational resources, motivating the development of efficient alternatives. This study addresses the problem by implementing a streamfunction–vorticity formulation using MATLAB which removes the pressure term and simplifies computation. The goal is to achieve numerical stability, accurate vortex resolution, and smooth convergence toward steady-state flow. The simulation focuses on the lid-driven cavity at a Reynolds number of 1000, employing the Runge–Kutta (RK2) time integration and the Successive Over-Relaxation (SOR) method for the Poisson equation. The final objective is to provide a reliable and extensible CFD solver suitable for both academic research and educational purposes.

3.    Mathematical approach:

The present study employs the two-dimensional incompressible Navier–Stokes equations in streamfunction–vorticity formulation to simulate the lid-driven cavity flow. Thom’s vorticity boundary formulation provides an accurate approximation for solid wall boundaries in two-dimensional flows [13]. The governing equations are derived from the fundamental laws of fluid mechanics, namely, the continuity and momentum conservation equations. The continuity equation for an incompressible fluid is expressed as:

 ∂u/∂x + ∂v/∂y = 0

where u and v represent the velocity components in the x and y directions, respectively.To eliminate the pressure term and ensure mass conservation, a scalar stream function (ψ)  is introduced such that:

u = ∂ψ/∂y and v = −∂ψ/∂x

The introduction of the streamfunction automatically satisfies the continuity equation. The vorticity (ω) is defined as the curl of the velocity field:

ω = ∂v/∂x − ∂u/∂y

Combining these definitions, the vorticity transport equation for two-dimensional incompressible flow is obtained as:

∂ω/∂t + u(∂ω/∂x) + v(∂ω/∂y) = ν(∂²ω/∂x² + ∂²ω/∂y²)

where ν is the kinematic viscosity and t denotes time. The relationship between stream function and vorticity is governed by the Poisson equation:

∇²ψ = −ω

where ∇² denotes the Laplacian operator, given by:

∇²ψ = ∂²ψ/∂x² + ∂²ψ/∂y²

To obtain the numerical solution, both equations are discretized using the finite difference method (FDM) on a uniform grid. The vorticity transport equation is solved explicitly using the second-order Runge–Kutta (RK2)  scheme for time advancement, ensuring temporal accuracy. The Poisson equation for streamfunction is solved iteratively using the Successive Over-Relaxation (SOR)  method to achieve rapid convergence at each time step. Boundary conditions are applied based on the physical geometry of the lid-driven cavity. The top wall moves with a uniform velocity (U = U_lid) while the side and bottom walls are stationary (U = 0, V = 0). The vorticity boundary conditions are evaluated using Thom’s formula , which relates wall vorticity to the stream function and known wall velocities.

Table 1: Mathematical Formulation Parameters.

This mathematical framework forms the foundation of the MATLAB-based CFD solver developed in this article, allowing accurate computation of stream function, vorticity, and velocity fields for a wide range of Reynolds numbers.

4.    Methodology:

The numerical methodology adopted in this research is based on the finite difference method (FDM) applied to the two-dimensional incompressible Navier–Stokes equations in the streamfunction–vorticity form. The computational domain is a square cavity discretized using a uniform Cartesian grid with equal spacing in both directions to ensure second-order spatial accuracy. Time advancement of the vorticity transport equation is achieved using an explicit second-order Runge–Kutta (RK2) method, which enhances numerical stability and accuracy compared to a first-order Euler scheme. Ferziger and Perić [14] demonstrated that accurate discretization schemes can significantly enhance convergence in laminar flow simulations.  The Lattice–Boltzmann approach offers an alternative CFD framework for solving complex boundary-driven flow problems [15]. The streamfunction equation, representing a Poisson-type problem, is solved iteratively using the Successive Over-Relaxation (SOR) technique to efficiently obtain the velocity field at each time step. Boundary conditions are applied using Thom’s formula, which ensures correct enforcement of the no-slip condition and realistic wall vorticity. The top wall of the cavity moves at a constant velocity, generating circulation within the domain, while the other walls remain stationary. The time step size is determined by the minimum of the advective and diffusive stability limits, ensuring the Courant–Friedrichs–Lewy (CFL) condition is satisfied. Each iteration involves computing the velocity components from the streamfunction, updating vorticity values, and solving for the new streamfunction field until steady state is reached. The computational loop continues until the changes in vorticity and streamfunction become negligibly small, indicating convergence.

Table 2: Simulation and Numerical Parameters.

The simulation is performed for a Reynolds number of 1000, where steady laminar flow behavior with distinct primary and secondary vortices is observed. MATLAB’s matrix-based environment is utilized to efficiently implement the iterative solver, allowing clear visualization of intermediate and final results. Six graphical outputs streamfunction contours, vorticity distribution, velocity vectors, horizontal and vertical velocity profiles, and a 3D streamfunction surface are generated to analyze flow characteristics. . Immersed boundary methods have also been developed to simulate fluid–structure interaction with moving walls [16].

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

Table 3: Numerical Scheme and Boundary Conditions.

The developed methodology ensures both numerical robustness and computational efficiency, serving as a reliable foundation for future CFD studies on more complex geometries or unsteady flow conditions.

5.    Design Matlab simulation and Analysis:

The MATLAB simulation developed for this research provides a robust numerical framework to solve the two-dimensional incompressible Navier–Stokes equations for the lid-driven cavity flow. The code is designed using modular programming, where each computational task such as vorticity calculation, velocity extraction, and Poisson equation solving is implemented as a separate function. The simulation begins by defining physical parameters, including Reynolds number, cavity dimensions, grid size, and boundary conditions. The computational grid is uniformly generated using MATLAB’s mesh utilities to ensure numerical consistency. The time-stepping process employs the second-order Runge–Kutta (RK2) method, which computes the vorticity field over each iteration based on convection–diffusion dynamics. The Poisson equation governing the streamfunction is solved iteratively using the Successive Over-Relaxation (SOR) technique ensuring convergence to a steady-state solution. Boundary vorticity is continuously updated using Thom’s formula to maintain accurate representation of wall shear effects. The simulation proceeds until the vorticity and streamfunction fields stabilize, indicating that a steady flow has been achieved. During execution, MATLAB’s plotting and contouring functions are used to visualize intermediate and final results, allowing real-time monitoring of flow development. The program outputs six key graphical results: streamfunction contours, vorticity contours, velocity vector plots, horizontal and vertical centerline velocity profiles, and a 3D streamfunction surface. Each plot reveals distinct flow characteristics such as primary vortex formation, corner eddies, and streamline curvature. The simulation results at Re = 1000 confirm the presence of a stable laminar circulation with secondary vortices, consistent with established benchmark data. Overall, the MATLAB-based simulation demonstrates the capability of high-level programming for efficient and accurate CFD modeling, serving as a foundation for further extension to transient or three-dimensional flows.

Figure 3 illustrates the streamfunction (( psi )) contours for the steady-state lid-driven cavity flow at a Reynolds number of 1000. The contour lines represent streamlines that indicate the flow paths within the cavity. The main feature is a large primary vortex occupying the central region, generated by the motion of the top lid. Near the bottom left and right corners, smaller secondary vortices are visible, signifying the recirculation zones caused by viscous effects and wall interactions. The smooth and continuous contour distribution confirms numerical stability and accurate convergence of the solver. The symmetry of the flow pattern about the cavity’s vertical centerline indicates that the discretization and boundary implementation are consistent. Overall this figure validates the streamfunction–vorticity formulation and provides a clear visualization of the steady laminar flow structure inside the cavity.

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

Figure 4 presents the vorticity contour map obtained from the numerical simulation of the cavity flow. The vorticity field highlights regions of intense rotational motion, especially near the moving lid and cavity corners. High vorticity magnitudes appear along the top wall, indicating strong shear caused by the lid motion. As the fluid circulates, secondary vortices form in the lower corners showing local rotational effects. The smooth contour distribution confirms numerical stability and convergence of the solution. The vorticity gradients also reveal the transition from strong shear to core flow regions. Overall, this figure validates the ability of the streamfunction–vorticity formulation to capture key flow dynamics accurately.

Figure 5 illustrates the velocity vector field that visualizes the flow circulation inside the cavity. The arrows represent the local direction and relative speed of the fluid particles, with denser arrows near the top indicating higher velocities due to the moving lid. The main vortex occupies the central region of the cavity, while smaller secondary vortices appear near the bottom corners. The smooth transition of velocity vectors demonstrates numerical accuracy and proper enforcement of boundary conditions. The clockwise rotation pattern clearly reflects the influence of the top-driven wall motion. This visualization effectively captures the flow structure and energy distribution within the confined domain.

Figure 6 presents the vertical velocity profile plotted along the horizontal centerline of the cavity. The curve illustrates how the vertical component of velocity changes from the left wall to the right wall. Near both sidewalls, the velocity approaches zero due to the no-slip boundary condition while the central region exhibits negative and positive peaks corresponding to the circulating motion of the main vortex. This variation indicates fluid rising near one wall and descending near the other consistent with the clockwise rotation of the cavity flow. The smooth and symmetric profile confirms the stability and accuracy of the numerical scheme. This figure is crucial for comparing simulation outcomes with benchmark CFD results.

Figure 7 illustrates the contour plot of the pressure distribution within the lid-driven cavity at steady state. The high-pressure region is observed near the upper left corner where the lid motion first interacts with the stationary wall, creating a compression zone. Conversely, a low-pressure area forms near the upper right corner, indicating fluid acceleration due to the moving lid. The central region exhibits a gradual pressure transition corresponding to the main vortex circulation. The smooth pressure gradient confirms the numerical stability and physical accuracy of the simulation. These results validate that the applied boundary conditions and numerical discretization correctly capture the incompressible flow behavior inside the cavity.

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

Figure 8 shows a 3D surface representation of the streamfunction field, providing a visual understanding of the flow circulation strength and pattern inside the cavity. Peaks and valleys on the surface correspond to regions of strong and weak rotational flow, respectively. The smooth curvature indicates steady-state convergence of the simulation. The highest gradient regions occur near the top wall, where the lid motion induces strong shear forces. The surface shape clearly depicts the dominant primary vortex at the center and weaker secondary vortices in the lower corners. This visualization effectively highlights the overall energy distribution and the dynamic behavior of the cavity flow.

6. Result and Discussion:

The simulation results for the two-dimensional lid-driven cavity flow at Re = 1000 demonstrate excellent agreement with established computational fluid dynamics (CFD) benchmarks. The streamfunction contours reveal a well-defined primary vortex occupying the central region of the cavity, accompanied by smaller secondary vortices near the lower corners, which are characteristic of high-Reynolds-number flows. The numerical approach employed here aligns with Chorin’s projection method for enforcing incompressibility constraints [17]. Turbulence modeling theories, such as those presented by Wilcox, provide the basis for extending this framework to high-Reynolds-number flows [18]. The vorticity distribution highlights intense shear layers near the top moving wall, confirming the strong rotational behavior induced by the lid motion. Velocity vector fields further illustrate the smooth and continuous circulation pattern, ensuring numerical stability and consistency throughout the domain. The horizontal and vertical velocity profiles along the centerlines align closely with literature data, indicating that the numerical discretization and boundary conditions were properly implemented. The SOR-based Poisson solver successfully maintained convergence within acceptable tolerance levels. The flow exhibits a steady-state structure with minimal numerical oscillations, validating the robustness of the RK2 integration approach. Pozrikidis emphasized the importance of verifying grid independence and numerical consistency in CFD studies [19]. Moreover, the 3D surface of the streamfunction provides an insightful visualization of flow intensity and vortex dynamics. Overall, the results confirm that the streamfunction–vorticity formulation, coupled with MATLAB implementation, provides a reliable and accurate framework for simulating incompressible laminar cavity flows.

7.    Conclusion:

The study successfully simulated the two-dimensional incompressible lid-driven cavity flow using the streamfunction–vorticity formulation in MATLAB. The developed numerical model efficiently captured the main flow structures, including the primary and secondary vortices, with high accuracy and numerical stability. Courant, Friedrichs, and Lewy [20] established the CFL condition, which dictates the stability limit for explicit time-stepping schemes used in this simulation. The results confirm that the explicit RK2 scheme, combined with the SOR Poisson solver, provides reliable convergence for steady-state solutions. The velocity and vorticity distributions match expected physical behavior and established CFD results. The approach demonstrates that MATLAB can serve as a powerful research tool for advanced computational fluid dynamics (CFD) investigations. Furthermore, the methodology is extendable to higher Reynolds numbers, transient analyses, or complex geometries. Overall, this work contributes a solid foundation for Ph.D.-level CFD research focusing on flow dynamics and numerical stability in confined domains.

8.    References:

[1] G. K. Batchelor, An Introduction to Fluid Dynamic, Cambridge University Press, 2000.

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

[3] J. H. Ferziger and M. Perić, Computational Methods for Fluid Dynamics, 3rd ed., Springer, 2002.

[4] P. Moin, Fundamentals of Engineering Numerical Analysis, 2nd ed., Cambridge University Press, 2010.

[5] S. V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, 1980.

[6] J. C. Tannehill, D. A. Anderson, and R. H. Pletcher, Computational Fluid Mechanics and Heat Transfer, 2nd ed., Taylor & Francis, 1997.

[7] G. D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford University Press, 1985.

[8] C. Pozrikidis, Introduction to Theoretical and Computational Fluid Dynamics*, Oxford University Press, 2011.

[9] H. K. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, 2nd ed., Pearson, 2007.

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

[11] D. F. Griffiths and D. J. Higham, Numerical Methods for Ordinary Differential Equations: Initial Value Problems, Springer, 2010.

[12] C. Hirsch, Numerical Computation of Internal and External Flows, Vol. 1, Elsevier, 1990.

[13] H. Schlichting and K. Gersten, *Boundary-Layer Theory, 9th ed., Springer, 2017.

[14] S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” J. Comput. Phys., vol. 103, no. 1, pp. 16–42, 1992.

[15] J. Kim and P. Moin, “Application of a fractional-step method to incompressible Navier–Stokes equations,” J. Comput. Phys., vol. 59, no. 2, pp. 308–323, 1985.

[16] U. Ghia, K. N. Ghia, and C. T. Shin, “High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method,” J. Comput. Phys. , vol. 48, no. 3, pp. 387–411, 1982.

[17] A. J. Chorin, “A numerical method for solving incompressible viscous flow problems,” J. Comput. Phys., vol. 2, no. 1, pp. 12–26, 1967.

[18] D. C. Wilcox, Turbulence Modeling for CFD, 3rd ed., DCW Industries, 2006.

[19] M. Griebel, T. Dornseifer, and T. Neunhoeffer, Numerical Simulation in Fluid Dynamics: A Practical Introduction, SIAM, 1998.

[20] M. P. Kirkpatrick, S. W. Armfield, and J. H. Kent, “A representation of curved boundaries for the solution of the Navier–Stokes equations on a staggered three-dimensional Cartesian grid,” J. Comput. Phys., vol. 184, no. 1, pp. 1–36, 2003.

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

Keywords: Lid-driven cavity flow, Computational Fluid Dynamics (CFD), streamfunction–vorticity formulation, incompressible Navier–Stokes equations, MATLAB simulation, Successive Over-Relaxation (SOR), Runge–Kutta method (RK2), Poisson solver, vorticity dynamics, laminar flow, boundary conditions, numerical stability, grid resolution, flow visualization, High Reynolds number simulation.

Do you need help with CFD Analysis of Incompressible Flow in a Lid-Driven Cavity in MATLAB? Don’t hesitate to contact our Tutors to receive professional and reliable guidance.