NUMERICAL METHODS IN FLUID MECHANICS – AN OVERVIEW

The article presents in a review way the most important numerical methods used in modern fluid mechanics. The individual chapters discuss Finite Difference Method, Finite Volume Method, Lattice Boltzmann Method, Discrete Element Method and Smoothed Particle Hydrodynamics. The aim of the article is to familiarize the reader with the most important concepts, features and mathematical equations used in particular methods. The article is intended mainly for people who want to get acquainted with the current possibilities of numerical modelling in the field of broadly understood fluid mechanics. The material is intended to facilitate the decision on how to implement the planned research.


Introduction
In the field of science, regardless of the specific discipline, numerical modelling has become one of the essential methods for investigating or predicting the behaviour of various systems.The primary challenge is the extensive knowledge required to effectively utilize numerical modelling.Over the years, a wide range of techniques has been developed to address various classes of problems.
A typical numerical modelling process consists of the following stages: -preparation of the numerical model for calculations (pre-processing): a) determining the geometry of the computational domain, b) discretizing the computational domain, c) selecting an appropriate mathematical model, d) defining a closure set, e) specifying properties of materials used, f) setting boundary conditions, g) establishing initial conditions; -performing calculations (solving): a) defining numerical techniques, b) determining monitoring parameters for the calculation process, c) configuring recording settings, d) defining termination criteria for the calculations, e) executing simulations and monitoring the progress; -processing and analysis of the resulting data (post-processing): a) visualization of the simulation results, b) in-depth analysis of the data, c) data processing to extract meaningful information and insights.
In the context of numerical analysis, a computational domain signifies a virtual representation of the geometry, encompassing either its entirety or a partial depiction of a real-world system.This representation is customtailored to suit the specific demands of numerical simulations.The term "closure" in numerical modelling refers to a mathematical model that elucidates a specific phenomenon or process.It augments the fundamental set of mathematical equations, enabling them to be effectively solved.The 13 stages of numerical modelling may seem straightforward at first glance, but each stage entails a wealth of extensive and intricate knowledge.Consider, for instance, the process of discretizing the computational domain.This process encompasses various methods, including the utilization of nodes, pixels, voxels, elements, or cells.These objects can form structured or unstructured meshes and may be node-centered, cell-centered, multi-nodal, or multi-cellular.Moreover, the selection of discretization can be further categorized into structured, unstructured, hybrid, multi-block, multi-zone, overlapping, adaptive, movable (with rotary or reciprocating movement), and various other techniques.Each type of domain discretization necessitates distinct generation techniques and algorithms, each characterized by its unique aspects.
Another example is the algorithm required for solving sets of linear equations, a common necessity during the solving stage.While developing numerical methods, it became apparent that classical methods typically described in textbooks were inadequate.New, significantly more intricate techniques had to be devised, such as relaxation methods, methods founded on Krylov space, or multi-grid methods.Another challenge was determining the most effective way to implement these methods.The implementation process, in turn, drove the advancement of parallel programming and computational techniques that could harness the power of graphics processors.It is evident that contemporary numerical modelling represents a fusion of knowledge from various fields, including mathematics, physics, algorithmic, and programming.
It is important to acknowledge that various numerical methods were developed during different eras, each with its own set of technological circumstances, including operating systems and programming languages.Over the years, programming languages like FORTRAN and C were commonly utilized, and later on, C++ gained popularity.In the present day, Python has also emerged as a widely-used language.Python can be employed directly in the implementation of numerical methods or indirectly as an environment that connects various software or libraries, often written in FORTRAN, C, or C++.These languages continue to dominate the field of numerical modelling due to their high-speed computational capabilities, which make them well-suited for complex mathematical operations.
All of the factors mentioned above, when considered collectively, make it challenging for a single researcher to master all the numerical techniques currently in use.Typically, researchers tend to specialize in one particular numerical method, especially when their research is focused on a specific type of matter, such as solids, fluids, or granular materials.However, in cases involving multiphase systems that demand the integration of different numerical techniques, researchers may employ multiple methods.This high degree of specialization often results in limited knowledge of alternative approaches, even those that could be potentially applied within the same research area.It is important to note that this is a generalized perspective and may not account for unique situations or exceptions.
Having a broader understanding allows for more strategic planning and the more efficient allocation of resources to achieve specific objectives.The primary goal of the research described in this paper is to introduce selected numerical methods, which are frequently employed in the realm of Computational Fluid Dynamics (CFD).This introduction encompasses the fundamental principles underlying these methods, the primary equations they employ, their key characteristics, areas of application, as well as their strengths and limitations.The content is designed to cater to two primary groups of readers: those who are new to numerical modelling and those who possess some experience in one or two specific areas and wish to expand their knowledge in this domain.
It should be emphasized that the starting point for the considerations presented in the article is the Finite Volume Method, which is described in the most detail for this reason.Other methods either serve as potential alternatives (such as the Lattice Boltzmann Method and the Smoothed Particle Hydrodynamics) or complement it (such as the Discrete Element Method).The Finite Difference Method is described because it is historically the first numerical method used in fluid mechanics and is still used there in some areas today.Some of the methods described in the article, especially the Discrete Element Method and the Smoothed Particle Hydrodynamics, can be used to model fluid-solid systems or materials located at the interface between solids and fluids.In this area, they often serve as better tools than methods typically dedicated to fluids.
A key to understanding the choice of numerical methods described in the article is presented further in Figure 5, summarizing the strategies for modelling multiphase media.All the numerical methods described in the article can be used independently or in combination to implement one of the visible variants.It should be noted that the aforementioned figure does not include solid matter, which explains why the Finite Element Method is not described in the article.Of course, there are issues where this method is coupled with a method designed for fluids.Examples include fluid flow in flexible pipes or shells, fluid interaction with deformable membranes, etc.It should be added that the Finite Element Method can also be used to directly solve flow problems, as evidenced by the capabilities of programs such as SolidWorks, COMSOL, or Abaqus, but it has never been the focus of the authors, who consider the Finite Volume Method to be more suitable for modelling fluid flows than the Finite Element Method.The article only describes numerical methods with which the authors had direct practical experience, not just theoretical knowledge.For the same reason, many other methods or their variants, such as the Discontinuous Galerkin Method, the Boundary Element Method, or the Immersed Boundary Method, are not described.For those interested in the application of the Finite Element Method in the analysis of fluid flow problems, it is recommended to explore publications such as HugHes and ZienkiewicZ (1979), weinan and Liu (2000) or ZienkiewicZ et al. (2005).

