Assessment of Numerical Accuracy and Parallel Performance of OpenFOAM and its Reacting Flow Extension EBIdnsFoam

OpenFOAM is one of the most widely used open-source computational fluid dynamics tools and often employed for chemical engineering applications. However, there is no systematic assessment of OpenFOAM’s numerical accuracy and parallel performance for chemically reacting flows. For the first time, this work provides a direct comparison between OpenFOAM’s built-in flow solvers as well as its reacting flow extension EBIdnsFoam with four other, well established high-fidelity combustion codes. Quantification of OpenFOAM’s numerical accuracy is achieved with a benchmark suite that has recently been established by Abdelsamie et al. (Comput Fluids 223:104935, 2021. https://doi.org/10.1016/j.compfluid.2021.104935) for combustion codes. Fourth-order convergence can be achieved with OpenFOAM’s own cubic interpolation scheme and excellent agreement with other high-fidelity codes is presented for incompressible flows as well as more complex cases including heat conduction and molecular diffusion in multi-component mixtures. In terms of computational performance, the simulation of incompressible non-reacting flows with OpenFOAM is slower than the other codes, but similar performance is achieved for reacting flows with excellent parallel scalability. For the benchmark case of hydrogen flames interacting with a Taylor–Green vortex, differences between low-Mach and compressible solvers are identified which highlight the need for more investigations into reliable benchmarks for reacting flow solvers. The results from this work provide the first contribution of a fully implicit compressible combustion solver to the benchmark suite and are thus valuable to the combustion community. The OpenFOAM cases are publicly available and serve as guide for achieving the highest numerical accuracy as well as a basis for future developments.


Introduction
Numerical simulations have proven to be a helpful and reliable tool in science and engineering. One of the most widely used open-source computational fluid dynamics (CFD) libraries is OpenFOAM (Weller et al. (1998), The open surce CFD toolbox www. openf oam. com). However, there is no systematic assessment of OpenFOAM's numerical accuracy and computational performance in the context of chemically reacting flows.
Due to its use of unstructured computational meshes with arbitrarily shaped cells, OpenFOAM is most commonly used for Reynolds-averaged Navier-Stokes (RANS) and large eddy simulations (LES) involving complex geometries of real engineering applications, where the numerical accuracy is typically limited to second order at best. However, if well-conditioned or structured meshes are used, OpenFOAM offers built-in higher-thansecond-order discretization schemes. Because of this, OpenFOAM has been used as a tool for direct numerical simulations (DNS) across various fields of science (Komen et al. 2014;Habchi and Antar 2015;Addad et al. 2015;Lecrivain et al. 2016;Chu and Laurien 2016a, b;Tufano et al. 2016;Vo et al. 2016Vo et al. , 2017Wang et al. 2017Wang et al. , 2018Komen et al. 2017;Bricteux et al. 2017;Zhong et al. 2017;Zheng et al. 2018;Tufano et al. 2019;Zheng et al. 2019;Zirwes et al. 2020). In this context, the term quasi-DNS is sometimes used to highlight that OpenFOAM's spatial discretization schemes are essentially limited to fourth-order accurate interpolation schemes and temporal discretization is limited to second order accuracy. Most publications assessing OpenFOAM's accuracy focus on newly developed numerical tools based on the OpenFOAM library, e.g. Hassanaly et al. (2018); Komen et al. (2020); Gärtner et al. (2020), but rarely on its built-in capabilities (Lee (2017); Noriega et al. (2018)). It is therefore important to not only assess OpenFOAM's numerical accuracy, but also to provide a direct comparison with other, well established high-fidelity codes for chemically reacting flows.
Even though OpenFOAM is often used for the simulation of combustion phenomena, it has some important shortcomings. Firstly, OpenFOAM does not possess the ability to compute detailed molecular transport coefficients from kinetic gas theory for multi-component mixtures. The accurate modeling of diffusion coefficients and heat conductivity is however important to capture the correct internal structure of flames (Poinsot and Veynante 2001). Secondly, detailed simulations of three-dimensional flames are computationally very expensive and therefore demand an efficient numerical implementation. The reasons being the large range of time and length scales governing combustion processes (Peters 1999), as well as the the large number of chemical reactions and intermediate chemical species that have to be accounted for. Even the oxidation of simple fuels like methane is described by hundreds of chemical reactions (Smith et al. 1999) forming a numerically stiff system of equations (Law 2010). The tools offered by OpenFOAM for these types of problems lead to high simulation times . To overcome these shortcomings, the extension EBIdnsFoam has been developed for OpenFOAM, which implements detailed transport coefficients, provides several performance optimizations that can reduce total simulation times by up to 70 % for common combustion cases ) and achieves excellent parallel scalability for large-scale simulations.
Because there are no standardized numerical benchmarks for reacting flow solvers providing a full verification and validation chain, Abdelsamie et al. have recently proposed a benchmark suite (Abdelsamie et al. 2021) based on the well-known Taylor-Green vortex (TGV). This numerical test case has been used for the validation of numerical schemes of non-reacting flow solvers in the past decades (Orszag 1974;Brachet 1991;Shu et al. 1 3 2005), including a predecessor of OpenFOAM (Drikakis et al. 2007). The new combustion benchmark suite consists of four cases (here referred to as steps), from two-dimensional incompressible flow with an analytic solution to a hydrogen flame interacting with the three-dimensional Taylor-Green vortex. The benchmark is originally designed for low-Mach number codes and three high-fidelity codes were presented in the original publication. The benchmark suite has also been used to validate a hybrid Lattice-Boltzmann finite-difference solver (Hosseini et al. 2020).
The aim of this work is to systematically assess OpenFOAM's numerical accuracy and computational performance by providing extensive direct comparison with four well-established, high-fidelity solvers for reacting flows. The main questions addressed in this work are: • How do OpenFOAM's built-in higher-than-second order schemes perform compared to classical linear interpolation schemes? And how do they compare to dedicated highfidelity codes? • What is the parallel efficiency of OpenFOAM for both simple flow cases as well as reacting flows when using state-of-the-art chemistry acceleration techniques implemented into EBIdnsFoam? • What is the current state of code validation for combustion codes? Do different simulation codes yield the same results for flames developing within the TGV flow field?
The basis for this is the benchmark suite by Abdelsamie et al. the new results by Open-FOAM from this work provide an important contribution to the benchmark from the first compressible and fully implicit code. Throughout this work, three different flow types are considered: • Incompressible denotes the simulation of flows assuming constant density; • Low-Mach number treatment of flows refers to the use of a spatially homogeneous thermodynamic pressure computed from the domain-averaged temperature (Abdelsamie et al. (2021)), which is used in the computation of density, thermo-physical properties and chemical reaction rates; • Compressible treatment means that at each location of the computational domain, density and pressure are fully coupled so that the (thermodynamic) pressure field is inhomogeneous in space and allows for the generation of pressure waves.
Although all presented cases fall into the regime of low-Mach flows according to their Mach number, the numerical treatment of either low-Mach flows with spatially constant thermodynamic pressure or fully compressible flows has an impact on the reacting flow cases, later discussed in Sect. 4.4. The paper is structured as follows: Sect. 2 details the numerical methods available in OpenFOAM relevant for this work and explains the implementation details of the EBIdns-Foam extension. A brief description of the high-fidelity reacting flow codes used for comparison is given as well. Simple validation cases are presented in Sect. 3 that assess the convergence order of OpenFOAM and validate the velocity-pressure coupling with a pulsatile pipe flow setup. Two simulations of canonical combustion setups, i.e. zero-dimensional (0D) auto-ignition and one-dimensional (1D) freely propagating flame, are compared to the open-source library Cantera (Goodwin et al. 2021). The four steps of the new benchmark suite are conducted in Sect. 4. The first step is an incompressible simulation of the two-dimensional Taylor-Green vortex with analytic solution and the second step is 1 3 an incompressible simulation of the three-dimensional TGV where a pseudo-spectral DNS reference solution is available. The third step superimposes an inhomogeneous temperature field onto the TGV velocity field as well as a spatially varying mixture composition of the gas, thus simulating the mixing of an inhomogeneous gas mixture including diffusion of mass and energy. The fourth step simulates a flame interacting with the TGV and thus additionally considers chemical reactions. Section 5 summarizes the findings from this work.

