Our research activities have a focus on the numerical modeling and simulation of fluid dynamic problems in nature and technology. More precisely, we are concerned with

  • computational fluid dynamics (CFD) simulations of
    • multiphase systems (sprays, gas-solid flows),
    • turbulent combustion,
  • low-dimensional stochastic modeling of combustion processes and multiphase flows,
  • model development for reactive and non-reactive fluid dynamic systems,
  • high-performance computing.

More detailed information on our current research is provided below.

Quadrature-Based Moment Methods (QBMMs) for Dispersed Multiphase Systems

Dispersed Multiphase Systems

Dispersed multiphase flows are dynamic systems of two or more phases, at least one of which consists of fine particles in any state of matter dispersed in a fluid. They are ubiquitous in nature and technology. A few examples of their occurence in nature are – from large to small scales – weather patterns, ocean dynamics and blood flow. On the other side, propulsion systems, power generation, heat exchange and manufacturing technologies highlight the importance of multiphase flows in technical and industrial applications. Numerical modeling of such systems is particularly challenging in terms of the mathematical description of the dispersed phase and the associated physical phenomena such as turbulence and various mechanisms of mass, momentum and energy exchange.

Modeling Approaches

The mathematical representation of the dispersed subsystem can be classified as microscopic, mesoscopic or macroscopic modeling approach:

  • Microscopic models aim to capture physics on the particle scale. Naturally, such methods provide the highest level of detail, though at the price of computational requirements that may be prohibitive in many physical applications.
  • Mesoscopic models aim at a statistical representation of the local population in terms of a number density function (NDF) of a vector of relevant particle properties (internal coordinates). Its evolution is governed by a transport equation, the so-called population balance equation (PBE). Due to its generally high dimensionality the numerical solution requires enormous computational resources.
  • Macroscopic models are often preferable since the quantities of interest are in many cases macroscopic quantities such as concentration, mass, momentum or kinetic energy.

Our research is primarily focused on a particular macroscopic modeling approach, namely the quadrature-based moment methods (QBMMs).

Quadrature-Based Moment Methods (QBMMs)

QBMMs are based on the macroscopic description in terms of a set of moments associated with the NDF. The governing moment transport equations can be derived from the underlying PBE. Presuming a single internal coordinate \(\xi\) in the presence of advection and diffusion with drift function \(\mathcal{A}(\xi)\) and diffusivity \(\mathcal{D}(\xi)\), the evolution of the NDF \(n\) over time \(t\) is governed by the PBE

            \(\frac{\partial n(\xi, t)}{\partial t} = - \frac{\partial }{\partial \xi} \left[ \mathcal{A}(\xi) \, n(\xi, t) \right] \; + \; \frac{\partial^2 }{\partial \xi^2} \left[ \mathcal{D}(\xi) \, n(\xi, t) \right]\).

Applying the definition of the \(k\)th moment \(m_k = \int_\Omega \xi^k n(\xi) \mathrm{d} \xi\) to the PBE yields the \(k\)th moment transport equation

            \(\frac{\mathrm{d} m_k(t)}{\mathrm{d} t} = - k \int_\Omega \xi^{k-1} \, \mathcal{A}(\xi) \, n(\xi, t) \mathrm{d} \xi \; + \; k(k-1) \int_\Omega \xi^{k-2} \, \mathcal{D}(\xi) \, n(\xi, t) \mathrm{d} \xi\).

Thus, without knowledge of n , the moment equations are unclosed. This closure problem is solved by QBMMs employing quadrature rules, i.e.

            \(\frac{\mathrm{d} m_k(t)}{\mathrm{d} t} \approx - k \sum\limits_j w_j \, \xi_j^{k-1} \, \mathcal{A}(\xi_j) \; + \; k(k-1) \sum\limits_j w_j \, \xi_j^{k-2} \, \mathcal{D}(\xi_j)\).

If the quadrature nodes \(\xi_j\) and weigths \(w_j\) are those of the Gaussian quadrature associated with the known moment sequence this method is referred to as the quadrature method of moments (QMOM) [1], which is the basis of various derived methods, see [2]. All have in common that the quadrature formula can be computed solely from moment information. We extended the QMOM making use of the anti-Gaussian quadrature formulas due to [3].

Results: Droplet Breakup

