Numerical Methods for Incompressible Viscous Flows: A Comparative Review of Boundary Element, Domain Discretization, and Meshfree Approaches

1. Introduction

The numerical simulation of incompressible viscous flows governed by the Navier–Stokes equations remains a central problem in fluid mechanics and engineering. These equations describe the conservation of mass and momentum in Newtonian fluids and form the foundation of computational fluid dynamics (CFD), which has become an indispensable tool for the analysis of complex flow phenomena [1,2,3,4].

Over the past several decades, a wide range of numerical methods has been developed to approximate the Navier–Stokes equations. Among these, domain discretization methods, including finite difference (FDM), finite element (FEM), and finite volume (FVM) approaches, have emerged as the dominant frameworks for general-purpose CFD simulations due to their ability to handle nonlinearities, complex geometries, and turbulence modeling [5,6,7]. In parallel, boundary integral formulations, particularly the boundary element method (BEM), provide efficient alternatives for linear flow problems and unbounded domains by reducing the dimensionality of the problem [8,9,10].

More recently, meshfree methods have been introduced as flexible numerical techniques that avoid predefined mesh connectivity and instead rely on particle or node-based approximations. These methods are particularly attractive for problems involving large deformations, moving boundaries, and multiphase flows [11,12,13].

A fundamental distinction in fluid mechanics is between external flows, occurring around bodies in unbounded domains, and internal flows, occurring within confined geometries such as pipes and channels [14,15]. The choice of numerical method depends strongly on this classification, as well as on the Reynolds number and the associated flow physics.

The objective of the present review is to provide a systematic comparison of finite difference, finite element, boundary element, and meshfree methods for incompressible viscous flows of Newtonian fluids. Particular emphasis is placed on their mathematical formulation, computational characteristics, and compatibility with turbulence modeling approaches.

2. Historical Development of Numerical Methods

The development of numerical methods for incompressible viscous flows has progressed in parallel with advances in computational science and applied mathematics. Early computational approaches in the 1950s and 1960s focused on finite difference approximations of the Navier–Stokes equations. A key milestone was the marker-and-cell (MAC) method, which enabled stable numerical simulation of incompressible flows [16].

Subsequent developments addressed the coupling between velocity and pressure fields, leading to projection and fractional-step methods that remain widely used in modern CFD [17,18]. These advances were supported by the establishment of theoretical frameworks for stability and convergence in numerical approximations of partial differential equations [19,20].

During the 1970s and 1980s, the introduction of finite element and finite volume methods significantly expanded the capabilities of CFD. Finite element formulations provided flexibility in handling complex geometries through variational approaches, while finite volume methods ensured strict conservation properties [21,22,23,24]. These developments established domain-based methods as the primary tools for engineering simulations.

In parallel, boundary integral methods were developed, leading to the formulation of the boundary element method. These approaches exploit fundamental solutions to reduce the dimensionality of the problem and proved particularly effective for Stokes flows and unbounded domains [25,26,27].

More recently, meshfree methods have emerged as an alternative paradigm motivated by challenges in mesh generation. Particle-based and kernel-based approaches have enabled the simulation of flows involving large deformations and complex interfaces [28,29,30].

The historical evolution of CFD thus reflects a continuous interplay between mathematical formulation, computational capability, and the physical complexity of fluid flows.

In Table 1 a brief timeline regarding the key developments in CFD methods in the period 1950–2020 is presented

From a historical standpoint, the current predominance of domain discretization methods in CFD should therefore be understood not as a rejection of boundary-based techniques, but as the outcome of a natural co-evolution between physical modeling requirements and computational capability. Boundary element methods remain valuable within their domain of applicability; however, the expansion of CFD into turbulence-dominated and industrially relevant regimes has firmly established domain discretization methods as the principal computational paradigm.

Overall, the historical progression of CFD methodologies reflects a continuous effort to balance mathematical rigor, computational efficiency, and physical realism.

3. Mathematical Foundations of CFD Formulations

3.1. Governing Equations and Formulations

The physical aspects of any fluid flow are generally governed by the following fundamental principles:

(1)

Mass conservation

(2)

Energy conservation

(3)

Newton’s second law for momentum conservation

(4)

Chemical species conservation

These aforementioned principles are represented in terms of basic mathematical equations which in their most general form are either integral or partial differential equations. In this framework, CFD aims at discretizing these equations to algebraic formulations, which are then solved to obtain numbers of the flow-field values, at discrete points in space and time, which evidently correspond to the dependent variables of the fluid flow. Numerical methods differ primarily in how these equations are reformulated and discretized. A key distinction is between the strong form, in which the equations are satisfied pointwise, and the weak form, obtained by integrating the equations against test functions over the domain. Weak formulations reduce regularity requirements and form the basis of finite elements and many meshfree methods.

Numerical methods differ primarily in how these equations are reformulated and discretized.

In fact, all governing conservation equations can be written out in the following general form:

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho \cdot \phi \right)+div\left(\rho \cdot u\cdot \phi \right)-\Gamma \phi \cdot grad\phi =S\phi$$

(1)

where ρ is the density of air or mixture (if chemical species are present), t is the time, u is the velocity vector, Γφ is the diffusion coefficient of the transported quantity, Sφ is the source or sink term and φ is the general transported quantity. In general, these methods can be divided into three categories [1,2,3,4,5,6,7,8,9,10]: Finite-Difference Method (FDM). This method is widely used in CFD applications due to its computational and conceptual simplicity. According to the FDM a flow-field variable is discretized in terms of Taylor’s series expansion about a (grid) point defined by the variable’s either-side derivatives. Finite-Volume Method (FVM). The algebraic difference equations are formulated by applying the conservation laws directly to a small control volume (grid-cell). In opposition with FDM, it has the advantage of the flexibility of the boundary shape of the computational cells, as the conservation laws apply irrespective of shape.Finite-Element Method (FEM). The method is formulated so that it minimizes the weighted residual (Galerkin residuals) over the flow field. As this method uses unstructured meshes consisting of cells of different topologies (triangles or quadrilateral for two dimensional meshes), very complex geometries can be handled adequately. A direct comparison of these methods is illustrated in Table 2.

3.2. Strong and Weak Formulations

A fundamental distinction exists between strong formulations, in which the governing equations are satisfied pointwise, and weak formulations, obtained through weighted residual or variational approaches. Weak forms reduce regularity requirements and underpin finite element and many meshfree methods [31,32].

3.3. Domain Discretization Methods

Finite difference methods approximate derivatives directly on structured grids and are efficient for simple geometries [19,20]. Finite element methods rely on variational formulations and provide geometric flexibility [21,31], while finite volume methods enforce conservation laws over control volumes [23]. On the other hand, the establishment of Finite Volume Method has been made on the basis of the following fundamental steps:

(1)

The calculation field is divided into a number of control volumes each enclosing a node.

(2)

The differential equation describing the transfer phenomenon is completed on a control volume.

(3)

The control volume consists of 4 surfaces (fronts).

(4)

The discretization of the transport equation is obtained according to the following assumptions:

a)

Uniform distribution of the different quantities in the control volume.

b)

Uniform distribution of different quantities on the control volume fronts.

c)

First-order approach with ascending differences in the term of time derivative.

On the other hand, referring to the Finite Element Method (FEM), we have to elucidate that a Delaunay triangulation is the straight-line dual of a Voronoi diagram [33]. In the sense of this rigorous mathematical statement, a one to one mapping between a Voronoi diagram and a Delaunay triangulation constitutes a Sine Qua Non Condition. This fact guarantees that the resultant truss after the triangulation of a surface cannot behave as an infinitesimal mechanism from the Kinematics viewpoint [34]. Actually, according to the concepts of FEM, Delaunay Triangulation and Voronoi Diagrams are the geometric pillars used to discretize a continuous domain into a computational mesh. Indeed, their influence is primarily seen in mesh quality, numerical stability together with the development of advanced element types. In this context, it can be said that both Voronoi diagrams and Delaunay triangulations are mathematical duals and they are linked to each other, in order to "smooth" a mesh.