Numerical Codes
The simulations presented throughout this work have been conducted with OpenFOAM. Additionally, results published in Abdelsamie et al. (2021) from three high-fidelity codes, DINO, Nek5000 and YALES2 are presented together with the OpenFOAM results. Lastly, new results from the compressible high-order code KARFS are analyzed in Sect. 4.4. These codes are briefly introduced in this section.

OpenFOAM
OpenFOAM provides a large number of different solvers for many types of applications and employs the finite volume method (FVM) on collocated meshes to discretize the governing equations. It can use unstructured, non-conforming meshes with arbitrary polyhedral cell shapes. The solution of the governing equations is treated implicitly, which makes OpenFOAM the only code with implicit time discretization discussed in this work. Two second-order time discretization schemes are available, the Crank-Nicolson method and the second-order backward (BDF2) scheme.
Because OpenFOAM applies the finite volume method, the spatial discretization consists of an interpolation from cell center values to the cell faces and a subsequent integration of the flux over the cell faces. For the cell face integration, a single integration point is used, which in general is second order accurate. For the interpolation to the cell faces, OpenFOAM offers a number of interpolation methods including the first-order upwind scheme and the second-order linear interpolation scheme. A fourth-order scheme is available as well, which is called "cubic" in OpenFOAM. It uses a cubic interpolation by taking into account the cell centroid values of the quantity as well as its explicitly computed gradients n f normal to the cell faces, according to The variable f is the value of interpolated to the cell face f, and C and N are the values in the cell centers of the current (C) and neighbor (N) cell. The parameter = is the normalized distance between the cell face and centers. The polynomial of degree three in Eq. (1) is rearranged so that the first part can be evaluated implicitly from the standard central difference scheme and the second part constitutes an explicit correction. Because of the explicit contribution, the scheme tends to be less stable and is only viable on well-conditioned meshes and with low Courant-Friedrichs-Lewy (CFL) numbers. (1) Another higher-than-second-order discretization method is given by the recently developed third-party weighted essentially non-oscillating (WENO) scheme library for Open-FOAM (Gärtner et al. 2020) based on the work of Martin and Shevchuk (2018). It implements the WENO concept, which introduces a polynomial representation of the cell values. This library is used for step 1 and step 2 of the benchmark suite for the convective term of the momentum balance equation with a polynomial of degree four (WENO4).
The time step of the simulation is adjusted according to the Courant-Friedrichs-Lewy (CFL) and Fourier (Fo) number, where the CFL number is defined as with u being is the local convective velocity, Δt the time step width and Δx the length scale of the computational cell. The Fourier number is evaluated from where a is the thermal diffusivity. It should be noted that even for the compressible solvers, the CFL number is the convective CFL number and thus evaluated with the local flow velocity, not the speed of sound. For additional discussion on this, see Sect. 4.4.
All OpenFOAM simulations presented in this work have been performed with the versions v1712 and v2006. More specifically, all simulations with EBIdnsFOAM are performed with OpenFOAM v1712 as well as the incompressible simulations in step 1. The simulation applying the WENO scheme in step 1 as well as all simulations in step 2 are run with OpenFOAM v2006. The simulations from step 1 and 2 have been performed additionally with the latest versions v2106 and 9, but no significant differences in terms of the numerical accuracy have been observed.