We formulated moment transport equations based on a modified version of the widely used Reitz-Diwakar model for droplet breakup [4] and applied the QMOM for closure. Figure 2 shows a comparison with highly-resolved Monte Carlo simulations of the PBE in terms of the total number concentration \(m_0\) as well as the Sauter mean diameter \(m_3/m_2\). Using the QMOM with more than two quadrature nodes, a high degree of accuracy can be achieved with significantly lower computational costs.

Results: Turbulence-Induced Diffusion

We applied our proposed extension of the QMOM, the Gauss/anti-Gauss QMOM (GaG-QMOM), to a spatially homogeneous turbulent system. The strongly nonlinear turbulence-induced terms in the PBE present a particular challenge to QBMMs. As illustrated in Figure 3 and Figure 4, the GaG-QMOM results in significantly enhanced accuracy in all cases.


[1]R. McGraw. Description of aerosol dynamics by the quadrature method of moments. Aerosol Sci. Tech., 27(2):255–265, 1997.
[2]D. L. Marchisio and R. O. Fox. Computational Models for Polydisperse Particulate and Multiphase Systems. Cambridge University Press, 2013.
[3]D. Laurie. Anti-Gaussian quadrature formulas. Math. Comput., 65(214):739–747, 1996.
[4]R. D. Reitz and R. Diwakar. Effect of drop breakup on fuel sprays. SAE Tech. Pap., 860469,1986.
Efficient modeling and simulation of reactive gas dynamics on adaptive Cartesian grids

Overview and aim of the research

In cooperation with the Geophysical Fluid Dynamics group of the FU Berlin, we are investigating the theoretical framework of Shockless Explosion Combustion (SEC), which represents a completely new operating concept for gas turbines. Core of the SEC process is the controlled quasi-homogeneous auto-ignition of an air-fuel mixture, which should lead to a dynamic approximation of an energetically more favorable constant volume combustion (CVC) [2]. The numerical investigation of the flows occuring here involves very small time and length scales, which also requires significant effort from the computational fluid dynamics (CFD) solver. Therefore, we are developing our own CFD solver and pay special attention to efficient numerical methods and parallelization of the simulations. In the following, we give a rough overview of the numerical methods we use and focus especially on methods for adaptive mesh refinement (AMR).

Numerical Methods

The general basis for the numerical solution of problems we discuss here are the multi-species reactive Euler Equations. These are given in their three-dimensional conservative Form as
            \frac{\partial \mathbf{U}}{\partial t}
            + \frac{\partial \mathbf{F}(\mathbf{U})}{\partial x}
            + \frac{\partial \mathbf{G}(\mathbf{U})}{\partial y}
            + \frac{\partial \mathbf{H}(\mathbf{U})}{\partial z}
            = \mathbf{S}(\mathbf{U}) \,, \label{eq:eulerEquation} \\
             \mathbf{U} = \begin{bmatrix}
                \rho \\ \rho u \\ \rho v \\ \rho w \\ \rho E \\ \rho Y_s
                \end{bmatrix} ,\,
             \mathbf{F} = \begin{bmatrix}
                \rho u \\ \rho u^2+p \\ \rho uv \\ \rho uw \\ u(\rho E+p) \\ \rho u Y_s
            \end{bmatrix} ,\,
             \mathbf{G} = \begin{bmatrix}
                \rho v \\ \rho uv \\ \rho v^2+p \\ \rho vw \\ v(\rho E+p) \\ \rho v Y_s
            \end{bmatrix} ,\,
             \mathbf{H} = \begin{bmatrix}
                \rho w\\ \rho uw \\ \rho vw \\ \rho w^2+p \\ w(\rho E+p) \\ \rho w Y_s
            \end{bmatrix} ,\, \label{eq:eulerComponents}
