A pressure-based solver for low-Mach number flow using a discontinuous Galerkin method

doi:10.1016/j.jcp.2020.109877

Journal of Computational Physics

Volume 425, 15 January 2021, 109877

https://doi.org/10.1016/j.jcp.2020.109877 Get rights and content

Highlights

•
We present a numerical method for low-Mach number flows with a time-splitting scheme.
•
A very simple change to the time discretization of the enthalpy guarantees stability and increases the order of accuracy.
•
A set of variable-property manufactured solutions demonstrate second-order accuracy in time, without the need for iterations.
•
The spatial discretization is a discontinuous Galerkin method, both equal-order and mixed-order.

Abstract

Over the past two decades, there has been much development in discontinuous Galerkin methods for incompressible flows and for compressible flows with a positive Mach number, but almost no attention has been paid to variable-density flows at low speeds. This paper presents a pressure-based discontinuous Galerkin method for flow in the low-Mach number limit. We use a variable-density pressure correction method, which is simplified by solving for the mass flux instead of the velocity. The fluid properties do not depend significantly on the pressure, but may vary strongly in space and time as a function of the temperature.

We pay particular attention to the temporal discretization of the enthalpy equation, and show that the specific enthalpy needs to be ‘offset’ with a constant in order for the temporal finite difference method to be stable. We also show how one can solve for the specific enthalpy from the conservative enthalpy transport equation without needing a predictor step for the density. These findings do not depend on the spatial discretization.

A series of manufactured solutions with variable fluid properties demonstrate full second-order temporal accuracy, without iterating the transport equations within a time step. We also simulate a Von Kármán vortex street in the wake of a heated circular cylinder, and show good agreement between our numerical results and experimental data.

Introduction

Several low-speed flows of practical importance are compressible, that is, the velocity is not divergence-free. This can occur due to mixing, or due to a temperature-dependent density near a heat source. An example is heat transfer in low-Mach number flows of supercritical fluids, where all fluid properties vary strongly with the temperature, but do not depend significantly on the pressure. Most flow solvers use either a pressure-based approach and assume a divergence-free velocity field, or a fully compressible (density-based) formulation. Neither of these methods is directly applicable to compressible flows in the low-Mach number limit.

Density-based solvers can be used to simulate zero-Mach flows by approximating the flow with a low, non-zero Mach number (e.g., [46], [11]). This has often been used for heat transfer in supercritical fluids at low speeds (e.g., [14], [4]). This is expensive for several reasons. First, the temporal discretization needs to resolve acoustic effects, and the resulting linear systems tend to be very stiff. Second, the system of transport equations is solved in a coupled way, which is more expensive than using a time-splitting method, though the performance may be improved with suitable preconditioning [31]. Finally, the fluid properties are evaluated as a function of two thermodynamic variables (usually the density and the volumetric enthalpy), so that a spline interpolation costs far more memory, thus complicating massively parallel calculations [14].

There is also substantial experience with discontinuous Galerkin (DG) discretizations for incompressible flows. These are either based on the introduction of artificial compressibility (e.g. [12]), or they solve for the pressure (e.g. [41], [8], [37], [35]). The artificial compressibility method can be more than second-order accurate in time, though it requires the system of transport equations to be solved in a coupled manner (e.g., [6], [5]). By choosing entropy variables as the unknowns, the DG method can also be formulated in a general way for both compressible and incompressible flows, at the cost of great complexity (e.g., [34]). There is, however, almost no literature on solving the low-Mach number equations with a pressure-based discontinuous Galerkin method.

The only previous work of which we are aware is by Klein et al. [24], [24], [25], who used a SIMPLE scheme to march the transport equations forward in time, iterating the equations within each time step. This required under-relaxation in order for the iteration to converge. They solved for the velocity, so that a predictor for the density is needed in the temporal derivative of the momentum equation.

