An OpenFOAM solver for computing suspended particles in water currents

The science of Hydroinformatics includes numerical modelling of water flow in lakes, rivers, canals and ducts (Maduka & Li 2021; Olsen 2022; Paschmann et al. 2022). Three-dimensional water flow fields are often found by solving the Navier–Stokes equations using a turbulence model in a science called computational fluid dynamics (CFD). The set of equations are complex to solve for a general geometry. The most used CFD programmes may have a source code with millions of lines that have taken many decades to make. There are several issues related to the complexity of codes with a large number of lines. Industry average experience is about 1–25 bugs pr. 1,000 lines of code (McConnell 2004). Also, the transfer of coding knowledge to new staff is very time-consuming. The risk of inexperienced staff introducing new bugs in the programme is considerable. Companies making CFD software are very dependent on specialised personnel with coding experience, and problems arise when key staff retire or move to other companies. Also, it is difficult for researchers outside the company to get access to the code, both for legal reasons and for the large amount of knowledge that needs to be transferred.

A solution to these problems is to use a high-level object-oriented computing language with open source access. The current article describes a solver made in OpenFOAM, which is based on the object-oriented C ++ programming language. OpenFOAM is a library of different solvers that the users can extend by adding equations and thereby creating new solvers for specific flow problems. In the current study, it was possible to create a solver for computing particle transport in water flow with only around 70 lines of new code. The rest of the code is in libraries that are written, maintained and tested by very experienced programmers. The new solver takes the fall/rise velocity of the particles into account by using a drift-flux approach. It is therefore called sediDriftFoam, as it is expected that it will mostly be used for computing sediment transport. However, it can also be used for computing transport of plastic particles. The code for the new solver is explained in detail in the current article, together with the equations that are solved. The source code for sediDriftFoam can be freely downloaded from the link provided in the Supplementary Material. Example input files are also given. The solver can thereby easily be used and further modified.

Earlier work using this approach includes Sattar et al. (2017), who programmed a sediment solver in OpenFOAM with a moving grid, and computed bed movements in a channel bend. Other examples in the field of sediment transport are SediFoam (Sun & Xiao 2016) and sedFoam (Chauchat et al. 2017). SediFoam models particles using a Lagrangian approach, where a large number of particles are released in the flow domain and their paths are tracked individually. The movements of the particles are a function of their acceleration, computed from Newton's second law of motion. The forces on the spherical particles are due to gravity, lift, drag, collision and lubrication. Several formulas with empirical coefficients were used for each force. The Lagrangian approach to modelling particles is popular within the field of fluid mechanics (De Marchis & Milici 2016; Schwarzmeier et al. 2023).

sedFoam (Chauchat et al. 2017) is an Eulerian two-phase model, where fractions or concentrations of particles are computed for each cell in the grid. The model solves momentum equations both for the water and for the particles. Fluid stresses include both Reynolds stresses and granular stresses, including fluid-particle interactions on the grain scale. The model is fairly complex, with 47 equations describing momentum, continuity and different turbulence modelling options.

Both SediFoam and sedFoam focus on fundamental research on fluid mechanics and particle movements. Much simpler models are often used to solve hydraulic engineering problems, for example Telemac, Delft3D, SSIIM and Flow3D. These models assume a one-way coupling where the particles will not affect the water velocity and turbulence, except when erosion and deposition changes the geometry. An Eulerian approach is used where the particles are computed as a passive concentration. Well-known formulas for particle pick-up rates are used as boundary conditions at the channel bed. These computer programmes have been used to model most sediment transport phenomena, including local scour (Baranya et al. 2014), turbidity currents (Basani et al. 2014), reservoir sedimentation (Mool et al. 2017), reservoir flushing (Harb et al. 2014), formation of braided rivers (Olsen 2021) and bedform modelling (Goll & Kopmann 2012).

The current work describes a particle model with a one-way Eulerian approach, implemented in the OpenFOAM library. This simpler solver uses a drift-flux approach to take the rise or fall velocity into account when computing the dispersion of the particles. It also includes particle pick-up rates for sediment flow over a fixed bed. The model is therefore called sediDriftFoam. The idea behind the drift-flux model is to superpose two velocity vector fields to produce the particle movements: one water velocity field from the solution of the Reynolds-averaged Navier–Stokes equations, and another field from the particle fall velocities. The algorithms are the same as used by Olsen & Skoglund (1994), but the current article provides much more details about the computational programme and its integration with OpenFOAM.