Finite Difference Method
The Finite Difference Method (FDM) (courant et al. 1928, causon, MingHaM 2010, Langtangen, Linge 2017) is an Eulerian method that involves approximating the derivative of a function using finite differences.It relies on a previously defined spatial discrete mesh of nodes.The starting point is the Newton-Raphson method, where the derivative of the function at a point equals the slope of the tangent line to the function (curve) at that point (iyengar, Jain 2009, Burden, Faires 2011).This concept is illustrated in Figure 1.However, this definition is not unambiguous because it can be implemented in at least three ways: the so-called forward scheme (1), the backward scheme (2), or the central scheme (3): ′ () ≈ (  + Δ) − (  ) Δ (1) ′ () ≈ (  ) − (  − Δ) Δ (2) ′ () ≈ (  + Δ) − (  − Δ) 2Δ (3) In equations (1-3), only three nodes are utilized (i -1, i, i + 1).However, higher--order numerical schemes can also be structured, involving 5 or more nodes to approximate a derivative.The second derivative of the function f(x) at the point x 0 is obtained by applying the difference scheme three times.Depending on the schemes employed, there can be multiple variations of the second derivative.
It should be noted that if the equation includes a time-dependent term, time must also be discretized.The maximum allowable time step depends on the numerical scheme utilized and is sometimes constrained by additional parameters, primarily the CFL condition.T x -coordinate [m].
By applying the central difference scheme for the second derivative and the explicit Euler scheme for time discretization, equation (4) takes the following form (Langtangen, Linge 2017): where: n and n + 1 represent the current and the next time levels, respectively.
Equation ( 5) is valid for nodes from to N -1, where N is the number of nodes in the numerical grid.Boundary conditions need to be specified for the extreme nodes.For instance, if equation ( 4) is assumed to describe the heat flow through a single-layer homogeneous wall with Newton's law of cooling on both sides, then: where: T l -outside temperature on the left side of the wall [K], T r -outside temperature on the right side of the wall [K].
It soon became evident that the FDM is well-suited for solving relatively simple differential equations (examples shown in Fig. 2).However, the method proved to be inadequate for both solid and fluid mechanics applications.Its primary limitations included challenges in adapting the mesh to complex geometries, issues related to calculation stability, difficulties in handling boundary conditions, often resulting in a loss of accuracy, the necessity for uniform grid spacing, and the associated requirement for a large number of nodes.These challenges prompted the search for alternative numerical modelling methods that could better meet the specific needs of various physics departments.In modern Computational Fluid Dynamics, the FDM is rarely employed.Nevertheless, it holds a significant place in the history of numerical modelling, representing an important milestone in the development of simulation techniques.

