Development of Openfoam Solvers for Incompressible Navier–stokes Equations Based on High-Order Runge–kutta Schemes

Direct Numerical Simulations (DNS) and Large-Eddy Simulations (LES) not only require massively parallel computing facilities due to their high computational cost, but they also need accurate numerical methods. Particularly high-resolution space-time discretization schemes would be used to ensure a minimal numerical diffusion. An accurate numerics is typically adopted only in academic codes with a very limited dissemination to general public. Nevertheless open-source CFD codes provide suitable environments for implementation and testing low-dissipative algorithms needed for turbulence simulation. Moreover these codes produce a reliable tool to test high-fidelity numerics on unstructured grids, which are particularly appealing for industrial applications. Therefore in this work we have developed several solvers for incompressible Navier-Stokes equations (INS) based on low-dissipative schemes for time-integration.

In this work the finite volume code OpenFOAM (Open-source Field Operation And Manipulation), was selected as development environment. Very few research papers dealing with Runge-Kutta methods within OpenFOAM code have appeared in literature, [1,2]. With this research work we want to contribute in the OpenFOAM-RK methods development, with emphasis to incompressible flows, which are early investigated only by Vuorinen et al. [1].

It is very interesting to remark as the application of RK methods to the INS equations is not a trivial task because of the DAE nature of the discretized equations. Many CFD practitioners explicitly advance the velocity at each RK stage using standard techniques for ODE systems, thus the Poisson equation for the pressure is solved to guarantee a divergence-free velocity field. Anyhow, the impact of this projection-type approach on the temporal order of accuracy is not clear [3]. In many papers the accuracy of the velocity is often darkly assumed to be unaffected by the DAE nature of INS. Besides the temporal accuracy of the pressure is often not reported.

For the previous reasons we have considered explicit RK-projected-type schemes specifically developed by Sanderse and Koren [3], for index 2 DAE systems, i.e. Navier-Stokes equations. It is worth noting that the projected nature of the previous schemes allows to avoid the use of Rhie-Chow correction, thus the previous schemes result efficient and easy to implement. Furthermore L-stable Singly Diagonally Implicit Runge-Kutta (SDIRK) schemes were also adopted. In this specific case, as will be described in the next sections, an iterated PISO-like procedure was used for handling pressure-velocity coupling within each RK stage. Finally the developed OpenFOAM solvers were tested and compared computing the steady lid-driven cavity flow problem using uniformly spaced computational grids. On the other hand the unsteady flow past a cylinder at $\mathrm{Re}=140$ and the three-dimensional laminar flow past a sphere at $\operatorname{Re}=300$ were studied. In all the cases very satisfactory results were obtained.

The incompressible Navier–Stokes equations read

where u , denotes the velocity vector and $p=P / \rho$ is the pressure divided by the density. In order to discretize eqn (1) the computational domain is subdivided into $N_c$ finite volumes and each volume is bounded by an arbitrary number of cell faces $N_f$. Therefore the NavierStokes equations are re-written in integral form as regards each finite volume. The Gauss-Green divergence theorem is used to convert the volume integrals to surface integration over the mesh element boundary, thus the semi-discrete Navier-Stokes equations read

Extending the discretization process of eqn (2) to each discrete control volume we obtain an ordinary differential equation (ODE) system. The details about our solution algorithms adopted for time-integration are described in the following subsections.

In this research work we have used a class of explicit Runge-Kutta (ERK) methods for incompressible Navier-Stokes equations developed by Sanderse and Koren [3]. These methods are based on the following algorithm, expressed for the generic mesh element:

where $i$ is a generic Runge-Kutta substage.

Note that the terms $\tilde{\mathbf{U}}_i$ and $\tilde{p}_i$, appearing in eqn (3), are the entries of shifted vectors. Similarly also the coefficients $\tilde{a}_{i j}$ and $\tilde{c}_i$ are collected in the shifted coefficients matrix and vector [3]. The sequence reported in eqn (3) of first computing a tentative velocity, then the pressure, and finally correcting the tentative velocity is similar to fractional step methods. However all terms are handled explicitly (except the pressure) and consequently there is no splitting error involved as in fractional step approaches [3].

