Lagrangian Ordinary Differential Equation Integrators in OpenFOAM 5.x:


Recently, as I was writing my thesis, it became necessary to dig a little deeper into how OpenFOAM integrates the Lagrangian particle equations. The particles in my work represent coal and we use the ReactingMultiphaseParcel class, or a modified version of it anyways, to model them. Three Lagrangian ordinary differential equations (ODEs) are used to evolve the particles. There is one equation each to describe the evolution of:

  • the particle velocity
  • the particle temperature/enthalpy
  • the particle mass

For simplicity we will just consider the enthalpy equation as an example. The particle enthalpy equation is
$$ cp \cdot \rho \frac{dT}{dt} = S. $$
Where \( cp \), \( \rho \) and \( T \) are the specific heat capacity, density and temperature of the particle respectively. \( S \) is a general source term, that is the real magic happens and it is probably beyond the scope of this article to get into, if you really want to get into it then read my thesis. It suffices to say that in general this source term is a non-linear function of the the temperature itself, as are \( cp \) and \( \rho \). This is not an insurmountable issue (you would just need to use an iterative solver) but the complexity of these simulations needs to be limited at some point and so OpenFOAM will treat these coefficients/sources as constants over a small period of time. It calculates them at the beginning of every main solver time step but for the ode integration itself (which takes place over that time step) they are considered frozen constants. This is fantastic news for solving these things since it makes the equation non-linear, this same approach allows all three of the particle equations to be described by a linear first order ode like
$$ \frac{d \phi}{dt} = -\beta \phi + \alpha \tag{1} $$
Where \( \beta \) and \( \alpha \) are constants and \( \phi \) is the variable of interest, e.g. the particle temperature or velocity.