In Table 3, the role of Delaunay Triangulation along with Voronoi Diagrams in Finite Element Method (FEM) are signified:

3.4. Boundary Element Methods

Boundary element methods reformulate the governing equations into boundary integrals using fundamental solutions, reducing problem dimensionality. This approach is particularly effective for linear problems but is limited in handling nonlinear and turbulent flows [25,26].

3.5. Meshfree Methods

Meshfree methods approximate the solution using scattered nodes or particles, enabling flexibility in handling evolving geometries. These methods retain a domain-based character while avoiding mesh generation [28,29,30].

3.6. Consistency, Stability, and Convergence

Numerical methods must satisfy consistency, stability, and convergence requirements. The Lax equivalence theorem establishes that consistency and stability imply convergence for linear problems [19,20]. A scheme is consistent if its discretization error vanishes as the grid is refined, and stable if numerical errors remain bounded during computation. The Lax equivalence theorem establishes that, for linear problems, consistency and stability together imply convergence [20].

According to the above, one may elucidate that in the context of contemporary CFD practice, the choice between boundary-based and domain-based numerical methods should not be framed as a competition between equivalent alternatives. Instead, it reflects a problem-driven selection governed by the mathematical structure of the governing equations and the physical phenomena of interest.

For linear, low-Reynolds-number, or potential flow problems, BEM remains an efficient and highly accurate choice. For turbulent, nonlinear, and industrially relevant flows, domain discretization methods—particularly finite volume and finite element formulations—are not merely preferable but effectively indispensable.

4. Turbulence Modeling and Its Interaction with Discretization

Turbulence modeling plays a central role in CFD and is inherently linked to the underlying numerical discretization. Approaches such as Reynolds-averaged Navier–Stokes (RANS), large-eddy simulation (LES), and direct numerical simulation (DNS) introduce different levels of modeling and resolution requirements [35,36,37].

Domain discretization methods provide a natural framework for turbulence modeling, as turbulent transport processes are volumetric and nonlinear. Finite volume and finite element methods are widely used in RANS and LES simulations due to their conservation properties and flexibility [5,35].

In contrast, the application of turbulence modeling within boundary element frameworks is limited, since turbulent stresses cannot be represented solely through boundary integrals [8,9,10]. Meshfree methods have been extended to turbulence simulation, particularly through SPH-based approaches, although challenges remain in accurately resolving wall-bounded turbulence [11,28].

The advent of Large-Eddy Simulation, introduced conceptually by Smagorinsky, and the continued development of Direct Numerical Simulation (DNS) further emphasized the importance of accurate domain-wide discretization. DNS studies initially relied heavily on structured finite difference methods for their spectral-like accuracy, while LES benefited from both FDM and FEM, particularly with variational multiscale interpretations linking discretization error to subgrid modeling.

By contrast, the integration of turbulence modeling into classical BEM proved inherently limited. Since the nonlinear convective term and turbulent stress transport are volumetric, exact boundary-only reduction is not feasible. Although hybrid boundary–domain approaches exist, BEM has remained most competitive in laminar or creeping-flow regimes rather than high–Reynolds-number turbulence.

Meshfree methods have also been extended to turbulent flows, including weakly compressible SPH-based LES formulations. Their Lagrangian nature offers advantages for free-surface and strongly transient flows, yet accurate enforcement of incompressibility and wall-bounded turbulence remains an active research frontier.

The distinction between external and internal flows provides a physically meaningful framework for comparing numerical methods.

External flows involve fluid motion past bodies in effectively unbounded domains. At low Reynolds numbers, such flows reduce to Stokes problems for which boundary integral formulations are particularly natural and efficient. In creeping external flows, BEM offers high accuracy in computing hydrodynamic forces and stresses without discretizing the far field. However, as Reynolds number increases and boundary layers, separation, and wake turbulence emerge, full volumetric discretization becomes unavoidable. Domain-based methods thus dominate in aerodynamic and hydrodynamic simulations at moderate and high Reynolds numbers.