Finite Volume Method
The Finite Volume Method (FVM) (McdonaLd 1971, MaccorMack, PauLLay 1972, MoukaLLed et al. 2016) is designed for modelling fluid flows within a discretized computational domain (Ω) using the Eulerian approach.This method relies on a mesh-based approach, where a fundamental role is played by the concept of the Finite Volume or Control Volume (CV).A Control Volume is a volume V ∈ Ω enclosed by a closed surface S. The orientation of this surface is defined by the outward-directed vector  ⃗ (as depicted in Fig. 3).In numerical models, control volumes represent the cells of the computational grid.The FVM is based on two types of balance equations: surface and volumetric balance equations.where: Θ -any physical quantity, e.g.mass, momentum or energy,  ⃗ -fluid velocity,  ⃗ -surface directional vector.
An example of a quantity that requires surface balancing is heat flux.The resulting value is determined by summing up the flows entering (with a negative sign) and exiting (with a positive sign) through the surface S of a Control Volume (CV).
The volumetric balance equation describes the ability to change the value of Θ within the volume V: where: s Θ -source of the Θ quantity.
An example of a quantity that requires volumetric balancing is the heat generated in a combustion process.In this case, the source term has a positive sign, but in general, a source term may increase or decrease the value of the balance quantity in a control volume, depending on its physical meaning.After applying the surface and volume balances to the entire computational mesh within the domain Ω, the following set of balance equations is obtained (MoukaLLed et al. 2016): where: e -unitary energy (a sum of internal and kinetic energies) Equations (9-11) represent the Mass Balance Equation (MaBE), the Momentum Balance Equation (MoBE), and the Energy Balance Equation (EBE), respectively.These equations can be expressed in what is known as vector form: where the similarities between specific balance equations become much more apparent (please note that the terms containing pressure are now positioned on the left side).
The first part of equation ( 12) represents the time-dependent term, which describes the rate of change of the balanced quantity with respect to time.For steady-state flows this term is equal to zero.The second term in equation ( 12) is the so-called convection term, which describes the exchange of the balanced quantity between the current control volume (CV) and neighbouring cells.The components of this term are reversible, meaning they can be converted.An example is the transfer of pressure energy during compression and release during fluid expansion.The third term in equation ( 12) is the dissipative term, which describes irreversible phenomena.This encompasses all phenomena related to viscosity, turbulence, diffusion, radiation, etc.The final part of equation ( 12) is the source term, which accounts for the possibility of changing the balance quantity within a CV, such as due to the influence of external forces.
The FVM is based on cell-centered meshes, which are typically structured or unstructured.The mesh is generated within the domain, and its boundaries coincide with the domain boundary.However, in the literature another approach may be found, the so-called Immersed Boundary Method (Peskin 1977), in which the domain is immersed in the grid and does not conform with the boundary.In the immersed boundary method, special treatment has to be taken at the boundary to incorporate the boundary conditions.The main advantage of the Immersed Boundary Method is the use of structured meshes, which are easy to generate, and the calculations performed on them are faster.The greatest difficulty of the method is the imposition of boundary conditions.
The equations set (9-11) has to be supplemented by other mathematical dependencies (closures mentioned before).Two of them are in particular critical: the laminar stress tensor, in which the key role plays the choice of the appropriate fluid model (rheology is the branch of physics that deals with this issue) and the turbulent stress tensor, which can be described in many different ways.In fact, turbulence modelling is the most fundamental issue in the fluid mechanics still not fully solved.
In Figure 4, an overview of applied turbulence models is visible.Reynolds--Averaged Navier-Stokes (RANS) is a group of turbulent flow simulation methods based on the concept of decomposing the Reynolds velocity field, in which the -based on the concept of turbulent viscosity -in this approach, turbulent stresses are proportional to the deformation rate, analogous to viscous stresses (the turbulent stress tensor serves as a form of correction to the viscous stress tensor); -not based on the turbulent viscosity concept -in this approach, the components of the turbulent stress tensor are calculated directly from additional equations, algebraic (ASM), differential (RSM) or evolution.
The evolution equations are added to the system of equations ( 12) (they must have the same structure) and are solved using the same numerical techniques.The evolution equations are not used to determine a new balance quantity but a variable necessary to solve the mass, momentum, or energy balance equation.Typically, evolution equations are needed to determine the turbulent stress tensor.Equations from the RSM group are not included in the set (12) and are solved independently of them.
Direct Numerical Simulation (DNS), as proposed by Orszag in 1970, is a method for simulating turbulent flows.It involves the direct solution of the Navier-Stokes equations without any simplifications, which means that all turbulence scales are considered in the calculations.Calculations performed with the DNS method enable the correct reconstruction of the dynamics of all spatial and temporal scales of turbulence, which, unlike the RANS methods, are not modelled, but are the result of a numerical solution.The DNS method requires enormous computing power and for this reason is still not applicable in practical issues.
Large Eddy Simulation (LES), as introduced by Smagorinsky in 1963, is a method for simulating turbulent flows.It involves the separation of vortices into two scales through the use of specific filters.The large-scale fluid movement is calculated based on the Navier-Stokes equations, while the small-scale structures (smaller than the width of the filter used), in which kinetic energy is dissipated by the action of viscous forces, are modelled.The smaller the filter width, the closer the LES method becomes to DNS.
Detached Eddy Simulation (DES), developed by Spalart in 1997, is a method for simulating turbulent flows that combines elements of both Large Eddy Simulation and Reynolds-Averaged Navier-Stokes.The RANS approach is applied in the boundary layers (where large vortices are absent), while the LES approach is employed in the core of the flow.Achieving a seamless transition of fields between the LES and RANS regions in DES necessitates the use of high-quality numerical grids.
It is currently believed that the most versatile model suitable for engineering modelling of typical fluid flows is the k -ω SST (Menter 1993) model from the RANS group.This model accurately reflects phenomena both in the boundary layer and far away from it.In some areas, other RANS models are recommended, such as the Spalart-Allmaras model, specifically designed for aerospace applications.LES or DES are used in cases where various turbulence scales occur, such as large coherent structures in vortex flows, which cannot be adequately modeled with the RANS approach.
The equations set (9-11) is valid for single-component fluids (soBieski, grygo 2019).However, with appropriate modifications and extensions, it can also be used for a wide range of multiphase flows.In Figure 5, the one-fluid approach means that the occurrence of certain phenomena may be predicted from scalar or vector fields obtained for a single-phase flow.An example of this is the barotropic cavitation model, in which the cavitation phenomenon is predicted based on the pressure field.Another way to introduce a second phase is by incorporating some effects of its presence.Good examples of this approach are the Porous Media Model or the Porous Jump Model, in which flow resistance is modelled by a source term in the momentum balance equation within a chosen volume or on a surface, respectively.
In the homogeneous approach, the mixture is considered as a substance composed of any number of phases and is modelled using only one set of balance equations (Niedźwiedzka et al. 2016).Mass or volume fractions of each individual phase are calculated using additional balance equations, which increases the number of variables in the vectors seen in equation ( 12).Depending on further details, the phase interface may or may not be tracked.The first case, known in the literature as the Volume of Fluid model, is used in simulations of fluid flows with a free surface or in flows involving two non-mixing liquids.Each of the components can be a single fluid or a homogeneous mixture of any number of components.The homogeneous approach is typical for modelling selected mixing phenomena, flows involving evaporation and condensation, or flows with chemical reactions.
In the non-homogeneous approach, each component of the mixture has its own set of balance equations (soBieski 2009).Phase coupling occurs through pressure and interfacial coefficients related to mass, momentum, and energy exchange.The description of this interaction primarily depends on the type of flow, whether it is liquid-liquid, liquid-gas, liquid-solid, or gas-solid.This approach is primarily intended for flows with one background phase (liquid or gas) and one or more dispersed phases (bubbles, droplets, or solid particles).The fluid phase is described using the Eulerian approach, while the dispersed phase can be treated either in the Eulerian manner (as in the Eulerian Multiphase Model) or in the Lagrangian manner (as in the Discrete Phase Model).In the latter case, appropriate movement equations must be added to the model.In practice, this task is accomplished by coupling the FVM (or another method designed for fluids), responsible for managing the continuous phase, with the Discrete Element Method, utilized in this context for modelling the dispersed phase.The Discrete Element Method is likewise described in the article.
The FVM is highly popular, and examples of its application can be found in every field of modern science.Some examples from mechanical engineering are shown in Figure 6.The primary limitation of this method remains the insufficient computing power of currently available computers.