OpenFOAM comes with two solver classes for these Lagrangian ODEs, Analytical and Euler (look within the Lagrangian/intermediate/integrationScheme/ directory to find them in OpenFOAM version 5.x). They are actually really simple once you understand what they do and which equation they were implemented to solve (hint: it's the one above). Unfortunately when I first started looking into them I was unable to find any documentation in the code and I couldn't understand the arguments. After staring at it for a few hours I realized that it actually wasn't anything all that unusual or difficult and decided I should do a writeup that hopefully someone in the future can stumble upon and avoid wasting a few hours.

Before we jump into the code it would probably behoove anyone reading this to refresh their basic differential equation knowledge if for nothing else than just so that we can be on the same page going forward. I, for one, can't remember this stuff to save my life without looking it up every time so the next section will provide a brief overview of the necessities.

Quick math background

As I mentioned above there are two solution methods that OpenFOAM employs to solve Equation 1. Here we will briefly review the math behind both options.

Option #1: The analytical solution

Since the ode in question is linear and first order it is possible to actually write down the analytical solution, something so rare in differential equations that it's worth getting a little excited about. The solution for the Equation 1 above is
$$ \phi (t) = \frac{\alpha}{\beta} + (\phi_0 - \frac{\alpha}{\beta}) \exp^{-\beta t} \tag{2} $$
This makes the assumption that \( \phi_0 = \phi(0) \).

Option #2: The implicit Euler method

Another possibility is to tackle this thing numerically and use, e.g., an implicit Euler method. Since it's non linear an iterative solution method is not necessary and this will go pretty quickly. The implicit numerical solution does not provide a function like the analytical solution but instead determines the value of \( \phi(t_1) \) when the value of \( \phi(t_0) \) is known and the time step \( dt = t_1 - t_0 \) is provided.
$$ \phi_1 = \frac{\phi_0 + \alpha dt}{1 + \beta dt} \tag{3} $$

This approach is a little less powerful in general than the analytical solution since it doesn't just give us an analytical function to use like Option #1. We will see though that in the context of the the overall simulation solution procedure we will only be advancing solution a time step \( dt \) anyway so that difference at least is not really a big deal.

OpenFOAM implementation

Okay, enough math. As I mentioned earlier both integration methods, analytical and Euler implicit are implemented in OpenFOAM. Both are represented as classes and they share a base class named IntegrationScheme. This base class includes a pure virtual function named integrate with the psuedo-code definition

integrateResult integrate(phi, dt, alphaBeta, beta)

The first argument, phi, is the initial value of phi at the start of the integration time, this is analogous to \( \phi_0 \) in the math section above. The time step over which the function integrates is dt.

The constants are a little confusing. The alphaBeta parameter in the code is the same thing as the \( \alpha \) constant in the math above, i.e. alphaBeta is the additive constant. The beta parameter is analogous to \( \beta \), i.e. it is the multiplicative coefficient for \( \phi \).

The return type, integrateResult is another helper class that in addition to some access functionality houses an array with two values. The first is the integrated result phi_1 that is analogous to \( \phi(t_1) = \phi_1 \). The second is the average value of \( \phi \) over the integration time \( dt \), phi_av.

Analytical class implementation

The analytical integration is implemented in the Analytical class. The integrate function is implemented, here simplified and in C++ like pseudo-code, as

// Analytical.C
integrationResult Analytical::integrate
const scalar phi,
const scalar dt,
const scalar alphaBeta,
const beta  
  integrationResult retvalue;
  scalar expTerm = exp(- beta*dt);

  scalar alpha = alphaBeta/beta;
  retValue.average() = alpha + (phi - alpha)*(1 - expterm)(beta*dt);
  retValue.value() = alpha + (phi - alpha)*expTerm;

  return retValue;

There are three things to see here:

  1. The solution is just as we saw in Equation 2. The trick is to not get bogged down in the alphaBeta/beta non-sense, it is the same thing as what we, on the math side, were calling \( \alpha/\beta \). That is to say that the variable alpha in the code is analogous to \( \alpha / \beta \) in our math equations above. With that clear you can see that the equations are identical.

  2. We also swap a \( t \), on the math side, for a \( dt \) on the code side. This is cleared up when you assume that for each of these integration periods the time resets. So if we were applying our math equation here we would always start from \( t=0 \), then the time interval \( dt = t - 0 = t \).

  3. The mean is calculated by integrating the analytical solution (Equation 2) with limits of \( 0 \) to \( dt \) and then dividing that integral by the time over which it was integrated, i.e. \( dt \). Some real Calculus 1 coming back to haunt me.

Euler class implementation:

Just like the analytical solution the implicit Euler method is implemented in its own class Euler. I'll reproduce the salient parts here again somewhat simplified and in a pseudo-code:

// Euler.C
integrationResult Analytical::integrate
const scalar phi,
const scalar dt,
const scalar alphaBeta,
const beta  
  integrationResult retValue;

  retValue.value() = (phi + alphaBeta*dt)/(1.0 + beta*dt);
  retValue.average() = 0.5 * (phi + retValue.value());

  return retValue;

Hopefully this is pretty straight forward:

  1. The value is calculated using an implicit Euler method just as in Equation 3.

  2. The average value here is just calculated as the arithmetic mean of the initial value, \( \phi_0 = \phi(t=t0) \), and the value at the end of the time step, \( \phi_1 = \phi(t=t_1) \), for \( t_1 - t_0 = 0 \). So \( \phi_{avg} = (\phi_1 - \phi_0)/2 \).


Hopefully this has helped you understand what these functions are doing and why, from a mathematical standpoint, they implement the equations that the do. I will say that the actual formulation of the equations in ReactingMultiphaseParcel is still a little confusing to me and is often further complicated by their attempts to linearize source terms. But at least with this knowledge you can hopefully back out what the terms in the equations are and if need be adjust them for your own purposes. Please comment if you have any questions or if you catch a mistake in my work.

I am also not honestly sure why both of these options exist. I used both in my work because the OpenFOAM default settings, in the coalChemistryFoam case, are to use one for the velocity and the other for the temperature (I can't remember which). It seems odd that you wouldn't just want to use the analytical solution all the time, but maybe the implicit Euler is just a little bit quicker to evaluate. If anyone knows the benefits of one of these options over the other please comment.