Internal flows occur within confined geometries—pipes, channels, cavities, and process equipment—where pressure gradients drive volumetric motion. Even in laminar regimes, the physics are inherently domain-wide. At high Reynolds numbers, wall-bounded turbulence and secondary flows further necessitate volumetric resolution. Finite difference methods remain efficient in canonical geometries, while FEM and related methods provide flexibility in complex internal configurations. Meshfree approaches offer potential advantages when boundaries deform or move, but their competitiveness in canonical wall turbulence remains under evaluation.

Hence, suitability of a numerical method depends not solely on geometric classification but on the interplay between geometry, Reynolds number, linearity, and the need to resolve volumetric transport processes.

Domain discretization methods gained prominence not merely due to computational convenience, but because their mathematical structure aligns naturally with the nonlinear, volumetric character of turbulent Navier–Stokes flows. Boundary element methods retain distinctive advantages in linear, low–Reynolds-number, and unbounded-domain settings. Meshfree methods provide geometric flexibility and Lagrangian advantages in strongly transient problems.

Thus, the interaction between turbulence modeling and discretization is a key factor in determining the applicability of numerical methods.

5. Classification of Numerical Methods for Incompressible Flows

Numerical methods for incompressible flows can be classified into three main categories:

• Domain discretization methods (FDM, FEM, FVM)
• Boundary-Element methods (BEM)
• Meshfree methods

This classification reflects differences in mathematical formulation and discretization philosophy. Domain methods discretize the entire flow field, boundary methods reduce dimensionality, and meshfree methods eliminate mesh constraints while retaining volumetric representation [6,25,28].

Additionally, methods can be classified according to flow regime, including laminar versus turbulent flows and internal versus external flows [14,15].

6. Comparative Assessment of Numerical Methods

Each class of numerical methods exhibits distinct strengths and limitations. Domain discretization methods provide general applicability and are well suited to nonlinear and turbulent flows. Finite volume methods ensure conservation, while finite element methods offer geometric flexibility [5,21].

Boundary element methods are highly efficient for linear problems and unbounded domains but are limited in handling nonlinearities and turbulence [25].

Meshfree methods offer flexibility for complex and evolving geometries but face challenges in achieving high accuracy and computational efficiency in wall-bounded flows [28,29,30].

The choice of method depends on the interplay between physical modeling requirements, geometric complexity, and computational resources.

This comparison underscores that no single method is universally optimal. Instead, the choice of method depends on the interplay between flow physics, geometric complexity, and computational constraints. Table 4 illustrates a mathematical comparison of principal CFD simulations.

7. Recent Advances in Numerical Methods of CFD

In the last decades, we have witnessed substantial progress in the development and application of numerical methods for incompressible viscous flows, driven by advances in high-performance computing, algorithmic design, and increasing demands for predictive accuracy in complex flow environments. These developments have further clarified the relative strengths of finite difference, finite element, boundary element, and meshfree approaches, while also highlighting emerging trends toward hybrid and high-order methodologies.

Recent developments in computational fluid dynamics have further clarified the relative strengths of domain discretization, boundary-based, and meshfree methods across a wide range of flow regimes [37,38]. High-fidelity simulations based on finite volume and finite element formulations remain the dominant tools for turbulent flow prediction, particularly in large-eddy simulation (LES) and direct numerical simulation (DNS), where scalability and conservation properties are critical [39,40,41]. In this context, modern high-order discretization schemes and adaptive mesh refinement techniques have significantly improved accuracy and efficiency in resolving complex flow structures.

Advances in smoothed particle hydrodynamics (SPH) and related methods have improved consistency, boundary treatment, and turbulence modeling capabilities. These developments have expanded the applicability of meshfree methods to challenging problems such as free-surface flows, multiphase systems, and flows with complex moving boundaries. Additionally, new kernel correction techniques and coupling strategies with grid-based methods have enhanced their robustness and accuracy [42,43].