pimpleFoam
For the simulation of incompressible flows in the following sections, OpenFOAM's standard solver pimpleFoam is used. It solves the Navier-Stokes equations assuming constant density and applies the PIMPLE algorithm, which combines the PISO (Issa (1986)) and SIMPLE (Caretto et al. 1973) algorithms, for pressure-velocity coupling. An extension of pimpleFoam to compressible flows is given by OpenFOAM's rhoPimpleFoam solver, which solves an additional energy equation and computes the density from an equation of state.

EBIdnsFoam
As described in the introduction, OpenFOAM's built-in reacting flow solvers do not support the computation of detailed diffusion coefficients. Because of this, the extension EBIdnsFoam (Solver for detailed and efficient simulation of reacting flows https:// vbt. ebi. kit. edu/ 594. php, Zhang et al. 2015Zhang et al. , 2016 has been developed at the Engler-Bunte-Institute. It couples OpenFOAM's finite volume capabilities for the solution of the fully compressible Navier-Stokes equations and the balance equations for the species masses and mixture energy together with Cantera's routines for computing detailed thermo-physical and transport properties. Cantera (Goodwin et al. 2021) is an open-source thermo-chemical library commonly used in the combustion community. Figure 1 illustrates this coupling. In this way, EBIdnsFoam behaves numerically in the same manner as the standard react-ingFoam solver but uses Cantera to provide the transport coefficients. The governing equations implemented into EBIdnsFoam are (Poinsot and Veynante (2001)): • conservation of total mass • conservation of momentum • balance of species masses • balance of energy Here, is the gas mixture density, t is time, ⃗ u is the bulk fluid velocity, p is the pressure, Y k is the mass fraction of species k, I is the identity tensor, is the dynamic viscosity of the gas mixture, ̇k is the reaction rate of species k and N the total number of species. The correction velocity ⃗ u c assures that the sum of all diffusive fluxes ⃗ j k is zero. In the balance equation for energy, c p is the isobaric heat capacity of the gas mixture and T the temperature. The enthalpy of formation h 0 k of species k and the sensible enthalpy of the mixture h s for ideal gases are computed from NASA polynomials. The diffusive heat flux is computed from where is the heat conductivity of the gas mixture and h s,k is the sensible enthalpy of the k-th species. EBIdnsFoam has been used in previous studies, including the simulation of partially premixed flames Hansinger et al. 2020;Soysal et al. 2019), Bunsen and slot burner flames Chen et al. 2021), flamewall interaction (Steinhausen et al. 2020;Zhang et al. 2021;Steinhausen et al. 2022;Zhang et al. 2022;Driss et al. 2022;Zhu et al. 2023), turbulent flames with instrinsic instabilities , combustion noise , spherically expanding flames Eckart et al. 2022), as well as ignition phenomena Häber et al. 2017;Wang et al. 2021Wang et al. , 2022, soot formation (Wang et al. 2023) and combustion under rocket engine conditions (Martinez-Sanchis et al. 2022;Martínez-Sanchis et al. 2022). For more information about the implemented governing equations and code details, we refer the reader to Zirwes et al. (2020); Zhang et al. (2015); .
Another focus in the development of EBIdnsFoam lies on computational performance, especially in the context of high-performance computing. One key optimization is the use of automatically generated routines containing highly optimized source code for the computation of chemical reaction rates ), which is able to reduce total simulation times by up to 70 % without affecting the accuracy of the simulation. Due to the high numerical stiffness generally introduced by the chemical reaction rates, an operator splitting approach is used in EBIdnsFoam. The open-source integrator CVODE by Sundials (Hindmarsh et al. (2005)) is used to integrate the chemical reaction rates over the simulation time step and provides time-step averaged source terms for the balance equations.
In the next sections, results of OpenFOAM's standard solvers and EBIdnsFoam are compared with other high-fidelity codes. These are described briefly in the following subsections. For a more detailed description, see the original benchmark suite publication (Abdelsamie et al. 2021) and the references given below. A brief overview of the discussed codes and their most important features is given in Table 1. DINO (Abdelsamie et al. 2016;Chi et al. 2017) is a three-dimensional, finite differences (FD) based computational fluid dynamics code that supports incompressible and low-Mach number flows. In the context of the Taylor-Green vortex benchmark suite, spatial derivatives are evaluated from a sixth-order central difference scheme and time integration is performed using an explicit fourth-order Runge-Kutta (RK) method. Thermophysical and transport properties are computed with Cantera 2.4.0.