The new solver was tested on two cases: Vertical distribution of positively buoyant plastic particles in a straight channel, and sediment deposition in a sand trap. The reason for choosing these two cases was the availability of input and verification data, and their relevance for current research in water engineering. Problems caused by plastic particles have received considerable attention over the last years. Micro and macro plastics are transported from rivers to the ocean, where they can be harmful to marine life. A number of studies have been conducted to model plastic particle dispersion in the ocean by numerical models. Most of these studies used a 2D depth-averaged model with either an Eulerian (Defontaine et al. 2020) or a Lagrangian (Liubartseva et al. 2016) approach. A more theoretical study on how plastic particles with different densities and sizes move in a field of turbulence was provided by Shamskhany & Karimpour (2022). They used a Lagrangian approach with Large Eddy Simulation (LES) to resolve the turbulence and its effect on the particles in detail. Shamskhany & Karimpour (2022) used an OpenFOAM solver.

None of these studies tested the numerical models against measured plastic concentrations. This is done in the current study. Such comparisons are important to verify if the numerical model is correct and reliable.

Another hydraulic engineering problem involving particle movement is how silt and sand particles affect irrigation and hydropower projects. It is estimated that around 1% of the total global reservoir water volume is lost annually due to sedimentation (Mahmood 1987). Hydropower projects in sediment-carrying rivers need to be designed to minimise the damages caused by the particles. A sand trap or a desilting basin is then often used to reduce the sediment concentration in the water before it causes wear on turbines or deposits in the channel system (Ivarsson et al. 2021). A typical design question is how large the sand traps need to be, to achieve sufficient particle removal. Not only the length and widths are important, but also the shape of expansion regions and the turbulence at the inflow section (Paschmann et al. 2017). Physical model studies have often been used to guide the design of the desilting basins. However, the scaling of the sediment particles is often challenging (Iyer & James 2022). For example, cohesive forces may occur when scaling down fine sand particles, causing practical problems in the laboratory model. A solution to this problem is to use a 3D numerical model that computes the water flow, turbulence characteristics and dispersion of the particles, as described in the current study.

The numerical model is based on an Eulerian approach, where the particles are computed as concentrations in a three-dimensional grid. The Reynolds-averaged Navier–Stokes equations are first solved on the same grid, together with a turbulence model. The standard k-ε model (Launder & Spalding 1974) was used in the current study, where k is the turbulent kinetic energy and ε is the dissipation rate of k. OpenFOAM provides a library where these equations can be solved by calling a function with four lines of code:

• #include "uEqn.H"

• #include "pEqn.H"

• laminarTransport.correct();

• turbulence- > correct();

The lines include standard functions taken from the simpleFoam solver. The two first lines compute the water velocity and pressure, while the last function computes the turbulence.

The velocities, U_j (j = 1, 2 and 3 to represent longitudinal velocity along x, lateral velocity along y and vertical velocity along z), and the turbulent diffusion Γ were obtained from the previously called Navier–Stokes solver from simpleFoam. The third term on the left side of Equation (1) is the drift-flux term. The fall or rise velocity of the particles is called w, and z is the coordinate in the vertical direction. F_R is a sediment pick-up rate given by van Rijn (1984). In the current numerical implementation, the velocity field from the Reynolds-averaged Navier–Stokes equations is altered by adding the fall or rise velocity. The resulting convection-diffusion equation is solved by a standard OpenFOAM solver.

The flux field from the fall or rise velocity is computed by interpolation from the cell centres to the cell surfaces. The flux field from the Reynolds-averaged Navier–Stokes equations is called phi. The flux field from the fall/rise velocity is called phi_c. This is computed with the following code:

• forAll(Umod,i) {

•  Umod[i].x() = 0.0;

•  Umod[i].y() = 0.0;

•  Umod[i].z() = fallvelocity;

•  }

• phi_c = linearInterpolate(Umod) & mesh.Sf();

The variable fallvelocity is w in Equation (1). Umod is a vector field with velocities at the centre of the cells. The function linearInterpolate interpolates the velocities from the cell centres to the cell surfaces.

A new flux function with the superposition of the two flux fields is coded as:

• phi_d = phi + phi_c

The OpenFOAM/C ++ libraries automatically wrap a loop around this line, so no explicit loop statement is needed.

The solution of Equation (1) is done with the following code:

• solve