We have also to point out that in eqn (3) the projection of $\tilde{\mathbf{V}}$ onto the space of diver-gence-free velocity fields is known as a projected Runge-Kutta approach. It leads to the classical order of accuracy for the velocity up to the fourth order, as showed by Sanderse and Koren [3]. Differently the pressure is only first-order accurate unless additional order conditions are satisfied. In our work we have adopted third order accurate schemes for the velocity that lead also to a second-order accurate pressure (the Butcher tableau referred to our implementations are reported in [3]). It worth emphasizing that the projected nature of the previous schemes allows to avoid the use of Rhie-Chow correction even with co-located grids, thus the previous schemes result efficient and easy to implement.

This class of methods allows to solve in each stage a linear system with rank equal to the spatial degrees of freedom. Such technique, from an implementation point of view, is very attractive since each stage resembles a BDF1 scheme with a source term. Indeed, the discretized momentum equation showed in eqn (2) can be expressed by

where $\mathbf{R}_{\mathrm{p}}$ is the residuals vector, hence the solution at each stage of SDIRK scheme can be written as:

Furthermore in the next time-step $\mathbf{U}_{\mathrm{P}}$ is updated by:

where $s$ is the total number of RK stages.

It is also worth noting that in our work, for each implicit RK stage, we have used a PISO algorithm in order to handle the pressure-velocity coupling consisting of: (i) discretized momentum equation solution (predictor step); (ii) Poisson equation for pressure solution; (iii) cell-center velocity, face velocity and mass-flux correction. We have to note that face residual vector, $\mathbf{R}_f^{(i)}$ is obtained by the discretized form of a SDIRK stage, as in [4]. KazemiKamyab et al. [4] showed that this technique allows to avoid temporal order reduction suffered by several implementation techniques analyzed in the present literature. Lastly in this paper we have used stiffly accurate SDIRK schemes having the following feature: $a_{s j}=b_j$. Hence the solution of the last stage is equal to the solution of the new time level: $\mathbf{U}_{\mathrm{P}}^{(n+1)}=\mathbf{U}_{\mathrm{P}}^{(s)}$. This condition allows to avoid the use eqn (6). In this research work we have considered two different SDIRK techniques having, respectively 3, [5], and 5, [6], stages. In particular the first approach has a third order of accuracy (SDIRK 3-3) while the second the fourth order (SDIRK 4-5).

In this section we present solutions of different flow problems in order to assess the reliability and robustness as well as the accuracy of the presented RK methods for INS equations.

For the reported cases a preconditioned bi-conjugate gradient method (PBiCG) with DILU preconditioner was used for the velocity. Differently for the pressure a preconditioned conjugate gradient method (PCG) with a diagonal incomplete-Cholesky preconditioner was adopted. In particular linear system for pressure was solved using a local accuracy of $10^{-7}$ differently the system for the velocity was considered converged when the residuals reached $10^{-9}$. All the solutions have been computed on a small Linux Cluster with 6 Opteron nodes for a total of 96 cores operating at 2.4 GHz .

The first aim of this test is to demonstrate that the schemes, described in section 2, allow to march in time preserving the expected temporal accuracy. A square computational domain with: $\Omega=[-1,1]$ was used. The Dirichlet boundary conditions as well as initial condition were deduced by the exact solution. The computations were run with $\mathrm{Re}=100$ using a uniformly spaced quadrilateral grid having $128 \times 128$ cells. The convergence rate, for velocity and the pressure error, of the proposed schemes was evaluated in $L_2$ norm. The results are summarized in Fig. 1, where it is easy to observe as the theoretical order of accuracy is satisfactorily retained.

As first test case, we consider the two-dimensional lid-driven cavity flow. The computational domain for this problem is a box with edges of unit length. The top side of the cavity slides with a constant imposed velocity, while no-slip Dirichlet boundary conditions are fixed on the remaining sides. Here we have considered the cases with $\operatorname{Re}=10^3$ and $\operatorname{Re}=2.5 \cdot 10^3$ using uniformly spaced computational grids having quadrilateral type cells. Grids with $N_c$ from $64^2$ to $256^2$ were used. The previously cited flow problems are steady ones. They were considered as preliminary test cases in order to verify the proposed methods ability to converge to steady-solutions.