Nek5000
The solver based on Nek5000 (2017), which is used in the Taylor-Green benchmark suite, has been developed at ETH Zürich and applies the spectral element method (SEM) to solve incompressible flows and low-Mach number chemically reacting flows. Its spatial discretization order is 7th-15th order, with a semi-implicit BDF3 scheme for temporal discretization.

YALES2
YALES2 (Moureau et al. 2011) simulates low-Mach number flows with the finite volume method, supporting unstructured meshes. It applies a fourth-order accurate node-based centered scheme for spatial discretization and uses a fourth-order explicit Runge-Kutta method for time integration.

KARFS
The KAUST Adaptive Reacting Flow Solver (KARFS) (Hernández Pérez et al. 2018) simulates fully compressible reacting flows. It uses the finite difference method with an 8thorder spatial discretization and explicit RK4 method for time integration. Its programming model allows work to be offloaded to GPUs and other hardware accelerators.

Validation Cases for pimpleFoam and EBIdnsFoam
Before the benchmark suite by Abdelsamie et al. (2021) is discussed, a few additional validation cases for both OpenFOAM's standard solvers and EBIdnsFoam are presented.

Convection-Diffusion Equation
A simple one-dimensional steady-state convection-diffusion equation, is used to quantify OpenFOAM's convergence order. Pe is the Péclet number and is an arbitrary scalar quantity. The boundary conditions are (x = 0) = 0 and (x = L) = 1 , and the analytic solution is then given by exact (x) in Eq. (9). The numerical mesh consists of equidistantly spaced cells with size Δx . The case has been run with OpenFOAM's scalar-TransportFoam solver. Figure 2 shows the L 1 error for different mesh resolutions, which is defined as where N is the number of grid points. Running the simulation with the first-order upwind scheme yields the expected convergence order of one. By using the linear and the cubic scheme (see Sect. 2.1), a second-order convergence is achieved. Even though the cubic scheme employs a fourth-order accurate cell face interpolation, the overall scheme is limited by the second-order face flux integration. However, the cubic scheme is less dissipative than the linear interpolation scheme, resulting in a lower L 1 error.

Pulsatile Pipe Flow
In the next setup, the pressure-velocity coupling from OpenFOAM's standard pimple-Foam and rhoPimpleFoam solvers are validated. The numerical setup is built to resemble an experimental setup for pulsatile laminar pipe flows Ünsal 2008). Air enters a pipe with a length of L∕D = 714 . The velocity profile at the inlet is set to a block profile. The mean mass flow rate is set according to a Reynolds number of Rēu = 2947 . The mass flow rate at the inlet oscillates harmonically with a prescribed amplitude m ⋆ A and frequency f, as shown in Fig. 3. The pipe is considered as a two-dimensional axisymmetric case, with a radial resolution of D∕Δr = 100 and an axial resolution of L∕Δx = 1000 . The numerical settings are summarized in Table 2.
As the air enters the domain and flows along the pipe, the velocity profile changes from the block profile given at the inlet to the parabolic profile for fully developed pipe flow. To validate the pressure-velocity coupling, the amplitude ratio of the mass flow rate ṁ and pressure gradient P = p∕ x , as well as the phase shift between mass flow rate and pressure gradient Δ are evaluated. Figure 4 shows an example of the time signal of the mass flow rate ṁ normalized by the mean mass flow rate ̄ṁ and the time signal of the pressure gradient P normalized by the mean pressure gradient P M . The has to increase in order to support the same normalized mass flow rate amplitude ṁ * A . At the same time, the phase shift between the two signals shifts from 0 • to 90 • with increasing frequency.
An analytic solution for fully developed, pulsatile pipe flows is available as function of the dimensionless frequency F, where R is the pipe diameter and the viscosity. The analytic solution reads (Lambossy (1952), Grace (1928), Ünsal (2008), Haddad et al. (2010), Ünsal et al. (2005))  where i is the imaginary unit and J n is the n-th Bessel function of the first kind. The amplitude ratio can then be obtained from and the phase shift from where ℜ(Ψ) is the real part and ℑ(Ψ) the imaginary part of Ψ.
The pulsatile pipe flow case has been run with OpenFOAM's incompressible flow solver pimpleFoam and its compressible equivalent rhoPimpleFoam for different F. The evaluated amplitude ratio of the mass flow rate and pressure gradient as well as the phase shift between mass flow rate signal and pressure gradient signal show excellent agreement with the analytic solution, as shown in Fig. 5. In this way, OpenFOAM's pressure-velocity coupling implementation is validated.

Canonical Combustion Cases
In order to validate the implementation of chemical reaction rates and detailed molecular diffusive fluxes in the EBIdnsFoam solver, two canonical cases for combustion solvers are presented.

0D Auto-Ignition
The first canonical case for reacting flows is the zero-dimensional (0D) auto-ignition of methane in air. A homogeneous, constant-volume batch reactor is filled with the unburnt gas mixture at a temperature above the auto-ignition temperature. In this case, the initial temperature is set to 1600 K and the initial pressure to 1 bar. Hydrogen is added to the gas mixture to speed up the ignition process. The initial conditions are shown in the table in Fig. 6. The reaction mechanism is GRI 3.0 (Smith et al. 1999). This case evaluates the accuracy of the temporal numerical integration of the chemical source terms. As the simulation starts, the mixture begins to ignite and transitions to the burnt state. The temporal evolution of species mass fractions obtained from EBIdnsFoam is depicted in Fig. 7 together with reference results from Cantera. The maximum numerical difference between the EBIdnsFoam and Cantera solutions are below 0.1 %.