where \(\rho\) is the density, \(\mathbf{u}=(u,v,w)\) the velocity, \(E=e+||\mathbf{u}||^2/2\) the total energy, \(e\) the internal energy, \(p\) the pressure and \(Y_s,\,s=1,\ldots,N_s-1\) are the mass fractions of the reactive species that are normalized to sum to unity. The equations above are closed by the caloric perfect gas approximation \(p=\rho(\gamma-1)e\), where \(\gamma=c_p/c_v\) is the ratio of specific heat. The term \(\mathbf{S}(\mathbf{U})\) is a placeholder for any possible source term that accounts, e.g. external forces, chemical reactions or any injection of mass, momentum, energy or species. To compute the numerical solution to the above equations, we utilize fully conservative finite volume schemes discretised on block-structured Cartesian grids. Due to the usage of the external C++ library AMReX [5] in our CFD solver, we can locally refine the grid in regions of interest through adaptive mesh refinement (AMR) techniques [1]. To ensure stability of any small or irregular cells that occur due embedding complex geometries on the structured grids, we use a conservative cut-cell method [3].

Adaptive Mesh Refinement (AMR)

Figure 1 a) shows a simplified sketch of an adaptive Cartesian grid with three levels \((l=0,1,2)\) in total, where the computational domain on each AMR level is decomposed into a union of rectangular domains. Bold lines represent grid boundaries and all cells of each AMR level \(l_i\) are refined with refinement ratio \(r\) in comparison to level \(l_{i-1}\). For efficient time integration on adaptive grids, we can use the so called subcycling-in-time approach [5], where each AMR level is advanced with it's own time step size \(\Delta t / r^l\) and synchronisation take place when level \(l_i\) and \(l_{i-1}\) have reached the same time point. This approach is shown in Figure 1 b) as an example for three levels. Two problems arise with this approach: (i) Coarse data on level \(l_{i-1}\) is not necessarily equal to the average of finer data on level \(l_i\). (ii) Due to inconsistent fluxes at the coarse-fine interface conservation is violated. Therefore, in each synchronisation step we must (i) fill the coarse grid data with averaged data from the overlaying fine grid and (ii) do a flux correction on the coarse-fine interfaces.


Here we present some selected results from simulations to study the diffraction of shock waves through non-quiscient media [4]. For this simulations we consider a one-dimensional tube traversed by a given mass flow, where after a certain time a shock is implemented. The tube is attached by a twodimensional plenum, where the shock will be diffracting. Here we limit the simulations to two dimensions, no species and add an additional geometrical source term:
      \mathbf{S}(\mathbf{U}) = - \frac{1}{y} \begin{bmatrix}
          \rho v ,& \rho uv ,& \rho v^2 ,& v(\rho E+p)
which accounts for cylindrical symmetry around the \(x\)-axis in the plenum. A sketch of our simulation domain with boundary conditions (BC’s) is shown in Figure 2, where the inset figure shows details of the embedded plenum in the structured adaptive grid with three levels in total.

Figure 3 shows contours of the non-dimensional velocity together with the adaptive grid. Initial Flow conditions are \(M_0=2.62\) and \(M_b=0.38\) for the shock Mach number and background Mach number, respectively. Comparing the two subfigures, the adaptive grid allows a high resolution of high gradients and discontinuities. In both subfigures it can be clearly seen that only the regions of the grid around embedded boundaries, the vortices formed by Kevin-Helmholtz instabilities in the shear layer of the background jet and those of the shock front have been refined. Whereas the other regions of the flow are only mapped on the coarsest grid.


[1]M. J. Berger and J. Oliger. Adaptive mesh refinement for hyperbolic partial differential equations. Journal of Computational Physics, 53(3):484-512, 1984.
[2]B. C. Bobusch, P. Berndt, C. O. Paschereit, and R. Klein. Shockless Explosion Combustion: An Innovative Way of Efficient Constant Volume Combustion in Gas Turbines. Combustion Science and Technology, 186(10-11):1680–1689, 2014.
[3]N. Gokhale, N. Nikiforakis, and R. Klein. A dimensionally split Cartesian cut cell method for hyperbolic conservation laws. Journal of Computational Physics, 364:186–208, 2018.
[4]B. S. Thethy, M. R. Haghdoost, B. S. Rhiannon Kirby, M. Nadolski, C. Zenker, M. Oevermann, R. Klein, K. Oberleithner, and D. Edgington-Mitchell. Diffraction of shock waves through a non-quiescent medium. Journal of Fluid Mechanics, 944:A39, 2022.
[5]W. Zhang, A. Almgren, V. Beckner, and J. B. et al. AMReX: a framework for block-structured adaptive mesh refinement. Journal of Open Source Software, 4(37):1370, 2019.