We avoid this by solving for the mass flux rather than the velocity. Another advantage of this approach is that the divergence term in the continuity equation does not have to be weighed by the density, so that the divergence matrix does not depend on the density. This makes the transport equations less tightly coupled, and it simplifies the pressure correction method, because the pressure matrix is constant for each time step.

The temporal discretization of the enthalpy equation requires special care. The convective term in the primitive enthalpy equation cannot be upwinded with a finite element method if the velocity is not divergence-free. We therefore discretize the enthalpy equation in conservative form, although we solve for the specific enthalpy, which is a primitive variable. This causes a complication in the temporal discretization, which features the unknown density at a new time step. A similar issue occurs in multispecies transport, and Najm et al. [28], [28] devised a wide-spread two-step iterative method to stabilize the temporal scheme. This has subsequently been adapted to handle the strong property variations in supercritical fluids [29].

We present an alternative method that does not use any predictor solves or iterations to handle the unknown density at a new time step. We show that the error in our approximation can be made negligible compared to the error in the finite difference scheme, and that the method can be made stable by offsetting the specific enthalpy with a constant.

The paper is structured as follows. Section 2 defines the mathematical problem. Section 3 and 4 explain the spatial discretization and the time-splitting scheme, paying particular attention to the polynomial solution spaces (Section 3.1) and the discretization of the viscous stress (Section 3.2.1), which differ slightly from previous literature. Section 5 and 6 treat the density predictor in the enthalpy equation in great detail, and illustrate the results with numerical examples; these sections do not depend on the spatial discretization, and they can be read independently. Section 7 verifies the numerical method and the implementation for the coupled transport equations with a series of progressively more challenging manufactured solutions, and validates our approach by reproducing experimental data for flow past a heated circular cylinder. Finally, we draw our conclusions based on the results in Section 8.

Section snippets

Governing equations

The transport equations in the low-Mach number limit are $\frac{\partial ρ}{\partial t} + \partial_{k} m_{k} = 0,$ $\frac{\partial m_{i}}{\partial t} + \partial_{k} (u_{k} m_{i}) = \partial_{k} τ_{i k} - \partial_{i} p + F_{i},$ $\frac{\partial ρ h}{\partial t} + \partial_{k} (m_{k} h) = - \partial_{k} q_{k} + Q .$ Here t is the time, ρ is the density, u is the velocity, $m : = ρ u$ is the mass flux, p is the pressure, h is the specific enthalpy, and F and Q are known external sources. Assuming a Newtonian fluid, the stress tensor is $τ_{i j} = μ (\partial_{i} u_{j} + \partial_{j} u_{i} - \frac{2}{3} \partial_{k} u_{k} δ_{i j}) .$ The heat flux is $q_{i} = - k \partial_{i} T = - \frac{k}{c_{p}} \partial_{i} h,$ where T is the temperature. k is the thermal conductivity, and $c_{p}$ is the specific heat capacity. The

Spatial discretization

This section details the spatial discretization of the enthalpy and mass flux transport equations on a domain Ω with outward normal n and space vector $r = [x, y] \in Ω$ . The domain is meshed into a set of elements $T$ . We denote by $F^{i}$ , $F^{D}$ , and $F^{N}$ the sets of internal faces, Dirichlet boundary faces, and Neumann boundary faces. The set of faces of an element $T \in T$ is $F_{T}$ , and $T_{F}$ is the set of neighbors of face F. Each internal face $F \in F^{i}$ has a normal vector $n^{F}$ with an arbitrary but fixed direction. The jump

Pressure correction method

We discretize the temporal derivatives in both the enthalpy and momentum transport equations with standard backward-difference formulae (BDF), thereby following previous DG literature (e.g., [41], [24], [26]). For the mass flux this is straightforward: for a constant time step size δt, $\frac{\partial m}{\partial t} \approx \frac{γ_{0}}{δ t} m^{n} + \sum_{i = 1}^{q} \frac{γ_{i}}{δ t} m^{n - i},$ where $m^{n}$ is the mass flux at time step n. The weights ${γ_{i}}_{i = 0}^{q}$ are listed in Table 1. The temporal discretization of the enthalpy equation requires more care, as we will explain later in