1D Freely Propagating Flame
The second case is a one-dimensional freely propagating flame.
A stoichiometric mixture of the unburnt methane and air is introduced at the inflow plane on the left of the computational domain. The inlet velocity is set to the laminar flame speed s L,0 to obtain a steady-state solution. The left half of the domain is initially filled with the unburnt gas mixture and the right half is filled with the combustion products at chemical equilibrium. The reaction mechanism is again GRI 3.0, and the diffusion model is the mixture-averaged model (Curtiss-Hirschfelder approximation) (Kee et al. 2005). The computational setup and boundary conditions are summarized in Fig. 8. Figure 9 shows the spatial profiles of species mass fractions at the steady-state solution. Again, the agreement with Cantera is very good, with largest differences below 0.5 %.

Taylor-Green Vortex Benchmark Suite
In this section, the four steps of the new benchmark suite for reacting flow solvers by Abdelsamie et al. (2021) are conducted with OpenFOAM's standard solver pimpleFoam (steps 1 and 2) and the custom solver EBIdnsFoam (steps 3 and 4). For steps 1-3, results are directly compared to the low-Mach codes DINO, Nek5000 and YALES2 to give a perspective on the accuracy of OpenFOAM. For step 4, the compressible code KARFS is additionally used to highlight the difference between the low-Mach and compressible

Step 1: 2D Incompressible Single-Component TGV
The first case is an incompressible, two-dimensional Taylor-Green vortex flow. The domain is a square with side lengths L and periodic boundary conditions. The velocity field is given by the analytic solution  Figure 11 shows the x-component of velocity along the centerline at y∕L = 0.5 (left) and the y-component of velocity along the centerline at x∕L = 0.5 (right) at t = 10 ref . The black line in the background denotes the analytic solution. The results from DINO and YALES2 are obtained using a mesh resolution of 64 2 cells or grid points, and the results of Nek5000 are reported for 8 2 spectral elements which is equivalent to 65 2 grid points. For the OpenFOAM (OF) solutions, a mesh with 65 2 cells was used so that the domain centerlines cross through the cell centers of the mesh. Two OpenFOAM simulations are conducted: one where all spatial derivatives are discretized with the cubic scheme (OF cubic), and one where the WENO4 scheme is used for the convective term of the momentum equation (OF WENO4). Visually, all simulation codes are in perfect agreement. Even the zoomed view on the left of Fig. 12 shows no visual differences. Because of this, Table 3 compares the maximum velocity at t = 10 ref from all codes with the analytic solution. Using the same mesh and the same time step, the peak velocity values of OpenFOAM lie within 1.8 × 10 −2 % and 1.1 × 10 −3 % of the analytic solution. In this way, the achieved accuracy is of the same order of magnitude as the one from DINO and YALES2. Appendix 29 provides additional velocity profiles at coarser mesh resolutions.
Before moving to the next step of the benchmark suite, the influence of various discretization schemes from OpenFOAM as well as the achieved convergence order are discussed. Figure 12 shows the same velocity profile as in Fig. 11 on the left, but all results shown in the Figure are now obtained from OpenFOAM using different spatial discretization schemes. It is  clearly visible that the first-order upwind scheme is too dissipative. By zooming into the peak velocity region, the solution obtained with the linear scheme also shows visible deviation from the analytic solution, whereas the cubic and the WENO4 schemes show excellent agreement with the analytic solution.
Due to the existence of an analytic solution, the spatial convergence order of OpenFOAM can be determined. Figure 13 shows the L 1 (see Eq. (10)) and L 2 error for the u x velocity profile at t = 10 ref for different mesh resolutions obtained with the cubic discretization scheme. The L 2 error is evaluated from For coarse meshes, the convergence order is approximately 2 (see also Sect. 3.1). However, at finer meshes, the convergence order approaches 4th order. This explains the overall good agreement with the high-fidelity codes shown above.

