A PDF projection method: A pressure algorithm for stand-alone transported PDFs

Abstract

In this paper, a new formulation of the projection approach is introduced for stand-alone probability density function (PDF) methods. The method is suitable for applications in low-Mach number transient turbulent reacting flows. The method is based on a fractional step method in which first the advection–diffusion–reaction equations are modelled and solved within a particle-based PDF method to predict an intermediate velocity field. Then the mean velocity field is projected onto a space where the continuity for the mean velocity is satisfied. In this approach, a Poisson equation is solved on the Eulerian grid to obtain the mean pressure field. Then the mean pressure is interpolated at the location of each stochastic Lagrangian particle. The formulation of the Poisson equation avoids the time derivatives of the density (due to convection) as well as second-order spatial derivatives. This in turn eliminates the major sources of instability in the presence of stochastic noise that are inherent in particle-based PDF methods. The convergence of the algorithm (in the non-turbulent case) is investigated first by the method of manufactured solutions. Then the algorithm is applied to a one-dimensional turbulent premixed flame in order to assess the accuracy and convergence of the method in the case of turbulent combustion. As a part of this work, we also apply the algorithm to a more realistic flow, namely a transient turbulent reacting jet, in order to assess the performance of the method.

Keywords:

• transient jet
• particle method
• Monte Carlo technique
• projection method
• method of manufactured solutions

Acknowledgements

The authors thank the reviewers for their insightful comments, which improved the final draft of this manuscript.

Notes

1. In general (non-uniform grids), the commutation of the volume-averaging operator on spatial derivatives is second-order accurate [63] (i.e. $\overline{\nabla ·\mathbf{U}}=\nabla ·\bar{\mathbf{U}}+O\left({\mathbf{x}}^2\right)$ ), and since the discretisation scheme used in this work is second-order accurate, this error is tolerated.

2. The mean pressure is not shown by $\bar{p}$ since pressure is not calculated at particle level and it is not a particle property. The mean pressure is calculated on a grid and interpolated onto the particle position whenever it is required.

3. In a numerical implementation, the operators in Equation (35(35) $$\begin{array}{ccc}& & \nabla ·{\bar{\nu }}^{n+1}\nabla P^{n+\frac{1}{2}}=\frac{1}{\Delta t}\left(\nabla ·{\bar{U}}^{\Pi }-{\bar{S}}^{n+1}\right)+\nabla ·{\bar{\nu }}^n\nabla P^{n-\frac{1}{2}},\hfill \\ & & \text{or}\hfill \\ & & P^{n+\frac{1}{2}}={\left(\nabla ·{\bar{\nu }}^{n+\frac{1}{2}}\nabla \right)}^{-1}\left[\frac{1}{\Delta t}\left(\nabla ·{\bar{U}}^{\Pi }-{\bar{S}}^{n+1}\right)+\nabla ·{\bar{\nu }}^n\nabla P^{n-\frac{1}{2}}\right].\hfill \end{array}$$(35) ) need to be discretised. Specifically, on each grid point g, discretised analogue to divergence and gradient are defined as, respectively, ${\mathbf{D}}_g={(\nabla ·)}_g$ and ${\mathbf{G}}_g={\left(\nabla \right)}_g$. If a discrete Laplacian at a grid point g is defined by ${\mathbf{D}}_g\nu {\mathbf{G}}_g$, the projection is called ‘exact projection’. However, for practical reasons as discussed in [34,60] instead of forming the discrete Laplacian operators based on a discretised divergence and gradient, the continuous operator (∇ · ν∇) is discretised directly to form a discrete Laplacian ${\mathbf{L}}_g={(\nabla ·\nu \nabla )}_g$ where ${\mathbf{L}}_g\ne {\mathbf{D}}_g\nu {\mathbf{G}}_g$. Consequently, rather than the discretised algebraic system being a projection, it ‘approximates’ the continuous projection.

4. Generally, a second-order accurate central difference can also be used; however, we have observed that using a fourth-order accurate implementation improves the results (although it does not increase the overall order of accuracy of the algorithm).

5. A symbolic manipulation package (e.g. Maxima${}^{\mathrm{TM}}$, Mathematica^®, Sage${}^{\mathrm{TM}}$, etc.) can be used to obtain the analytical derivatives and integrations. In this work, Mathematica^® is used.

6. This configuration is a 1D problem; therefore, a 1D formulation would be sufficient for the calculations and discretisation in the lateral direction is not required. Since the aim of the calculations here is to assess the accuracy of the code, we use the same 2D code which is developed for a more general flow setup.

7. According to the velocity model used in this work, this can be achieved by replacing $(\frac{1}{2}+\frac{3}{4}{C}_0)$ with $\left(\frac{3}{4}{C}_0\right)$ in Equation (5(5) $$\mathrm{d}{U}_i^{*}=-\left({\nu }^{*}{\frac{\partial ⟨p⟩}{\partial x_i}}^{*}\right)\mathrm{d}t-\left(\frac{1}{2}+\frac{3}{4}{C}_0\right)\left(U_i^{*}-\widetilde{U_i}\right)\Omega \mathrm{d}t+\sqrt{C_{0}k\Omega }\mathrm{d}{W}_i,$$(5) ).

8. Since the mean progress variable is a monotonic function of distance in the streamwise direction x [30], the profiles can be shown against $\tilde{c}$ rather than x.

9. A classical gradient diffusion assumption would yield negative values of $\widetilde{c^{\text{'}}{u}^{\text{'}}}$.

Additional information

Funding

This work was supported in part by Underwriter Laboratories Inc. Northbrook, IL; and in part by the Deutsche Forschungsgemeinschaft within the Research Group FOR 1447.