Linearizing the temporal derivative of the enthalpy

Solving for a primitive variable (h) with the enthalpy equation in conservative form slightly complicates the temporal derivative, because it is weighed by the temperature-dependent density. To study the stability and convergence of the time stepping scheme, we consider the space-independent enthalpy equation: $\frac{d (ρ h)}{d t} = - λ h + Q,$ where λ is a constant, and $Q = Q (t)$ . Using an implicit finite difference scheme, the enthalpy and the corresponding density can be estimated at a time step n by $\frac{γ_{0}}{δ t} {(ρ h)}^{n} + \sum_{i = 1}^{q} γ$

Test case for the space-independent Enthalpy Equation¹

Before solving the full system of transport equations, we clarify the theoretical results for the space-independent enthalpy equation in Section 5 with a numerical example that is based on a manufactured solution. Omitting the units of measurement, the exact temperature is $T^{ex} (t) = 0.5 + 0.1 \sin (2 π t)$ with $0 \leq t \leq 1$ . The equation of state is $ρ = ρ_{0} T + ρ_{1} (1 - T),$ and the specific heat capacity is kept constant, so that $h = c_{p} T - h_{0} .$ The required source term $Q (t)$ follows from Eq. (33). For the numerical tests we let $ρ_{0}$

Test cases for the coupled transport equations

The numerical method for the coupled transport equations and our computer implementation are verified and validated with the following test cases. We first investigate a series of manufactured solutions for the system of transport equations, ranging from constant properties (Section 7.2), to variable transport properties (Section 7.3.1), and a variable density (Section 7.3.2). As in Section 6, we leave out the units of the variables for these manufactured solutions. Section 7.4 then shows

Discussion and conclusions

We have presented a discontinuous Galerkin method for solving the system of transport equations in their conservative form for low-Mach number flows, while taking the mass flux m, the pressure p, and the specific enthalpy h as the unknowns.

Since the density is a function of the specific enthalpy, the temporal finite difference scheme requires that the volumetric enthalpy (ρh) be linearized in h, and we have analyzed two methods for doing this, both of which need a predictor for the enthalpy at

CRediT authorship contribution statement

Aldo Hennink: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing - original draft, Writing - review & editing. Marco Tiberga: Methodology, Software, Writing - review & editing. Danny Lathouwers: Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

We are grateful to the sCO₂-HeRo project that has received funding from the European research and training program 2014-2018 (grant agreement ID 662116) for funding the first author.

References (49)

F. Bassi et al.
An artificial compressibility flux for the discontinuous Galerkin solution of the incompressible Navier-Stokes equations
J. Comput. Phys.
(2006)
F. Bassi et al.
Linearly implicit Rosenbrock-type Runge–Kutta schemes applied to the Discontinuous Galerkin solution of compressible and incompressible unsteady flows
Comput. Fluids
(2015)
Lorenzo Botti et al.
A pressure-correction scheme for convection-dominated incompressible flows with discontinuous velocity and continuous pressure
J. Comput. Phys.
(2011)
Andrea Crivellini et al.
A Spalart-Allmaras turbulence model implementation in a discontinuous Galerkin solver for incompressible flows
J. Comput. Phys.
(2013)
J.-L. Guermond et al.
A splitting method for incompressible flows with variable density based on a pressure Poisson equation
J. Comput. Phys.
(may 2009)
C.M. Klaij et al.
Space-time discontinuous Galerkin method for the compressible Navier-Stokes equations
J. Comput. Phys.
(2006)
Benjamin Krank et al.
A high-order semi-explicit discontinuous Galerkin solver for 3d incompressible flow with application to dns and les of turbulent channel flow
J. Comput. Phys.
(2017)
G. Mengaldo et al.
Dealiasing techniques for high-order spectral element methods on regular and irregular grids
J. Comput. Phys.
(Oct. 2015)
Habib N. Najm et al.
A semi-implicit numerical scheme for reacting flow, I: stiff chemistry
J. Comput. Phys.
(Jul. 1998)
F. Nicoud
Conservative high-order finite-difference schemes for low-Mach number flows
J. Comput. Phys.
(2000)