Step 2: 3D Incompressible Single-Component TGV
Step 2 of the benchmark suite is the incompressible simulation of the three-dimensional Taylor-Green vortex. The computational domain is a cube with side length L = 2 L 0 and L 0 = 1 m. All faces of the cube are defined as periodic boundary conditions, and the initial velocity field is given by with u 0 = 1 m/s. The viscosity is again set to = 1 1600 m 2 /s, leading to a Reynolds number of Re = u 0 L 0 ∕ = 1600 . The characteristic time scale is ref = 1 s. All OpenFOAM simulations performed for the 3D TGV cases are run with outer PIMPLE iterations until initial residuals of all quantities are below 5 × 10 −10 or up to a maximum of 10 outer iterations. For this case, there is no analytic solution, but a reference pseudo-spectral DNS on a mesh with 512 3 grid points is available, which is denoted here as RLPK (Van Rees et al. 2011). This and all subsequent 3D TGV cases are run without sub-grid turbulence models. While the resolution of 512 3 grid points can be seen as DNS resolution (Van Rees et al. 2011), simulations at lower resolution would correspond to an underresolved DNS or LES without explicit turbulence model. However, since the decaying flow field of the TGV stays symmetrical throughout the whole simulation time considered here (pseudo-turbulence), the terms underresolved LES or implicit sub-grid model can be misleading and are omitted from the description of all cases in this paper. Figure 14 shows iso-surfaces of the Q-criterion at the beginning of the simulation (left). After about 10 s, the flow decays into pseudo-turbulence, as shown on the right of Fig. 14.
The time step of the simulations is dynamically adjusted to ensure CFL ≤ 0.3 . The computational mesh for the OpenFOAM simulations consists of 513 3 equidistant cells to ensure that the centerlines pass through the cell centroids, while DINO and YALES2 use 512 3 cells or grid points. Figure 15 shows the velocity profiles of the x (left) and y (right) components of velocity along the centerlines at t = 12.11 ref from the different codes. The reference solution from RLPK is not available for this specific data set. The simulation with OpenFOAM is again conducted with purely cubic discretization schemes (OF cubic) and a separate simulation is performed using the WENO4 scheme for the convective term of the momentum equation. Visually, the velocity profiles from OpenFOAM agree very well with the three high-fidelity codes. The zoomed region in Fig. 15 on the left shows very good agreement between DINO and Nek5000, while YALES2 and the OpenFOAM results show small deviations. The zoomed region to the peak velocity on the right of Fig. 15, shows a good agreement between all codes, with small deviations by DINO. Again, the results by OpenFOAM compare very well to the high-fidelity codes.
Another important validation is the temporal evolution of the volume-averaged kinetic energy and dissipation rate within the computational domain, defined as These quantities are shown in Fig. 16 together with the reference solution from the pseudospectral DNS (RLPK). Both the temporal evolution of the kinetic energy and of the dissipation rate show very good agreement between all codes (see also appendix B for additional data). The only deviation is observed by YALES2, as shown in the zoomed region on the right of Fig. 16. The influence of the choice of discretization scheme in OpenFOAM is shown in Fig. 17. All simulations are run with 513 3 cells, but with different discretization schemes. When looking at the volume-averaged dissipation rate, it is immediately clear that the first-order upwind scheme is not sufficient to describe the dynamics of the TGV. The second-order central difference scheme (CD) is less dissipative, but also shows significant deviations from the reference solution (zoomed region). With the cubic and WENO4 schemes, the reference solution is recovered. Since the dissipation rate in Fig. 17 is computed from the stress tensor, which in turn is computed from the gradient of the velocity field (Eq. (22)), Fig. 16 Normalized volume-averaged kinetic energy (left) and dissipation rate (right) over time for the 3D incompressible TGV from step 2

Fig. 17
Volume-averaged dissipation rate over time on a mesh with 513 3 grid cells and different OpenFOAM discretization schemes for the 3D incompressible TGV from step 2 simulations run with schemes having larger numerical dissipation lead to smoother velocity gradients and therefore to lower predicted physical dissipation rates.
Finally, the influence of the mesh resolution is discussed. Figure 18 shows the volumeaveraged dissipation rate obtained by OpenFOAM with its cubic discretization scheme as well as DINO with its 6th-order spatial discretization for 64 3 , 128 3 and 256 3 cells or grid points. For 256 3 points, DINO's solution agrees very well with the reference solution and is less dissipative compared to that from OpenFOAM. In order to test the solution with less grid points, an 8th-order explicit filter was activated in DINO, which is not needed in OpenFOAM. Looking at the results in Fig. 18, it can be observed that OpenFOAM with 64 3 , and 128 3 cells shows less dissipative solutions compared to the ones from DINO with filter. In this way, the comparison with the high-order code DINO again shows the good numerical performance of OpenFOAM's built-in cubic scheme.

Step 3: 3D Multi-component Non-reacting TGV
In steps 1 and 2, the incompressible flow solver pimpleFoam from OpenFOAM has been used. It achieved excellent agreement with both the analytic solution in step 1 as well as with the three high-fidelity codes in steps 1 and 2. In step 3, the EBIdnsFoam solver is used, as the three-dimensional TGV velocity field is superimposed by a temperature field, so that an incompressible or constant-density approach is no longer valid. The concentration profiles for N 2 , O 2 and H 2 are set according to where Y k are the mass fractions of species k. The shape function Ψ is defined as The characteristic length scale is set to L = 2 mm and the velocity to u 0 = 4 m/s. The temperature at each position is set to the temperature of the chemical equilibrium according to the local species mixture at 1 atm and 300 K. The thermo-physical properties of the chemical species are given by the thermodynamic quantities associated with the reaction mechanism by Boivin et al. (2011). Figure 19 shows the initial profiles of H 2 , O 2 and T along the centerline. The simulations are performed on a mesh with 256 3 cells. The time step is dynamically adjusted to ensure CFL ≤ 0.25 and Fo ≤ 0.15 . Chemical reactions are turned off in this step and the mass diffusion of the species is treated using constant but non-equal Lewis numbers. For more information, see Abdelsamie et al. (2021).
Figure 20 (left) shows the temporal evolution of the maximum temperature inside the simulation domain. Since chemical reactions are disabled for this step, the mixing leads to 1 3 a decrease of the peak temperatures. The OpenFOAM-based EBIdnsFoam solver as well as the three high-fidelity codes display again a good agreement. On the right of Fig. 20, the spatial profile of the temperature along the y centerline at t = 0.5 ms is depicted. By zooming into the central region, minor differences can be observed between the codes, where the solutions from YALES2 and EBIdnsFoam overlap but are slightly lower than the results from the high-order codes DINO and Nek5000. The peak temperatures on the other hand agree within 1 % between all codes, as shown in Table 4. Figure 21 shows the u x and u y velocity profiles along centerlines at t = 0.5 ms.
The u x velocity profile shows again remarkable agreement while smaller deviations by about 3 % of the u y profile of EBIdnsFoam are visible near the peak values.  A comparison between the spatial profiles of species mass fractions of O 2 and H 2 is given in Fig. 22. While the oxygen mass fraction profile yields close agreement with the three high-fidelity codes, there are considerable deviations between EBIdnsFoam and the other three codes for the profile of the hydrogen mass fraction. One major difference between EBIdnsFoam and the three other codes considered so far is that EBIdnsFoam is a compressible code and the other three codes solve for low-Mach number flows. Because of this, the simulation of the non-reactive mixing case has been repeated with the mixtureaveraged diffusion model instead of the constant Lewis number approach in order to compare with results from the high-order KARFS code, which is also fully compressible and employs the mixture-averaged diffusion approach. Figure 23 shows again the mass fraction profile of hydrogen with the same values from DINO as before for reference, but this time with results from EBIdnsFoam using the mixture-averaged diffusion model (mix-avg) and the KARFS code using the mixture-averaged model as well. This time, the agreement for the hydrogen profile between KARFS and EBIdnsFoam is very good. This is further discussed in the next subsection.