In Fig. 2 are showed the obtained results for $\operatorname{Re}=10^3$. The proposed solvers produce very close results; for the sake of clearness in Fig. 2 only the results obtained with SDIRK 4-5 are reported. With the $128^2$ cells grid the results show a very close match with the reference data of Erturk et al. [7], for the velocity, Fig. 2a. Very good agreement is obtained with Botella and Peyret [8], benchmark results for the pressure distribution.

A similar situation is evidenced for the $\mathrm{Re}=2.5 \cdot 10^3$ case. In particular in Fig. 3 are reported only the results related to our SDIRK 3-3 solver. However a $256^2$ cells grid is needed to reach a very close fitting with Erturk et al. [7] data, Fig. 2a. A final remark concerning to this benchmark case is concerned to the At size choice. The time-step size was fixed in order to obtain a maximum Courant number, $\mathrm{Co}_{\text {max }}$, lesser than 0.5 .

The first unsteady test case considered in this work is the flow past a circular cylinder at $\operatorname{Re}=140$. It is well known that this flow field exhibits the Karman vortex street, as represented in Fig. 4a. For this case we have built an O-type grid having $N_c=6 \cdot 10^4$ and the far field boundary was placed at 25 times the cylinder diameter. All the developed solvers were benchmarked on this flow case with different time-step sizes. Note that a report of the considered cases is disclosed in Table 1. More precisely in Table 1 the solutions are compared on the basis of the average drag coefficient: $\left\langle C_D\right\rangle$, force coefficients oscillation amplitudes $\left(\Delta C_D=\left(C_{D, \max }-C_{D, \min }\right) / 2\right.$ and $\left.\Delta C_L=\left(C_{L, \max }-C_{L, \min }\right) / 2\right)$ and the Strouhal number $S t=f D / U_{\infty, \text { where }} f$ the frequency of the vortex shedding. We have to put in evidence that for the sake of compactness only solutions with smaller time-step size are reported in Table 1.

The mean drag coefficient was determined in the range $C_d=1.31 \div 1.32$. The values are in very good agreement with DNS results of Inoue and Hatakeyama [9] and Muller [10] high-order finite difference data. As regards the force coefficients amplitudes we have obtained $\Delta C_L=0.47 \div 0.49$ and $\Delta C_D \cdot 10^2=2.16 \div 2.21$ which are in good agreement with reference solutions and experiments. The computed Strouhal number was 0.180 and also in this case the flow parameter exhibits a very good agreement with main literature references available in Table 1. Finally Fig. 4b shows that, for our finest solutions, the average pressure coefficient (on the cylinder wall) is in good agreement with literature, [9].

The flow past a sphere was also considered in this work in order to assess the solvers performance on three-dimensional configurations. The case with Re = 300, represented in Fig. 5, was selected as case study since it is markedly in the unsteady regime. The solutions were computed using a grid built on the half of the sphere adopting a symmetry plane with $N_c= 1.5 \cdot 10^6$. The far field boundary was placed at 15 times the sphere diameter. The flow features are clearly reflected on the time histories of the lift coefficient $C_L$ and drag coefficient $C_D$ reported in Fig. 6.

In Table 2 is showed a report of the performed computation for this test case together with literature references. It also put in evidence as our solutions show a good agreement with reference one and experimental data. Also for this benchmark we have used different $\Delta t$ to compute the flow field. However we have clearly noted the effect of the implicit approach. Indeed, the largest time-step size admitted by SDIRK solvers is significantly greater than ERK one. Lastly we have to remark that the SDIRK 3-3 and SDIRK 4-5 methods produce very similar results. On the other hand ERK technique underestimates the mean drag coefficient showing a very good agreement with other flow parameters.

In this work several OpenFOAM solvers for incompressible Navier-Stokes equations (NSE) based on high-order explicit and implicit RK schemes for time-integration were developed. The effectiveness of the proposed approaches was evaluated computing several standard benchkmarks for incompressible flows. The obtained results were very satisfactory. Note that the present solvers have also showed very good performance for transitional flows; these results are not presented for brevity. Future work will be devoted to turbulent flows computations and the obtained results will be presented in a future paper.