• (

•  fvm::ddt(Conc)

•   + fvm::div(phi_d,Conc)

•   + fvm::Sp(S_implicit,Conc)

•   =

•  S_explicit

•   + fvm::laplacian(turbulence- > nut(),Conc);

• );

The OpenFOAM/C ++ syntax in the code follows Equation (1). The first term including fvm::ddt is the transient term in Equation (1), with the time derivative. The second term with fvm::div(phi_d,Conc) is the convective term, which now is the sum of the two other terms on the left side of Equation (1). The last line with the Laplacian operator is the diffusive term, which is the first term on the right side of Equation (1).

There are two additional terms in the code above: an implicit and an explicit term. These two terms incorporate two different options of taking boundary conditions at the bed and the free surface into account. If the particle densities are larger than water and the bed shear stress is higher than the critical Shields value, then some particle resuspension will occur at the bed. This is given at the last term in Equation (1), F_R. The boundary condition at the bed follows the sediment equilibrium theory by Einstein (1950), implemented in the bed concentration formula by van Rijn (1984). This formula was coded in the solver and the resulting flux included in the S_explicit variable. The formula gives a maximum pick-up rate as a function of the bed shear stress. If the concentration exceeds this limit, deposition will occur. Lower concentrations will cause all the particles to resuspend, preventing deposition on the sand trap bed. This boundary condition is suitable for a situation with sediment transport over a fixed, non-erodible bed, as opposed to an alluvial channel where erosion can occur.

• scalar bedshear = bedvelocity * 0.4 / Foam::log(30.0 * dist_to_wall/roughness);

• bedshear = 1000.0 * bedshear * bedshear;

• scalar critshear = 0.047 * 1650.0 * 9.8 * particlesize;

• scalar Tstar = (bedshear - critshear)/critshear;

• if(Tstar < 0.0) Tstar = 0.0;

• scalar Dstar = particlesize * 25296.0;

• scalar vanRijn = 0.015 * particlesize/dist_to_wall * Foam::pow(Tstar, 1.5)/Foam::pow(Dstar,0.3);

• if(bedshear > critshear) { // settling particles

•  if(vanRijn > Conc[cellbed]) {

•   S_explicit[cellbed] = -fallvelocity /

•    (2.0 * dist_to_wall) * Conc[cellbed];

•  } else {

•   S_explicit[cellbed] = -fallvelocity /

•    (2.0 * dist_to_wall) * vanRijn;

•   }

•  } // end else if settling particles

Some modifications are required at the water surface when modelling plastic particles with a positive rise velocity. The drift-flux model needs to take into account that the particles will not move above the water surface. This was done by including an implicit source term in the convection-diffusion equation (Equation (1)) for the surface cells:

• if(fallvelocity > 0) // plastic

•  S_implicit[cellbed] = fallvelocity/(2.0 * dist_to_wall);

D is the water depth, k is the turbulent kinetic energy at the free surface and f is an empirical parameter in the damping function. Celik & Rodi (1988) recommended a value of f=0.43, which was used in the current study. This epsilon boundary condition was implemented using the codedFixedValue function in OpenFOAM, following Kadia et al. (2022). The code was given in the epsilon file used as boundary condition for sediDriftFoam.

The remaining details of the solver are given in the full source code of the programme, provided in the Supplementary Material. This also includes boundary conditions for the following cases.

Three grid arrangements were used to model the case. The coarse grid had 0.1 m long and 0.023–0.07 m wide cells, with 20 cells uniformly distributed in the vertical direction. The total number of cells was 65,800. The finest grid had 0.075 m long cells that were 0.016–0.05 m wide. Thirty cells were uniformly distributed in the vertical direction, giving the total number of cells as 159,300. A third grid was also made, with 0.1 m long and 0.032−0.1 m wide cells. In this grid, the number of cells in the vertical direction was varied between 23 and 32, depending on the water depth. The total number of cells in this grid was 68,540. The lines in the third grid were skewed towards the bed, following the paths of the particles, to reduce false diffusion. This grid is seen in Figure 1.

The sediment particles used in the laboratory model study were divided into four groups, according to Table 1. The numerical model was applied to each group and the concentrations were merged in the final result. The fall velocities in the table were taken from Rouse (1937).

Smooth wall laws were used on the side walls, and rough wall laws with a roughness of 1 mm were used on the bed. The distances from the bed cell centres to the bed were 1–3 cm, giving y+ values between 200 and 700 for the bed. Detailed boundary conditions for this case are given in the Supplementary Material.