Step 4: 3D Multi-component Reacting TGV
The most challenging numerical setup of the benchmark suite is step 4. Compared to step 3, chemical reactions are enabled and 9 species as well as 12 chemical reactions from the Boivin reaction mechanism for hydrogen combustion are considered. The initial conditions for the species are changed, as now not only the temperature is set according to the chemical equilibrium at each position, but also the species concentrations (see Abdelsamie et al. (2021) for reference). The reacting TGV case corresponds to two, initially planar, diffusion flame fronts. Due to the TGV flow field, the flame fronts are distorted. Initially, the local mixture composition is set to the chemical equilibrium at each grid point. Because of this, heat release rates, even within the hot temperature region, are initially zero. As the fuel and oxygen species start to diffuse into the high temperature region and the profiles are convected by the TGV flow field, an ignition event occurs and a net heat release rate establishes that forms the diffusion flame fronts. Again, the mesh consists of 256 3 cells with CFL ≤ 0.25 and Fo ≤ 0.15 . Because a flame propagates during the simulation, the gas inside the simulation domain expands, which leads to an increase of the pressure, as the periodic boundary conditions act as a closed domain for the global pressure. As a consequence, the low-Mach number codes DINO, YALES2 and Nek5000, which treat the thermodynamic pressure as homogeneous throughout the domain (but variable in time), show a different behavior as explored below. Figure 24 shows the temporal evolution of the maximum temperature inside the simulation domain. After about 0.1 ms, there are significant discrepancies depending on the employed simulation model: the blue dash-dotted line is the result from DINO, which represents the low-Mach number solvers using the constant Lewis number diffusion approach (const Le); the green dashed line is the compressible EBIdnsFoam solver, using constant Lewis number diffusion as well. However, these two results are not comparable due to the low-Mach vs. compressible approach. On the other hand, the results from KARFS and EBIdnsFoam are almost identical, as in this case both codes use the mixture-averaged diffusion model and utilize a compressible flow formulation. The spatial temperature profile along the y centerline at t = 0.5 ms is less sensitive to the choice of compressible or low-Mach solver and also to the choice of diffusion model, as shown in Fig. 25. Note that the maximum temperature within the domain, that is shown in Fig. 24 is not located on the line from Fig. 25 but at a different location inside the simulation domain.
Although Fig. 24 suggests that the two compressible codes EBIdnsFoam and KARFS should compare well when considering other variables for the reacting Taylor-Green vortex case, an evaluation of the velocity profile reveals noticeable differences. While the profile of u y along the centerline at t = 0.5 ms shows a good agreement in all cases on the right of Fig. 26 and the two EBIdnsFoam cases agree very well with DINO for u x (Fig. 26 on the left), the velocity profile of KARFS shows a considerable deviation. One difference is the treatment of the compressibility or velocity-density-pressure coupling. Since KARFS is an explicit, fully compressible code, its time step is based on the acoustic CFL number and thus resolves the evolution of pressure waves in the simulation. On the other hand, due to OpenFOAM's implicit treatment of the governing equations and the PIMPLE method, OpenFOAM's stability criterion depends on the convective CFL number, so that time steps are usually much larger than what would be needed to resolve pressure waves. The difference in time step size also affects the operator splitting for the coupling of the flow and chemical source terms. However, an additional simulation with OpenFOAM applying the acoustic CFL number as time step limiter yields the same results as with the convective CFL number. The difference in the velocity profile is also visible when comparing other quantities, like the mass fraction of OH and the heat release rate (HRR), as shown in Fig. 27. Although the results between EBIdnsFoam and KARFS are similar, a significant deviation remains. Since the velocity profiles of EBIdnsFoam and DINO match remarkably well, but the maximum temperature evolution between KARFS and EBIdnsFoam also matches remarkably well, it is not immediately clear which result, if any, is more accurate.
While the uncertainty reported in Abdelsamie et al. (2021) for the reacting TGV case between the three low-Mach number codes is about 1 %, it is much higher for compressible reacting flow solvers. A similar discussion has arisen during the preliminary presentation of benchmark results at the 17th International Conference on Numerical Combustion (May Fig. 26 Velocity profiles along the centerlines at t = 0.5 ms for step 4 1 3 6-8, 2019) in Aachen, Germany, without any final conclusion. It is therefore important, to include more compressible reacting flow solver results for this benchmark suite to aid the improvement of existing reacting flow solvers and to provide a reliable numerical benchmark for future codes. All results and setups from the OpenFOAM/EBIdnsFoam simulations are publicly available (OpenFOAM Data for the TGV benchmark suite for reacting flows https:// github. com/ g3bk47/ TGV-OpenF OAM-bench mark).

Parallel Performance
The computational and parallel performance of the different codes can be evaluated by tracking the time required for performing one simulation time step for the different steps of the benchmark suite in dependence on the number of parallel processes. These performance measurements are presented in Fig. 28 for the simulations from step 2-4. The numbers near the data points inside the figures indicate the parallel scaling efficiency E n where n 0 is the number of parallel processes or MPI ranks of the first data point, and t the wall-clock time per time step. The measurements have been performed on the following supercomputers:  Fig. 28), the time per iteration is almost an order of magnitude higher than the other codes included in this comparison. The reason for this is that OpenFOAM is a general purpose code with focus on unstructured meshes with arbitrary cell shapes and thus cannot be as efficient as e.g. structured mesh solvers. In the depicted range, the parallel scaling however is very good, being in a super-linear regime.
As briefly described in Sect. 2.1, EBIdnsFoam has been designed with a focus on computational performance, specifically targeting the computation of thermo-physical properties and chemical reaction rates. Therefore, the time required per time step for steps 3 and 4 lies well within the range of the tested high-fidelity codes. The second to last data point for OpenFOAM on the bottom figure in Fig. 28 for step 4 is the performance measurement with 9728 MPI ranks. With that many parallel processes, OpenFOAM is the only code that retains near-perfect parallel scaling (E 9728 = 0.93 ). This confirms that the EBIdnsFoam extension is well suited for large-scale parallel simulations.