Meshfree methods have experienced renewed interest, particularly in applications involving free-surface flows, multiphase systems, and fluid–structure interaction. Recent studies demonstrate that smoothed particle hydrodynamics (SPH) and related particle-based methods can successfully capture highly nonlinear and transient flow phenomena, although challenges remain in achieving consistent accuracy near solid boundaries and in wall-bounded turbulence [43,44,45,46].

At the same time, hybrid and immersed boundary approaches have emerged as powerful alternatives for handling complex geometries within domain-based frameworks. These methods combine the robustness of Eulerian discretizations with the flexibility of Lagrangian representations, enabling efficient simulation of moving and deforming bodies in viscous flows [38,40].

Boundary element methods continue to play an important role in specialized applications, particularly in low-Reynolds-number hydrodynamics and multiphysics problems involving coupled fields. Recent advances have focused on accelerating BEM through fast algorithms and coupling it with volumetric methods to overcome limitations associated with nonlinear flow regimes.

In the meanwhile, in Ref. [47] a novel method was developed to simulate incompressible flows for the purpose of predicting green water events. In particular this method is: meshless, Lagrangian, volume–conservative, second–order accurate in space, efficient, and suitable for coupling. In this context a set of novel spatial operators based on the weighted–least squares, which are used to describe and solve the Navier–Stokes equations in strong form. Volume–conservative Lagrangian advection was used, which naturally handles violent free–surface flows. Boundary conditions are conforming to moving geometry at each time step. This makes coupling with other mesh–based and structural solvers straightforward.

Moreover, in Ref. [48] a high-accurate meshless formulation for solving the incompressible Navier-Stokes equations with the weakly compressible approach was developed.

Further, in a valuable work presented in Ref. [49] a fully meshless approach was applied to the simulation of heat conduction problems. The code has been verified according to the Method of Manufactured Solutions, and it is proven to be reliable and suitable for practical applications, even with convective boundary conditions. Indeed, the results presented in this paper exceeded the most optimistic expectations, especially on complex domains.

Finally, in Ref. [50] a generalized multiscale finite element method was combined with a balanced truncation, towards the numerical treatment of a class of parameter-dependent parabolic problems.

Overall, one may deduce that the recent advances in CFD, definitely highlight a clear shift toward high-order accuracy, hybridization, and improved turbulence–discretization coupling, reflecting the growing complexity of fluid dynamics problems and the need for robust, scalable computational tools.

8. Future Perspectives: Motivations, Prospects, and Challenges

The continued evolution of computational fluid dynamics (CFD) for incompressible viscous flows is driven by the increasing demand for predictive accuracy, computational efficiency, and the capability to address complex, multiscale, and multiphysics phenomena. While finite difference, finite volume, finite element, boundary element, and meshfree methods have each reached a high level of maturity, significant challenges and opportunities remain. Future progress is expected to arise from advances in numerical analysis, high-performance computing, data-driven modeling, and the integration of hybrid methodologies.

8.1. Fundamental Motivations

A central motivation for future developments in CFD lies in the need to bridge the gap between mathematical modeling and physical realism. Specifically, in incompressible flows, this concept is particularly evident in:

• High–Reynolds-number turbulence, where fully resolved simulations (DNS) remain computationally prohibitive for most practical applications
• Multiphase and free-surface flows, involving interface dynamics and topological changes
• Fluid–structure interaction (FSI) and moving boundary problems
• Complex geometries and industrial-scale configurations
• These challenges demand numerical methods that are not only accurate and stable, but also scalable and adaptable across a wide range of flow regimes.

8.2. Prospects for Domain Discretization Methods (FDM, FVM, FEM)

a) 

Finite Difference Methods (FDM)

Future developments in FDM are expected to focus on:

• High-order and spectral-like schemes, improving resolution of turbulent scales
• Adaptive structured grids, enabling local refinement while preserving efficiency
• GPU-accelerated implementations, enhancing performance for large-scale DNS and LES
• Despite their geometric limitations, FDM will likely remain highly relevant in canonical turbulence studies and fundamental research.

b) 

Finite Volume Methods (FVM)

FVM is expected to continue as the dominant industrial CFD framework, with future directions including:

• Advanced flux reconstruction and high-order schemes
• Improved coupling with turbulence models (LES, hybrid RANS–LES)
• Adaptive mesh refinement (AMR) for multiscale flows
• Robustness in complex multiphysics environments
• A major challenge for FVM lies in achieving high-order accuracy on unstructured meshes while maintaining strict conservation properties.

c) 

Finite Element Methods (FEM)

FEM is likely to play an increasingly important role in:

• Multiphysics simulations (e.g., fluid–structure interaction, biofluid mechanics)
• High-order and isogeometric analysis, improving geometric fidelity and solution accuracy
• Variational multiscale (VMS) turbulence modeling frameworks
• Stabilized formulations for convection-dominated flows
• The primary challenges include computational cost, solver scalability, and the development of robust formulations for high-Reynolds-number turbulence.

8.3. Prospects for Boundary Element Methods (BEM)

Boundary element methods are expected to remain highly specialized but valuable in:

• Low-Reynolds-number and Stokes flow regimes
• Unbounded and exterior flow problems
• Multiphysics applications (e.g., electrohydrodynamics, microfluidics)
• Future advances will likely focus on:
• Fast algorithms (fast multipole methods, hierarchical matrices)
• Hybrid BEM–volume coupling strategies
• Extension to weakly nonlinear flows

However, a fundamental limitation persists: the inability of classical BEM formulations to naturally handle nonlinear and turbulent volumetric phenomena. Overcoming this limitation remains a major open challenge.

8.4. Prospects for Meshfree Methods

Meshfree methods are expected to expand significantly in areas where traditional mesh-based approaches face difficulties:

• Free-surface and multiphase flows
• Large deformation and fragmentation problems
• Moving boundary and interface dynamics
• Future research directions include:
• Improved consistency and convergence properties
• Accurate boundary condition enforcement
• Integration with turbulence modeling frameworks
• Coupling with grid-based methods (hybrid particle–mesh approaches)

A key challenge is achieving competitive accuracy and efficiency in wall-bounded turbulent flows, where mesh-based methods currently dominate.

8.5. Hybridization and Method Integration

One of the most promising directions in CFD is the development of hybrid numerical frameworks, combining the strengths of different methodologies:

• Immersed boundary methods bridging Eulerian and Lagrangian descriptions
• Particle–mesh coupling techniques
• BEM–FEM or BEM–FVM hybrid formulations
• Multiresolution and adaptive methods

Such approaches aim to overcome the limitations of individual methods by leveraging their complementary advantages. Hybridization is expected to play a central role in next-generation CFD solvers.

8.6. Turbulence Modeling and Discretization Coupling

Future progress in CFD will depend critically on the interaction between turbulence modeling and numerical discretization. Key directions include:

• Wall-resolved and wall-modeled LES for complex geometries
• Scale-aware and adaptive turbulence models
• Data-driven and machine-learning-assisted turbulence closures
• Consistent coupling between discretization errors and subgrid-scale models
• Ensuring compatibility between turbulence models and discretization schemes remains a central challenge, particularly in high-order and meshfree frameworks.

8.7. Emerging Trends and Open Challenges

Several overarching trends are shaping the future of CFD:

• Exascale computing and parallel scalability
• High-order accurate methods for complex geometries
• Data-driven and physics-informed machine learning approaches
• Uncertainty quantification and verification/validation frameworks
• At the same time, key open challenges remain:
• Achieving predictive accuracy in turbulent flows at realistic Reynolds numbers
• Developing robust and efficient hybrid methods
• Ensuring numerical stability and consistency across scales
• Bridging the gap between academic methods and industrial applications

8.8. Concluding Perspective

In conclusion, the future of incompressible CFD will not be defined by a single dominant numerical method, but rather by the integration of multiple approaches, each tailored to specific physical and computational requirements. Domain discretization methods will remain central for turbulent and general-purpose simulations, boundary element methods will retain niche but critical roles in linear and unbounded problems, and meshfree methods will continue to expand into complex, evolving flow regimes. The convergence of high-order discretization, hybrid methodologies, and advanced turbulence modeling is expected to define the next generation of CFD tools, enabling increasingly accurate and efficient simulations of complex fluid systems.