Lattice Boltzmann Method
The Lattice Boltzmann Method (LBM) (BHatnagar et al. 1954, sukoP, tHorne 2006) is used to model the behaviour of a lattice gas in a state of thermodynamic imbalance, employing a Lagrangian approach.The foundational concept in the LBM is the phase space, where each dimension corresponds to a distinct physical quantity.In classical mechanics, the phase space typically encompasses position, velocity or momentum, and time (as shown in Fig. 7 where: f -distribution function of gas density in the phase space,  ⃗ -position vector (here in Cartesian coordinates), t -time.
The function f( ⃗ ,  ⃗ , t) represents the probability of finding a given particle at a given location in the space  ⃗ and time t with a given velocity  ⃗ .If there are no particle collisions, then during the movement in the phase space, all those particles that, at time t + dt, are in the phase space element described by the coordinates  ⃗ +  ⃗ and  ⃗ +  ⃗ , were at time t within the phase space element with the coordinates  ⃗ and  ⃗ where:
If particle collisions occur, meaning that particles move not only away from a given point in space but also towards it, the derivative does not equal zero: The simplest model describing the collision term is the BGK linear model (BHatnagar et al. 1954): where: τ -the so-called relaxation time.This parameter determines how quickly the f function may obtain the equilibrium value, which depends on the assumed viscosity of the fluid.In turn, f eq means the distribution function of gas density in the equilibrium state: where: If the values of the distribution function f( ⃗ ,  ⃗ , t) in each i th direction are known, different macroscopic parameters may be calculated.For example the density and velocity of lattice gas may be determined with the following formulas (sukoP, tHorne 2006): where: n i -the number of spatial directions in the model (equal 8 in the D2Q9 model; 0 is referred to the current lattice node), e i -the unit vector of the i th velocity direction (depends on the chosen lattice model).
A notable characteristic of the LBM is that the computational domain must be represented using a binary geometry.This kind of geometry can be directly prepared, for example, using random techniques or indirectly converted from another type of geometry.A conversion scheme from a vector geometry to a binary geometry is depicted in Figure 8.In this process, the original vector geometry (a) is overlaid by a lattice of nodes (b), and each node is assigned a value of 1 (solid node) or 0 (fluid node) based on its location.This method of defining the computational domain enables the relatively easy preparation of even highly complex geometries.Indeed, it is the primary advantage of this method and indicates its fundamental areas of application.However, the major drawback is that the quality of the resulting numerical model strongly depends on the assumed lattice resolution.Increasing the node density improves the model's quality but also increases the computational power requirements and the time needed to obtain results with the same level of convergence.Another feature of the method is that the movement of the lattice gas can only occur along strictly defined directions, and the probability of selecting a given direction is determined by appropriate weights.The weight value is highest for the current node, meaning that a significant number of molecules always stay in the same place.The most popular lattice models are D2Q9, D3Q15 (limited range of stability, low computational effort), D3Q19 (a compromise between stability and computational effort), and D3Q27 (good stability, high computational effort).D2 or D3 represents the number of space dimensions, and Q9, Q15, Q19, or Q27 represents the number of possible directions of gas movement.The lattice models described here, along with the weight values of individual directions, are presented in Figure 9 (sukoP, tHorne 2006(sukoP, tHorne , wagner 2008)).
In Figure 10, the main steps of the LBM algorithm are presented.In stages a and b, formulas (20-21) and ( 19) are used in a discrete form, respectively.In step c, collisions between particles and between particles and walls are included.The collision process involves swapping values of the appropriate components of the ( ⃗,  ⃗, ) function.The specific details depend on the assumed collision model, whether it includes slippage or not.On this stage, periodic conditions may also be implemented.The symbol f * represents the f function after the collision stage but before the streaming procedure.In step d, the lattice gas is moved (streamed) along individual directions of the assumed lattice model, and this step includes the previously mentioned weights.
The model described in this chapter is applicable to single-phase laminar fluid flow.However, the LBM has now reached an advanced stage, enabling the modelling of a wide range of fluid flows, including the same fluid models, turbulence models, and multiphase flow concepts as in the FVM.The areas of application are also similar (e.g., Fig. 11a), but the LBM is particularly useful for modelling flows in highly complex geometries (e.g., Fig. 11b).Interesting examples of applications can be found in medicine (krause 2010, aBas 2016), chemical engineering (Xie et al. 2018), agricultural engineering (korn, HerLitZius The main drawback of the LBM is its significant computational power requirement, which limits its practical application in many cases.To expand its applicability, much attention is devoted to parallelization techniques or the use of graphical processing units (GPUs) instead of CPUs in computations.The resultant force can be expressed as the sum of two types of actions: where: In turn, contact forces consist of two components: where:

26, 2023
Numerical Methods in Fluid Mechanics -An Overview

205
In the case of the i th circular or spherical body in contact with other similar bodies, equations ( 22) and ( 23 where: F n,ji -normal force acting on the i th body in contact with the j th body [F], F t,ji -tangential force acting on the i th body in contact with the j th body [F], r i -radius of the i th body [m], f ij = f ji -distance between the direction of the resultant normal force acting between the bodies i and j, and the centre of rotation of the i th body [m]. If the contacting circular or spherical bodies are non-deformable, the contact zone is limited to a point, and the distance f ij = f ji is equal to zero (Fig. 12a).The distance between the centres of the particles i and j is equal to r i + r j .This type of contact in the literature is called hard contact.If the bodies are deformable (Fig. 12b), then instead of a point contact, a line (in 2D) or a surface (in 3D) contact will be created.In this case, the resulting point of action of the normal force depends on the distribution of normal stresses on that line or surface.Consequently, f ij = f ji > 0, and the normal force will generate an additional moment of force.The particles are deformable, therefore the distance between their centres is smaller than the sum of the radii r i and r j .In this way, the degree of particle overlap becomes a measure of their deformation.The contact of deformable particles in the literature is called soft contact.In summary, in the DEM, two main tasks can be distinguished: -searching for all contacts between objects; -calculating forces and moments of forces at all points of contact based on adopted contact models.
Contact models are based on different material models, and additional factors may also be introduced, such as forces resulting from the existence of a liquid layer on the particle surfaces.The most popular material models include: -elastic model -a model in which the strain is reversible and does not depend on the relative velocity between bodies.A distinction is made between a linear elastic model, according to Hooke's law, and non-linear models, e.g. the Hertz model (HertZ 1881): , (  ) =     (28) where: δ n -absolute deformation at the contact point in the normal direction [m], k n -modulus of elasticity in the normal direction [N/m α ], α -model exponent [-].
For mentioned Hooke and Hertz models is equal to 1 and 3/2, respectively; -viscoelastic model -a model in which the deformation of a body is assumed to be a partially irreversible process (due to internal damping).Linear and non-linear viscoelastic models are distinguished: where: ̇ -speed of absolute deformation at the contact point [m/s], β -model exponent [-]; -elasto-plastic model -a model in which it is assumed that in the first phase of loading, a body behaves elastically, and after a certain limit (yield point) is exceeded, the behaviour of the body becomes typical for a plastic material.This model can be described by various combinations of linear and nonlinear elastic or plastic models.
When determining the tangential force, two scenarios are usually considered.If the tangential force is relatively small, the tangential deformation is calculated.In the simplest approach, similar formulas are used as when calculating the normal force.If the tangential force exceeds the static friction force, sliding occurs with a constant dynamic friction force.The displacement in the normal direction is only dependent on material properties, not time.In contrast, if the colliding bodies are in a state of sliding, the relative displacement in the tangential direction depends on time.
Technical Sciences

26, 2023
Numerical Methods in Fluid Mechanics -An Overview

207
In Figure 13, the main steps of the DEM algorithm are presented.To simplify the diagram, only forces are shown, but moments of forces are calculated in an analogous way.
A distinctive feature of the DEM is that the computational domain has no explicit boundaries.This means that the domain size depends on the current location of the objects involved in the simulation.Alternatively, the maximum range of coordinates for a simulation may be defined.A typical area of application for the DEM is the modelling of loose materials' behaviour in various processes, including transportation, loading and unloading, tank filling and emptying, collection, spreading, separation, and many others.An example from this field is shown in Figure 14a.
The DEM also bridges the gap between granular mechanics and classical solid mechanics.By utilizing carefully selected contact models, it can simulate phenomena such as cracking, delamination, crushing, grinding, pelleting, and other processes in which the internal cohesion of materials is disrupted or formed.Additionally, the DEM has significant applications in modelling the mechanics of elastic and highly deformable bodies, such as ropes, membranes, foils, nets, fibres, rubber, silicone, and more.
In the context of CFD, the most interesting application of the DEM is its combination with FVM or LBM for modelling fluid-solid systems.Examples in this field include fluidized beds (FuLcHini et al. 2019, Liu, ZHao 2021), cyclones (cHu 2010), soil erosion (Zou et al. 2020), and more.The method can also be utilized to generate the geometry of granular porous media, which is subsequently employed, for instance, in the analysis of fluid flows within the pore space (soBieski et al. 2018).In the first of the areas mentioned here, the coupling of DEM with FVM or DEM with LBM is performed iteratively within a single computational loop.In the second case, the connection is sequential (Fig. 14b).The DEM can also be used independently, without combining it with other numerical methods, to model the behaviour of semi-liquid materials such as pastes (PLatZer, FiMBinger 2021), concretes (gao et al. 2023), 3D printer filaments (kHan, koç 2022), and so on.

Smoothed Particle Hydrodynamics
Smoothed Particle Hydrodynamics (SPH) (Liu, Liu 2003, Li, Liu 2008, Lind et al. 2020) is a mesh-less Lagrangian method for analysing the behaviour of solids or fluids in which continuous matter is represented by a set of particles.
Technical Sciences

26, 2023
Numerical Methods in Fluid Mechanics -An Overview 209 Each particle represents a small volume of the matter and contains individual parameters such as mass, density, pressure, velocity, etc.The starting point for SPH comes from the observation that a function operating on a domain may be expressed in integral form as follows: ( ⃗) = ∫ ( ⃗′) ⃗′ Ω (30) where:  ⃗ -position vector of the point for which the value of the function is sought,  ⃗′ -position vector of a point where the value of the function is known.
In equation ( 30), if the Ω space is reduced to a single point, then the value of the function f( ⃗ ) will be equal to the value of f(  ⃗′ ) at that point.Assuming that Ω space is continuous, a value will always be obtained for each  ⃗ location (see Fig. 15a).However, if Ω is a discrete space, two scenarios are possible: -if the  ⃗ location coincides with one of the  ⃗′ locations, then the correct value of the f( ⃗ ) function can be obtained; -in the other case, the f( ⃗ ) function will be zero (see Fig. 15b).This can be expressed in the following form: The main question is how to obtain a value for any location in a discrete system.To address this issue, the following approximation can be applied (as shown in Fig. 15c) where: ( ⃗ −  ⃗ ′ , ℎ) -a weighting function (also known as a smoothing function or kernel function) that determines the degree of influence of a given location  ⃗′ on the approximation result at point  ⃗ .h is referred to as the smoothing length, which is typically constant for all particles and is responsible for defining the extent of the approximation area.The smoothing length is often multiplied by a factor κ, typically falling in the range between 1 and 3.In the literature, the space within the circle or sphere κh is referred to as the support domain.
The smoothing function must satisfy several important conditions, the most significant of which are as follows: -normalization condition (eq.33): This condition ensures that the influence of discretization points within the support domain sums up to 1.In other words, each discretization point contributes a percentage share to the calculation of the approximated value; -Dirac function property (eq.34): This property dictates that when the support domain is reduced to a single point, the smoothing function behaves like the Dirac delta function; -compactness conditions (eq.35): These conditions ensure that the value of the smoothing function is zero outside the support domain.
The smoothing function must meet a number of conditions, the most important of which are: -normalization condition (eq. 33), what means that the influence of the discretization points in the support domain sums up to 1 (each of the discretization points has its percentage share in the calculation of the approximated value); -Dirac function property (eq. 34), what means that in the case of reducing the size of the support domain to a point, the smoothing function becomes the Dirac delta; -compact conditions (eq. 35), what means that the value of the smoothing function outside the support domain is zero: SPH is among others dedicated to the analysis of fluid flows.The fundamental equations of fluid dynamics can be expressed in the following form (Liu, Liu 2003, Li, Liu 2008): After applying the previously described method, the mass balance equation (39) will be discretized as follows (Liu, Liu 2003, gesteira et al. 2010): where: i -particle at location  ⃗ , j -particles at locations  ⃗′ , N -the number of particles belonging to the support domain of the i-th particle, W ij -a shorter notation for the smoothing function.
The discrete form of the momentum balance equation ( 40) can be expressed as follows (Liu andLiu 2003, gesteira et al. 2010 .In this case, the shear tensor is calculated based on physical viscosity.Another popular approach is to use artificial viscosity, where the ( ⃡, ) function is represented by the artificial viscosity, typically denoted as Π ij .The advantages of using artificial viscosity are to ensure the stability of calculations when particles move apart, ensure the correct direction of momentum transfer, and facilitate implementation.
To calculate the term f( p, ρ), the equation of state must be applied.For gases, the ideal gas law may be used (arcHaMBeau 2013, Liu, Liu 2003).In the literature, two forms of this equation may be found: Wojciech Sobieski, Božidar Šarler The time step is limited by conditions that take into account the intensity of mass forces and the intensity of viscous dissipation.
In Figure 17, the main SPH algorithm is presented.It is worth noting that nested loops are executed during a single time step: one loop iterates over all fluid particles, and another loop iterates over all particles belonging to the support domain of every fluid particle.When there are walls within the calculation domain, appropriate boundary conditions must be introduced.The three most popular variants are as follows: 1) Adding an artificial repulsive force to prevent liquid particles from penetrating the wall; 2) Introducing so-called ghost particles outside the wall with the same parameters as fluid particles but with the opposite direction of velocity.This results in velocity components on the wall being equal to zero; 3) Treating wall particles in the same way as fluid particles, requiring the definition of two layers of non-moving wall particles.Conditions 1 and 2 are commonly used.The third condition is sensitive to the initial distance between the fluid and the wall.When the initial distance is too great, the fluid may move towards the boundary, and when it is too small, the fluid molecule will experience suction as it moves away from the boundary.
The SPH method is primarily employed for modelling systems in which matter undergoes significant deformations.It finds applications in fracture mechanics (douiLLet-greLLier et al. 2016)

Summary
Over the past few decades, numerous specialized numerical modelling methods have been developed.These techniques have been refined to the extent that they now constitute fundamental research tools in almost every field of science.Initially, there was a clear distinction between methods designed for modelling the behaviour of solids, fluids, or dispersed matter.Often, choosing one approach made it impossible to incorporate other options.However, the progress of science and current trends are expanding the possibilities of simultaneously using multiple methods and their variants.In particular, simulations that combine two or more numerical methods are becoming more prevalent.For example, the FVM or LBM can be combined with the DEM (e.g., for modelling fluidization, sedimentation, etc.), or the FVM can be coupled with the Finite Element Method (not described in this paper, e.g., for modelling fluid flow in flexible conduits).This development allows researchers to address complex problems that require a more comprehensive and integrated approach.
It is worth emphasizing that the range of possible variants in a given numerical model of physical phenomena is vast.The palette of techniques and variations is also highly diverse, particularly evident in the case of multiphase systems or domains with moving boundaries.This diversity often leads to communication difficulties between researchers involved in experimental studies and numerical modelling.Even a minor change in assumptions, concepts, or research objectives may necessitate the use of entirely different mathematical frameworks or computational techniques.Since researchers specializing in simulation studies typically focus on specific problem classes, such changes can significantly delay the entire task due to the need to familiarize themselves with new theories, relevant literature, and the practical application of newly acquired concepts in selected software.Moreover, the time required for conducting a single simulation can extend to days or weeks.Considering that the development of new models often involves testing multiple variants (e.g., during the search for new closures), obtaining initial, qualitatively correct results may take months or even years.This situation often leads to misunderstandings, as experimenters or engineers usually expect rapid results.The necessity for laborious preparation when approaching modelling new problem classes is one of the factors limiting the widespread application of numerical modelling at the engineering level.
Another critical factor is the current inadequacy and insufficient computational power of today's computers.It can be speculated that a new chapter in numerical research will begin with the advent of long-awaited quantum computers.Furthermore, some currently employed approaches are expected to become less attractive: Reynolds-Averaged Navier-Stokes may be replaced by Direct Numerical Simulation; diffusion-convection equations may be superseded by full mass, momentum, and energy balance equations; macro-scale Porous Media Models may be substituted by micro-scale modelling of fluid flow in porous channels, and so on.It is also important to note that greater computational power will enable modelling on a significantly larger scale, such as simulating the entire flow system in a hydroelectric power plant rather than focusing solely on a single interblade channel.This development will hold paramount practical significance.

Fig. 1 .
Fig. 1.Illustration to the main idea of the FDM

Fig. 2 .
Fig. 2. Examples of using the FDM: a -resolving the diffusion equation in 2D space and modelling of the spread of smell in the layout of laboratory rooms, b -resolving the wave equation in 2D and modelling of waves propagation on a surface with random falling drops

Fig. 8 .
Fig. 8. Schematic representation of geometry conversion from a vector to a binary geometry; description in the text

Fig. 9 .
Fig. 9.The most popular lattice models in the LBM

Fig. 11 .
Fig. 11.Examples of using the LBM: a -modelling a flow around the NACA6412 airfoil, b -investigation of the hydraulic tortuosity in random generated pore structures

Fig. 13 .
Fig. 13.Main steps of the DEM algorithm

Fig. 14 .
Fig. 14.Examples of using the DEM: a -mixing particles in a rotating tank (own calculation of an example available in the YADE DEM code (The YADE code 2023)), b -predicting the flow resistance by the use of the FVM in a granular porous bed created by the DEM Source: own calculation of an example available in the YADE DEM code (2023).
Fig. 15.Determining a function value in: a -a continuous space, b -a discrete space based on Dirac function, c -a discrete space based on smoothing function

Fig. 18 .
Fig. 18.Example of using SPH to modelling movement of objects floating on the water surface near a coast Source: own calculation of an example available in the SPHYSICS_2D code (gesteira et al. 2010).
The publication was written as a result of the author's internship in the Ljubljana University in Slovenia, co-financed by the European Union under the European Social Fund (Operational Program Knowledge Education Development), carried out in the project Development Program at the University of Warmia and Mazury in Olsztyn (POWR.03.05.00-00-Z310/17).
) will take the form: ):