Summary and Conclusions
In this work we provide a thorough comparison and validation of one of the most used open-source CFD codes OpenFOAM. Although OpenFOAM is often used in the context of complex geometries that do not allow for high-order schemes, it supports fourth-order interpolation schemes. Specifically in the context of reacting flows, there are not many publications assessing the numerical accuracy of OpenFOAM. It is therefore important to provide well-designed benchmark cases and to perform direct comparisons with other well-established codes. Because of this, an extensive validation of OpenFOAM's built-in solvers and numerical schemes, as well as of external developments (a high-order WENO discretization library and the EBIdnsFoam extension for reacting flows) has been conducted in this work. For the first time, a full benchmark suite recently developed for reacting flow solvers has been performed with OpenFOAM in conjunction with direct comparison including four high-fidelity codes to systematically assess the numerical accuracy of OpenFOAM.
The main findings are summarized as follows: • If the simulation setup allows for hexahedral cells and low CFL numbers, it is strongly recommended to use OpenFOAM's cubic discretization scheme due to its higher accuracy, instead of defaulting to its "linear" scheme, which is commonly done in the Open-FOAM community. • The pressure-velocity coupling approach with OpenFOAM's PIMPLE algorithm has been validated by means of pulsatile pipe flows. Likewise, the implementation of diffusive fluxes and chemical reaction rates in EBIdnsFoam has been validated by comparison of canonical combustion cases with the open-source thermo-chemical library Cantera. • The numerical accuracy of OpenFOAM in the incompressible flow setups of steps 1 and 2 of the Taylor-Green vortex benchmark is excellent. Depending on the mesh reso-lution, a convergence order of up to four can be achieved. For the 2D setup, the deviation from the analytic solution is about 0.001 % and thus of the same order of magnitude as DINO and YALES2. Similar results are found for the 3D TGV case. • For the non-reactive mixing case in step 3, agreement is again very good. However, considerable deviations exist in the spatial profiles of the hydrogen mass fraction. This difference vanishes, when the compressible high-order code KARFS is used for comparison together with the mixture-averaged diffusion model. Therefore, these deviations might be attributable to the low-Mach vs. compressible approach. • The reactive TGV case in step 4 is sensitive to the flow treatment, i.e. low-Mach number or compressible flow. This, however, leads to additional uncertainties. While the temporal evolution of the maximum temperature agrees best between the compressible KARFS and EBIdnsFoam solvers, the spatial velocity profiles contained in the benchmark match best between the low-Mach DINO code and EBIdnsFoam. In general, the reacting 3D TGV case reveals different results between DINO, KARFS, OpenFOAM and YALES2, depending on which quantity is compared, i.e. velocity fields or mass fraction profiles. This highlights the need for more reliable benchmark cases, especially for compressible reacting flow solvers in the future. • The results presented in this paper constitute a valuable contribution to the benchmark repository, not only from one of the most widely-used open-source CFD codes, but also from the first fully implicit and compressible reacting flow code. All results and numerical setups for OpenFOAM used in this work are publicly available (OpenFOAM Data for the TGV benchmark suite for reacting flows 2023).