In February 2025, there were significant improvements to MULES in the development line of OpenFOAM (OpenFOAM-dev, packaged here) from The OpenFOAM Foundation. The links to the relevant code commits are provided below.
- Improvements to the iterative MULES algorithm to guarantee boundedness irrespective of limiter convergence. The log for that code commit was partially incorrect, but was clarified in the log of a small subsequent commit.
- A new control structure for MULES parameters.
- Optional controls based on limiter convergence for semi-implicit MULES.
- The control for sub-cycling the α-equation now uses a Function1 of the α-Courant number. This is applied both to sub-cycling in VoF interface-capturing and to sub-cycling in the multiphaseEuler solver.
Background to MULES
“MULES” is an acronym for Multi-dimensional Universal Limiter for Explicit Solution. The method is designed for equations whose solution variable, e.g. ψ, which has some lower and/or upper bound(s). Its purpose is to calculate accurate solutions for ψ which are also bounded, i.e. do not fall outside of the bound(s) of ψ. The method was developed primarily for interface-capturing using the volume of fluid (VoF) method, in order to maintain phase fraction α between its bounds of 0 and 1. It is also used within the Euler-Euler multiphase solver (multiphaseEuler) and is an option for general scalar transport (in the scalarTransport function object).
MULES works by limiting the size of the fluxes of the variable ψ (in units of ψ per second) through cell faces, calculated by interpolating ψ from cell centres to the faces. Upwind interpolation guarantees a bounded solution, both when fluxes are treated implicitly or explicitly, although a time step limit for Courant number (Co) = 1 is also required for boundedness of an explicit solution. Upwind, however, is generally too inaccurate to be practically useful, e.g. for calculating phase fraction α in the VoF method.
The calculation of fluxes therefore needs an interpolation scheme that delivers better accuracy, but such a scheme cannot guarantee boundedness. Instead, a limiter λ is introduced in the range 0 ≤ λ ≤ 1 which is applied as a multiplier to “correction” fluxes that include the accurate, unbounded scheme. Calculations split the fluxes into a sum of (upwind fluxes + limited correction fluxes) to maintain boundedness. In the worst case, when λ = 0 at a given face, the scheme reduces to upwind for least accuracy. In the best case, λ = 1 and the full correction will be applied, corresponding to the chosen accurate interpolation scheme, e.g. the van Leer scheme.
Original MULES (2006)
MULES provides the method to calculate λ at each face. More specifically, in cases with lower and upper bounds, it calculates two limiters λ+ and λ− , corresponding respectively to: positive fluxes out of a cell, which tend to decrease ψ towards its lower bound; and, negative fluxes which tend to increase ψ towards its upper bound. MULES was created by Henry Weller in 2006 and released in OpenFOAM v1.4 in 2007, providing a “significant improvement to robustness and efficiency of two-phase interface-tracking solvers”. MULES delivered on this promise, but the implementation in 2006 only provided an explicit solution. Generally it operated under a maximum Co = 0.5 with additional sub-cycling that effectively reduced the maximum Co for the α-equation to below 0.25. This tight restriction on time step made the original version of MULES prohibitively costly for large scale problems.
Semi-Implicit MULES (2014)
Weller overcame the time step limitation for MULES with a new “semi-implicit” option released in OpenFOAM v2.3.0 in 2014. This enabled cases to run at Co > 1, including many example cases with Co = 5. For steady flows such as some ship hydrodynamics, example cases use Co = 10 and even Co = 25. The new semi-implicit option exploited the splitting the fluxes into the linear combination of upwind fluxes and limited correction fluxes described above. Rather than solve an explicit solution on the entire set of fluxes, it performs an implicit solution using the upwind fluxes first (which is guaranteed to be bounded). The limited fluxes are then applied in an explicit correction step. There is no Co restriction on the initial implicit solution. The limiting also ensures no Co restriction on the explicit correction, but as Co increases, so does the limiting (i.e. decreasing λ+ and λ−), causing accuracy to reduce.
MULES in 2025
The calculation of λ+ and λ− in MULES is iterative. Until 2025, the code initalised λ+ and λ− to 1 and converged downwards towards values that would produce a solution within an acceptably small level of unboundedness, e.g. 0.1%. The recent change initialises λ+ and λ− to 0 so that successive iterations increase the values towards their converged level. The method now guarantees boundedness because the limiters cannot exceed their optimal values which produce a bounded solution.
MULES also includes a limiter convergence test that enables users to specify a tolerance for the limiter calculation. A tolerance control can be applied to stop the iterations automatically, rather than executing a fixed number of iterations whether they are needed or not.
The calculation of λ+ and λ− occurs first within each cell based on extrema (bounds) that are evaluated from adjacent cells. In certain circumstances, the values in a region of flow can contract within some local extrema, e.g. 0.2 to 0.8, which fall well inside the global/physical bounds, e.g. of 0 to 1. Solutions can become “stuck”, preventing them to return to the actual bounds. In order to overcome this problem, further controls are available to relax the local extrema to enable a solution to reach its global bounds again.
The control structure of MULES has been updated to include a MULES sub-dictionary to simplify these parameter names. In the example below, a limiter tolerance of 0.01 is used which is usually sufficient for α. The number of iterations, nIter, is then set as a high value, simply to stop the calculation in the event of non-convergence. The extremaCoeff setting of 0.1 relaxes the local extrema on α by 10% of the difference in global bounds. A value of 0 implies no relaxation, i.e. uses the local extrema; a value of 1 causes MULES to use the global/physical bounds.
solvers
{
"alpha.water.*"
{
...
MULES
{
nIter 10;
tolerance 1e-2;
extremaCoeff 0.1;
}
...
}
...
}