L. Pesch et al.

A discontinuous Galerkin finite element discretization of the Euler equations for compressible and incompressible fluids

J. Comput. Phys.

(2008)

Marian Piatkowski et al.

A stable and high-order accurate discontinuous Galerkin based splitting method for the incompressible Navier-Stokes equations

J. Comput. Phys.

(2018)

Sander Rhebergen et al.

A space-time discontinuous Galerkin method for the incompressible Navier-Stokes equations

J. Comput. Phys.

(2013)

Khosro Shahbazi et al.

A high-order discontinuous Galerkin method for the unsteady incompressible Navier-Stokes equations

J. Comput. Phys.

(2007)

Lee Shunn et al.

Verification of variable-density flow solvers using manufactured solutions

J. Comput. Phys.

(2012)

O. Axelsson et al.

Numerical solution of the time-dependent Navier-Stokes equation for variable density-variable viscosity. Part I

Math. Model. Anal.

(2015)

Satish Balay et al.

Efficient management of parallelism in object oriented numerical software libraries

Satish Balay et al.

PETSc users manual

(2018)

R. Barney et al.

Fully-implicit, high-order, reconstructed discontinuous Galerkin method for supercritical fluid flows

Ian H. Bell et al.

Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library CoolProp

Ind. Eng. Chem. Res.

(Feb. 2014)

K. Andrew Cliffe et al.

Adaptive discontinuous Galerkin methods for eigenvalue problems arising in incompressible fluid flows

SIAM J. Sci. Comput.

(Jan. 2010)

Bernardo Cockburn et al.

An equal-order DG method for the incompressible Navier-Stokes equations

J. Sci. Comput.

(2009)

S. Scott Collis

Discontinuous Galerkin methods for turbulence simulation

C. Ross Ethier et al.

Exact fully 3D Navier-Stokes Solutions for Benchmarking

Int. J. Numer. Methods Fluids

(March 1994)

Cited by (13)