9. Discussion

The comparative assessment of numerical methods for incompressible viscous flows highlights fundamental differences in their mathematical formulation, computational characteristics, and applicability across flow regimes. In particular, finite element, boundary element, and meshfree methods represent three distinct yet complementary approaches within the broader CFD landscape.

The finite element method (FEM) has become one of the most versatile and widely adopted techniques for the simulation of Newtonian flows, owing to its rigorous variational foundation and flexibility in handling complex geometries. Its formulation in weak form enables the natural incorporation of boundary conditions and facilitates the use of unstructured meshes, which is particularly advantageous in internal flows and multiphysics problems. Moreover, the development of stabilized formulations has significantly improved the robustness of FEM for convection-dominated and incompressible flows. In the context of turbulence modeling, FEM has been successfully applied in Reynolds-averaged Navier–Stokes (RANS) and large-eddy simulation (LES) frameworks, particularly when combined with high-order discretizations and adaptive refinement strategies.. However, the method may involve higher implementation complexity and computational cost compared to structured-grid approaches.

In contrast, the boundary element method (BEM) is based on integral representations of the governing equations and offers a fundamentally different computational paradigm. By reducing the dimensionality of the problem, BEM provides significant advantages for flows in unbounded domains and for problems governed by linear operators, such as Stokes flow. In such regimes, BEM enables accurate computation of hydrodynamic forces without discretizing the entire flow domain. Nevertheless, the applicability of BEM to general incompressible Navier–Stokes flows remains limited due to the nonlinear and volumetric nature of convective and turbulent transport processes. Recent developments have focused on hybrid approaches and fast algorithms, including fast multipole methods, to improve computational efficiency and extend applicability to more complex problems. Despite these advances, BEM remains primarily a specialized tool for low-Reynolds-number and multiphysics applications.

Meshfree methods provide an alternative framework that combines aspects of domain discretization with enhanced geometric flexibility. By eliminating the need for predefined mesh connectivity, these methods are particularly well suited to problems involving large deformations, moving interfaces, and free-surface flows. Particle-based approaches such as smoothed particle hydrodynamics (SPH) have been successfully applied to a range of fluid dynamics problems, including multiphase and transient flows. Their Lagrangian nature allows natural tracking of interfaces and complex boundary motion. However, challenges remain in achieving high accuracy near solid boundaries and in modeling wall-bounded turbulence, which limits their widespread adoption in high-Reynolds-number engineering applications. Ongoing research aims to address these limitations through improved kernel functions, boundary treatments, and hybrid formulations.

From a broader perspective, the comparison of these methods underscores that no single numerical approach is universally optimal. Instead, the suitability of a given method depends on the interplay between flow physics, geometric complexity, and computational requirements. FEM provides a general and robust framework for complex, nonlinear flows; BEM offers high efficiency in linear and unbounded-domain problems; and meshfree methods introduce flexibility for problems involving evolving geometries.

A notable trend in recent CFD research is the increasing development of hybrid methodologies, which seek to combine the strengths of different approaches. Examples include coupled boundary–volume formulations and immersed boundary techniques that integrate Lagrangian and Eulerian descriptions. Such approaches reflect the broader evolution of CFD toward flexible, multi-method frameworks capable of addressing the growing complexity of modern fluid dynamics problems.

On the other hand, Artificial intelligence (AI), and in particular machine learning (ML), is emerging as a transformative tool in computational fluid dynamics, with the potential to fundamentally enhance traditional numerical methods such as finite difference, finite volume, finite element, boundary element, and meshfree approaches. One of the most promising directions is the development of data-driven turbulence models, where neural networks are trained on high-fidelity datasets (e.g., DNS) to improve or replace classical closure models. These approaches aim to overcome longstanding limitations of RANS and LES modeling by learning complex, nonlinear relationships that are difficult to represent analytically. Importantly, such models can be embedded within existing discretization frameworks, enabling a hybrid paradigm in which physics-based solvers are augmented by data-driven components.