The computations with the sediDriftFoam solver used the first-order upwind scheme for the convective terms in all equations. The solution converged to a steady state with sufficiently low residuals (≤10⁻⁴) after 900–6,600 time steps. The computational time of one particle size for the skewed grid was 30 min on one core of a PC with an i7 processor.

Figure 6 shows the turbulent kinetic energy distribution along the centreline of the sand trap. Note the higher values at the start of the expansion region, especially close to the bed. The turbulent kinetic energy at the boundary is linearly proportional to the bed shear stress (Rodi 1980). The high shear stress at this location causes no sediment deposition here. This can also be seen in the photograph of the expansion region from the laboratory study on the web page http://pvv.org/∼nilsol/cases/doris.

The science of Hydroinformatics aims at improving computational tools for water engineering and research. Engineers today often use commercial CFD programmes to solve practical problems in hydraulics. The programmes have intuitive and friendly user interfaces that reduce the time spent on projects. Open source alternatives, like OpenFOAM, requires more training for the users. But an advantage of the open source alternative is the possibility to learn and understand all the details on how the solution algorithms work. Commercial CFD programmes do normally not give the users access to the source code. Does a programme actually solve the exact same equations as given in the Users Manual? Are there undocumented stabilisation algorithms that introduce systematic errors? In an open source programme, the user is able to see all algorithms that are coded.

Another main advantage with open source software is the ability to change the code. The current article describes a new solver that has been made from the libraries of a high-level CFD language. The relatively few lines of code make the solver easy to learn for new users. Also, it should not be too complicated to modify the solver further. OpenFOAM contains several libraries for changing the grid, so it should be possible to extend sediDriftFoam with algorithms that compute bed level changes over time as a function of sediment erosion and deposition. It would also be possible to compute particles that have different characteristics than sediments or plastics. An example is gas bubbles formed by air entrainment in hydropower plants, causing mortality for fish.

The two examples in the current article show that it is possible to use the new solver for computing plastic particle concentrations in a flume and sand deposition in a sand trap. Tests with different grid sizes show that there is some variation in the results. However, Figures 3 and 7 show that the variations in the measurements are much larger. The largest deviation between the computed and measured values is seen in Figure 3, section 9.25 m, about 0.8 m above the bed. Improved results for this section were found by Olsen & Skoglund (1994) when modifying a coefficient in the k-ε turbulence model. In the current study, the RNG (Re-Normalized Group) and Realisable k-ε turbulence models were tested, but this did not give a converged solution. A converged solution for the sand trap was also only possible with first-order upwind discretization schemes. A second-order discretization scheme is often recommended to reduce false diffusion in CFD computations. However, false diffusion is also reduced by grid refinement. The relatively small variations in the results when using different grids is an indication that false diffusion did not have a large effect on the results in the current case.

Convergence problems for the sand trap may be due to the physics of the flow field. The upstream boundary has a small cross-sectional area with relatively high velocity. A jet is formed, and it is difficult to predict its location in the sand trap. Almeland et al. (2019) found very different flow fields in a sand trap, depending on grids, discretization schemes and turbulence models. The flow field in the sand trap may not be steady, and then it is very difficult to get a converged solution. Almeland et al. (2019) even found that the second-order scheme gave the incorrect flow field for one particular sand trap case, whereas the first-order scheme gave the better result.

A new open-source CFD solver for calculation of suspended particles in water currents is presented. The programme is made in the OpenFOAM framework. It has been tested for deposition of sand in a three-dimensional sand trap and for the vertical distribution of plastic particles with a positive rise velocity in a flume. The computed results show good agreement with vertical concentration profiles obtained previously from the two laboratory studies. In the sand trap case, the particles deposit in the expansion zone and the concentration decreases towards the water surface and in the downstream direction. The rise velocity of the plastic particles creates an opposite trend with the higher concentrations located near the free surface. The developed code is relatively simple and can therefore be applied and modified by a wide range of users, including students and researchers working on sediment transport cases.

The main limitation of the model is its one-way coupling between the particles and the water. The water will affect the particles, but the particles will not affect the water flow field. This assumption is only valid for low particle concentrations, including the examples shown in the current article.

All relevant data are available from an online repository or repositories: http://pvv.org/~nilsol/sediDriftFoam.

The authors declare there is no conflict.