A numerical benchmark for modelling phase change in molten salt reactors
2023, Annals of Nuclear Energy
The design of a molten salt reactor is largely based on CFD simulations. Phase change plays an important role in the safety of the reactor, but numerical modelling of phase change is particularly challenging. Therefore, the knowledge of the margin of error of CFD simulations involving phase change is very important. Relevant experimental validation data is lacking. For this reason, a numerical benchmark designed after the freeze valve is proposed. The benchmark consists of five stages, where with each step more complexity is added. The step-wise addition of complexity allows for pinpointing potential sources of discrepancy. Results were obtained with three different codes: STAR-CCM+, OpenFOAM, and DGFlows. The results were found to be largely consistent between the codes, however the addition of conjugate heat transfer introduced some discrepancies. These results indicate that careful consideration is needed when coupling conjugate heat transfer solvers with solid–liquid phase change models.
Code-verification techniques for the method-of-moments implementation of the combined-field integral equation
2023, Journal of Computational Physics
Code verification plays an important role in establishing the credibility of computational simulations by assessing the correctness of the implementation of the underlying numerical methods. In computational electromagnetics, the numerical solution to integral equations incurs multiple interacting sources of numerical error, as well as other challenges, which render traditional code-verification approaches ineffective. In this paper, we provide approaches to separately measure the numerical errors arising from these different error sources for the method-of-moments implementation of the combined-field integral equation. We demonstrate the effectiveness of these approaches for cases with and without coding errors.
Code-verification techniques for the method-of-moments implementation of the magnetic-field integral equation
2023, Journal of Computational Physics
For computational physics simulations, code verification plays a major role in establishing the credibility of the results by assessing the correctness of the implementation of the underlying numerical methods. In computational electromagnetics, surface integral equations, such as the method-of-moments implementation of the magnetic-field integral equation, are frequently used to solve Maxwell's equations on the surfaces of electromagnetic scatterers. These electromagnetic surface integral equations yield many code-verification challenges due to the various sources of numerical error and their possible interactions. In this paper, we provide approaches to separately measure the numerical errors arising from these different error sources. We demonstrate the effectiveness of these approaches for cases with and without coding errors.
Analysis of the Molten Salt Fast Reactor using reduced-order models
2021, Progress in Nuclear Energy
In this paper, we present a reduced-order modeling approach to study the Molten Salt Fast Reactor (MSFR). Our approach is nonintrusive and based on the proper orthogonal decomposition method. We include adaptivity in selecting the sampling points both in time and parameter space. Steady-state and transient analysis were both performed using the developed models. In the steady-state analysis, we capture the effect of 30 model parameters on the spatial distributions of fission power and temperature, and on the multiplication factor. The dimensionality of the fission power was reduced from the 104288 nominal dimensions in the physical space to 10 dimensions in the reduced space, whereas the temperature was reduced from 220972 dimensions to 3. The reduced model was then used for uncertainty and sensitivity study of the maximum temperature in the reactor and the multiplication factor. In the transient analysis, the reduced model captured the effect of perturbations in the flow rate of salt in the intermediate circuit on the fission power density and temperature. The reduced models were successfully tested on a set of points that were not part of the snapshots used during the construction stage.
Preliminary uncertainty and sensitivity analysis of the Molten Salt Fast Reactor steady-state using a Polynomial Chaos Expansion method
2021, Annals of Nuclear Energy
In this work, we present the results of a preliminary uncertainty quantification and sensitivity analysis study of the Molten Salt Fast Reactor (MSFR) behavior at steady-state performed by applying a non-intrusive Polynomial Chaos Expansion (PCE) approach. An in-house high-fidelity multi-physics simulation tool is used as reactor reference model. Considering several thermal-hydraulics and neutronics parameters as stochastic inputs, with a limited number of samples we build a PCE meta-model able to reproduce he reactor response in terms of effective multiplication factor, maximum, minimum, and average salt temperatures, and complete salt temperature distribution. The probability density functions of the responses are constructed and analyzed, highlighting strengths and issues of the current MSFR design. The sensitivity study highlights the relative importance of each input parameter, thus providing useful indications for future research efforts. The analysis on the whole temperature field shows that the heat exchanger can be a critical component, so its design requires particular care.
1 Incremental Pressure Correction Method for Subsonic Compressible Flows
2024, SSRN

View all citing articles on Scopus

View full text

A pressure-based solver for low-Mach number flow using a discontinuous Galerkin method

Highlights

Abstract

Introduction

Section snippets

Governing equations

Spatial discretization

Pressure correction method

Linearizing the temporal derivative of the enthalpy

Test case for the space-independent Enthalpy Equation1

Test cases for the coupled transport equations

Discussion and conclusions

CRediT authorship contribution statement

Declaration of Competing Interest

Acknowledgements

J. Comput. Phys.

Comput. Fluids

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

Numerical solution of the time-dependent Navier-Stokes equation for variable density-variable viscosity. Part I

Math. Model. Anal.

Efficient management of parallelism in object oriented numerical software libraries

PETSc users manual

Fully-implicit, high-order, reconstructed discontinuous Galerkin method for supercritical fluid flows

Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library CoolProp

Ind. Eng. Chem. Res.

Adaptive discontinuous Galerkin methods for eigenvalue problems arising in incompressible fluid flows

SIAM J. Sci. Comput.

An equal-order DG method for the incompressible Navier-Stokes equations

J. Sci. Comput.

Discontinuous Galerkin methods for turbulence simulation

Exact fully 3D Navier-Stokes Solutions for Benchmarking

Int. J. Numer. Methods Fluids

Test case for the space-independent Enthalpy Equation¹