Beyond turbulence modeling, AI is increasingly being used to enhance the numerical performance and efficiency of CFD methods. For example, machine learning techniques have been applied to accelerate iterative solvers, optimize mesh generation and adaptation, and construct reduced-order models for rapid prediction of flow fields. In finite element and finite volume contexts, AI-assisted adaptivity can guide mesh refinement based on learned error indicators, while in meshfree methods, it can improve particle distribution and kernel selection. Similarly, in boundary element methods, data-driven surrogates may help approximate Green’s functions or accelerate the evaluation of integral operators, thereby reducing computational cost.

Despite these promising developments, significant challenges remain in integrating AI with classical CFD methodologies. Ensuring physical consistency, generalization across flow regimes, and interpretability of learned models is critical, particularly in high-stakes engineering applications. Moreover, the coupling between data-driven components and numerical discretization must be carefully designed to preserve stability and convergence properties. Nevertheless, the convergence of AI and CFD represents a powerful emerging direction, with the potential to redefine how complex fluid flows are modeled, simulated, and understood in the coming decades.

10. Conclusions

This review has presented a comprehensive examination of numerical methods for the simulation of incompressible viscous flows, with particular emphasis on domain discretization approaches, boundary element methods, and meshfree techniques. By tracing their historical development, mathematical foundations, and computational characteristics, a unified perspective has been established that highlights both their individual capabilities and their relative limitations.

A central conclusion is that domain-based methods remain the most general and widely applicable class of numerical techniques for incompressible Navier–Stokes equations. Finite difference, finite volume, and finite element methods provide robust frameworks for handling nonlinearities, complex geometries, and turbulence modeling. Their compatibility with established simulation strategies such as Reynolds-averaged Navier–Stokes, large-eddy simulation, and direct numerical simulation has solidified their role as the primary tools in both academic research and industrial applications. However, their effectiveness is often constrained by the need for high-quality mesh generation and the significant computational cost associated with high-fidelity simulations.

Boundary element methods, in contrast, offer a fundamentally different approach based on boundary-only discretization. Their principal advantage lies in the reduction of problem dimensionality and their ability to achieve high accuracy for linear flow problems, particularly in unbounded domains. Despite these strengths, their applicability to general incompressible flows is inherently limited by the nonlinear nature of the governing equations. As a result, their use is largely confined to specialized problems or hybrid formulations, rather than serving as a general-purpose CFD framework.

Meshfree methods represent an evolving class of numerical techniques that address some of the intrinsic limitations of mesh-based approaches. By eliminating the need for predefined mesh connectivity, they provide exceptional flexibility for problems involving large deformations, moving boundaries, and free-surface flows. These features make them particularly attractive for emerging applications in multiphase flows and fluid–structure interaction. Nevertheless, issues related to numerical stability, consistency, and computational efficiency continue to limit their widespread adoption, and further development is required for them to reach the level of maturity of traditional methods.

The comparative analysis conducted in this work underscores that the selection of an appropriate numerical method is inherently problem-dependent. Factors such as flow regime, geometric complexity, boundary conditions, and computational resources must all be carefully considered. No single method offers a universally optimal solution; rather, each approach provides a different balance between accuracy, efficiency, and flexibility.

An important overarching theme is the increasing relevance of hybrid and interdisciplinary approaches. The integration of different numerical methodologies, as well as the incorporation of emerging technologies such as data-driven modeling and high-performance computing, is reshaping the landscape of computational fluid dynamics. These developments suggest that future progress will not be driven by a single dominant method, but rather by the strategic combination of complementary techniques.

In conclusion, numerical simulation of incompressible viscous flows remains a dynamic and evolving field, characterized by both mature methodologies and ongoing innovation. Continued advances in numerical analysis, algorithm design, and computational technology are expected to further enhance the accuracy, efficiency, and applicability of CFD methods. The insights provided in this review aim to support informed methodological choices and to contribute to the development of next-generation computational tools for fluid dynamics.