key: cord-0045962-oidnevdg
authors: Parekh, Jigar; Verstappen, Roel
title: Intrusive Polynomial Chaos for CFD Using OpenFOAM
date: 2020-05-25
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50436-6_50
sha: cdcf8d6076c5461dbe759460d4628da1587be4ee
doc_id: 45962
cord_uid: oidnevdg

We present the formulation and implementation of a stochastic Computational Fluid Dynamics (CFD) solver based on the widely used finite volume library - OpenFOAM. The solver employs Generalized Polynomial Chaos (gPC) expansion to (a) quantify the uncertainties associated with the fluid flow simulations, and (b) study the non-linear propagation of these uncertainties. The aim is to accurately estimate the uncertainty in the result of a CFD simulation at a lower computational cost than the standard Monte Carlo (MC) method. The gPC approach is based on the spectral decomposition of the random variables in terms of basis polynomials containing randomness and the unknown deterministic expansion coefficients. As opposed to the mostly used non-intrusive approach, in this work, we use the intrusive variant of the gPC method in the sense that the deterministic equations are modified to directly solve for the (coupled) expansion coefficients. To this end, we have tested the intrusive gPC implementation for both the laminar and the turbulent flow problems in CFD. The results are in accordance with the analytical and the non-intrusive approaches. The stochastic solver thus developed, can serve as an alternative to perform uncertainty quantification, especially when the non-intrusive methods are significantly expensive, which is mostly true for a lot of stochastic CFD problems.

In simulating a physical system with a model, uncertainties may arise from various sources [15] , namely, initial and boundary conditions, material properties, model parameters, etc. These uncertainties may involve significant randomness or may only be approximately known. In order to enhance the predictive reliability, it is therefore important to quantify the associated uncertainties and study its non-linear propagation especially in CFD simulations.

In order to reflect the uncertainty in the numerical solution, we need efficient Uncertainty Quantification (UQ) methods. Broadly there exists two classes of UQ methods, the intrusive method, where the original deterministic model is replaced by its stochastic representation, and the non-intrusive method, where the original model itself is used without any modifications [11, 14] . Monte Carlo (MC) sampling is one of the simplest non-intrusive approaches. However, due to its requirement of a large number of samples, MC method is computationally expensive for application in CFD. As an alternative, we can use Generalized Polynomial Chaos (gPC) representations which has been proven to be much cheaper than MC [6, 15] . This approach is based on the spectral decomposition of the random variables in terms of basis polynomials containing randomness and the unknown deterministic expansion coefficients. In this paper, we focus mainly on the Intrusive Polynomial Chaos (IPC) method, where a reformulation of the original model is performed resulting in governing equations for the expansion coefficients of the model output [20] .

As the model code, we use OpenFOAM [2] , which is a C++ toolbox to develop numerical solvers, and pre-/post-processing utilities to solve continuum mechanics problems including CFD. OpenFOAM (a) is a highly templated code, enabling the users to customize the default libraries as needed for their applications, and, (b) gives access to most of the tensor operations (divergence, gradient, laplacian etc.) directly at the top-level code. This avails enough flexibility to implement the IPC framework for uncertainty quantification in CFD. To obtain the inner products of polynomials we use a python library called chaospy [3] , as a pre-processing step to the actual stochastic simulation.

First, the idea of generalized polynomial chaos is presented with a focus on the intrusive variant. A generic differential equation is used to explain the steps involved in IPC, leading to a simple expression for the mean and variance as a function of the expansion coefficients. Next, we present the set of deterministic governing equations followed by its stochastic formulation using IPC. In particular, a Large Eddy Simulation (LES) method is used to model turbulence, which includes an uncertain model parameter. Thereafter, we discuss the algorithm and implementation steps required for the new stochastic solver in OpenFOAM. The stochastic version of the Navier-Stokes equations has a similar structure to the original system. This allows reusing the existing deterministic solver with minimal changes necessary.

The stochastic solver developed so far is tested for various standard CFD problems. Here we present two cases, the plane Poiseuille flow with uncertain kinematic viscosity, and the turbulent channel flow with uncertain LES model parameter. The results are found to be in accordance with the non-intrusive gPC method.

The Generalized Polynomial Chaos approach is based on the spectral decomposition of the random variable(s) f , in terms of basis polynomials containing randomness ψ i (known a priori) and the unknown deterministic expansion coefficients f i , as f (x, q) = Non-intrusive Polynomial Chaos (NIPC). In IPC, a reformulation of the original model is performed resulting in governing equations for the PC mode strengths of the model output, while in NIPC, these coefficients are approximated using quadrature for numerical evaluation of the projection integrals.

The level of accuracy of these methods can be associated with the degree of gPC. To attain the same level of accuracy, particularly for a higher dimensional random space, IPC requires the solution of a much fewer number of equations that needed for NIPC. Moreover, for such a random space, the aliasing error resulting from the approximation of the exact gPC expansion in the NIPC method can become significant. This suggests that, for a multi-dimensional problem, the IPC method can deliver more accurate solutions at a much lower computational cost than the NIPC method [18] .

Since the current work is based on IPC, we would introduce here the important features of the intrusive variant and we refer to the literature [5, 16, 19] for more details about NIPC and gPC in general.

In order to demonstrate the application of IPC, we first consider a general stochastic differential equation

where L is usually a nonlinear differential operator consisting of space and/or time derivatives, v(x, t, ω) is the solution and S(x, t, ω) is the source term. The random event ω represents the uncertainty in the system, introduced via uncertain parameters, the operator, source term, initial/boundary conditions, etc. The complete probability space is given by (Ω, A, P), where Ω is the sample space such that ω ∈ Ω, A ⊂ 2 Ω is the σ-algebra on Ω and P : A → [0, 1] is the probability measure on (Ω, A). We now employ the Galerkin polynomial chaos method, which is an IPC method for the propagation of uncertainty [16] . It provides the spectral representation of the stochastic solution and results into higher order approximations of the mean and variance. Galerkin polynomial chaos method is a non-statistical method where the uncertain parameter(s) and the solution become random variables. These random variables are approximated using the polynomial chaos (polynomial of random variables) as follow [5] v(

(2)

It is worth noting that the expansion (2) is indeed the decomposition of a random variable into a deterministic component, the expansion coefficients v i (x, t) and a stochastic component, the random basis functions (polynomial chaoses) ψ i (ξ(ω)). Here, ξ(ω) is the vector of d independent random variables {ξ 1 , ..., ξ d }, corresponding to d uncertain parameters. Based on the dimension of ξ (which here is d) and the highest order n of the polynomials {ψ i }, the infinite summation has been truncated to P + 1 = (d + n)!/(d! n!) terms.

An important property of the basis {ψ i } is their orthogonality with respect to the probability density function (PDF) of the uncertain parameters, ψ i ψ j = ψ 2 i δ ij . Here, δ ij is the Kronecker delta and ·, · denotes the inner product in the Hilbert space of the variables ξ, f (ξ)g(ξ) = f (ξ)g(ξ)w(ξ)dξ . The weighting function w(ξ) is the probability density function of the uncertain parameters. Such polynomials already exist for some standard distributions which can be found in the Askey scheme [19] , for example, a Normal distribution leads to Hermite-chaos, while Legendre-chaos corresponds to a Uniform distribution. For other commonly used distributions or any arbitrary distribution, one can for example use Gram-Schmidt algorithm [17] to construct the orthogonal polynomials.

Substituting (2) in the general stochastic differential Eq. (1), we obtain

In order to ensure that the truncation error is orthogonal to the functional space spanned by the basis polynomials {ψ i }, a Galerkin projection of the above equation is performed onto each polynomial {ψ k },

After using the orthogonality property of the polynomials, we obtain a set of P + 1 deterministic coupled equations for all the random modes of the solu-

Following the definition, the mean and the variance of the solution are given by

As the coefficients v i (x, t) are known, the probability distribution of the solution can be obtained.

We first discuss the governing equations in the deterministic setting. The Navier-Stokes equations for an incompressible flow is given by

where u is the velocity, p is the pressure and ν is the kinematic viscosity. In Large Eddy Simulation, the reduction in the range of scales in a simulation is achieved by applying a spatial filter to the Navier-Stokes Eqs. [12] . This results into ∂u ∂t

where u is the filtered velocity, p is the filtered pressure and τ = uu − u u is the so-called subgrid-scale (SGS) stress tensor. The subgrid-scale stress tensor represents the effect of the small (unresolved) scales on the resolved scales, and has to be modeled in order to close the filtered Navier-Stokes equations. A popular class of SGS models is the eddy-viscosity models. In order to take account of the dissipation through the unresolved scales, the eddy-viscosity models, locally increases the viscosity by appending the molecular viscosity with the eddy viscosity. Mathematically, these models specify the anisotropic part of the subgrid-scale tensor as

where ν t is the eddy-viscosity and S = (∇u + (∇u) T )/2 is the resolved strain tensor. Substituting into the filtered momentum Eq. (6), we obtain

where the incompressibility constraint is used to simplify the equation. The pressure here is altered to include the trace term of Eq. (8). Smagorinsky model [13] is one of the oldest and most popular eddy-viscosity SGS model. The eddy viscosity of the Smagorinsky model is expressed as

where C s is the Smagorinsky coefficient, Δ is the LES filter width and |S| = √ 2S : S. It should be noted that the coefficient C s must be known prior to the simulation and is usually adapted to improve the results [12] . For example, C s = 0.2 is used for isotropic homogeneous turbulence, while a value of C s = 0.1 is used in case of channel flow. Similar values (C s 0.1 − 0.12) are realized from the shear flow studies based on experiments [10] .

Let us consider the Navier-Stokes Eqs. (6) with some uncertainty in the system. The sources of the uncertainty considered here are boundary conditions, material properties and model parameters. We employ the IPC method (see Sect. 2) by presuming the dimensionality (d) and probability density function of the uncertain random variables {ξ 1 , ξ 2 , ...ξ d } to be known, allowing us to construct the finite set of orthogonal polynomial basis {ψ i }.

In order to obtain a rather generic formulation, unless specified otherwise, we consider uncertainty in all sources listed above. Thus, the associated polynomial chaos expansion (PCE) for kinematic viscosity is given by

Note that the coefficients ν i are assumed to be known. The dependence of the flow variables, i.e. velocity and pressure, on the stochastic variables is expressed by the following PCEs

where u i and p i are the unknown polynomial chaos mode strengths of velocity and pressure fields, respectively. For deterministic boundary conditions, u i and p i are all zero for i = 1, 2, ..., P . In case of uncertain boundary conditions with a known (or modeled) probability density function, u i and p i , for i = 0, 1, ..., P , can be estimated. Introducing these expansions in Eqs. (6) and taking a Galerkin projection onto each polynomial {ψ k } while using the orthogonality of the polynomial chaos and finally dividing by ψ k ψ k , results, for k = 0, 1, ..., P , into the following set of deterministic equations:

where

Note that the original system of Eqs. (6) is transformed into a system of (P +1) divergence-free constraints on velocity modes and (P +1) coupled equations in velocity and pressure modes. A detailed discussion on the solution procedure adopted for this large system of equations is deferred to Sect. 4. Similarly, for the filtered momentum Eq. (7) with Smagorinsky model for turbulence, applying the above steps, results in

where

Note that |S| 2 = 2S : S, and applying the IPC steps to this identity -using polynomial chaos expansion and projecting on each basis polynomial, we obtain

where S i is the resolved strain tensor based on i th velocity mode. The above corresponds to system of (P + 1) non-linear equations in the unknown expansion coefficients of |S|. This system is solved using Picard iterations with |S| i = 2S i : S i as the initial guess.

OpenFOAM uses the finite volume method (FVM) for the disretization of partial differential Eqs. [4] . Among the various fluid dynamic solvers offered by Open-FOAM, we choose the solver called pimpleFoam [2], which allows the use of large time-steps to solve the incompressible Navier-Stokes Eqs. (6) . This solver is based on the PIMPLE algorithm for pressure-velocity coupling using Rhie and Chow type interpolation [7] . From the previous section, it can be realized that the system of governing Eqs. (14) for the evolution of the velocity and pressure modes u k , p k for k = 0, 1, ..., P , has a structure similar to the original deterministic Navier-Stokes Eqs. (6) . Due to the coupling via convection and diffusion terms, the size of this new system is P + 1 times its deterministic version. It can be observed that the divergence-free velocity constraints are decoupled and can be solved independently. Based on this observation, a fractional step projection scheme has been previously implemented [8] . In the first fractional step, the convection and diffusion terms are integrated followed by enforcing the divergence-free constraints in the second fractional step.

Our approach of the stochastic solver is based on the development of the existing deterministic solver (pimpleFoam) such that it can accommodate and solve (P + 1) coupled Navier-Stokes like systems in u k , p k for k = 0, 1, ..., P . We solve each of these systems sequentially, by using the initialized/updated velocity and pressure modes, and repeat until convergence. Figure 1(a) , provides a graphical representation of this approach. Depending on the type of flow, the value of P and a few other parameters; it usually takes about 3-6 explicit iterations (I e ) to realize convergence at every time-step. The default value of I e is set to P + 1. Following the conventions of OpenFOAM, we call this solver, gPCPimpleFoam. In contrast to the fractional step scheme, this approach admits better stability, stronger coupling and faster convergence; with an efficient data management and minimal changes in the exiting solver.

In Fig. 1(b) , we highlight the most important steps needed to develop gPCPimpleFoam, over the existing solver, pimpleFoam. In contrast to the deterministic solver, two nested loops are introduced. The first loop is over the explicit iterations (I e ), which updates the mode strengths between the two consecutive time-steps. For every explicit iteration, the second nested loop solves the Navier-Stokes like system (u k , p k ) for each mode strength k, while employing the existing, however modified, structure of the PIMPLE scheme. The modifications are inevitable due to the summations in the convection and the diffusion terms of the stochastic equations. It is realized that the mean mode (u 0 , p 0 ) changes slowly as compared to the other modes. Thus, in order to increase the stability, we start solving the last system (u P , p P ) first, and updating all the modes before solving the first system (0 th mode) representing the mean.

In addition to the exiting modules, we require to either modify or create some completely new routines for turbulence, pre-and post-processing etc. Restricting the verbosity, we attempt to provide an overview of the major implementation steps: (a) Creation of new variables (vectors) for the list of mode strengths of all the uncertain parameters, flow variables and derived variables for post-processing, (b) For the solver to read these inner products (obtained using chaospy library [3] ), only a small routine is added. Another similar lines of code are added to read in the values of d, n, I e , etc, (c) Very subtle changes are needed in the transport model in order to read in the transport properties for all the mode strengths. To accommodate for reading and initializing the turbulence models, a few minor changes are made in the LES model library. Significant changes are required specially for the Smagorinsky model which allows automatic reading of all the mode strengths of C s from the input file, and (d) To estimate the mean and the variance of flow variables and other derived quantities, a separate routine is added to the post-processing step of the solver. We use the IPC steps to calculate the derived quantities (like Reynolds Stresses) using the resolved expansion coefficients of the flow variables and the known coefficients of other parameters.

The order of accuracy and the convergence rate are governed by the large variety of space and time disretization schemes and iterative solvers offered by the OpenFOAM library.

We first consider a 2D steady laminar flow in long rectangular channel (with a height of 2δ) in the absence of any external forces. For a given average inlet velocity u avg , the fully developed flow has an analytical solution known as the Hagen-Poiseuille solution.

Therefore, the velocity field is independent of the viscosity and at y = 0, u = u max = 3 2 u avg . We assume the boundary condition to be deterministic and presume a known PCE for the uncertain kinematic viscosity, ν = P i=0 ν i ψ i (ξ). Also, we consider a Gaussian random variable to model the viscosity, for which the associated polynomial chaoses ψ i (ξ), are the Hermite polynomials. Using PCE of pressure gradient and random viscosity, for i = 0, 1, ..., P , we obtain

The use of polynomial chaos for incompressible laminar flow in a 2D channel has been previously investigated [8] , and as a validation case, we carry out a similar study with the IPC solver developed using OpenFOAM.

A uniform velocity is used at the inlet with no-slip boundary conditions at the top and bottom walls, and the gradient of velocity is set to zero at the outlet. Note that the use of deterministic boundary condition implies, for i = 1, ..., P , the unknown mode strength and/or their derivatives are by default set to zero at the boundaries. The Reynolds number, Re = 2δu avg /ν 0 , is set to 100. Fifthorder 1D Hermite polynomials are employed for all the PCEs, i.e. P = 5. The coefficient of variation (CoV = σ/μ) for the uncertain viscosity is set to ∼20%, with ν 1 /ν 0 = 2 × 10 −1 and ν 2 /ν 0 = 8 × 10 −5 . The remaining mode strengths of viscosity are assigned a value of zero. The simulation is performed in a domain with L/δ = 50 and a 250 × 100 mesh with a near-wall grading. Second-order disretization schemes are used both in space and time, and time-step size Δt = 10 −2 δ/u avg is specified. Figure 2 shows the profile of the mean and the standard deviation of velocity. The mean depicts the gradual transition in the flow along the channel length from a uniform inlet profile to a parabolic profile with maximum on the centerline. The uncertainty in velocity at inlet is zero, which is indeed the consequence of the deterministic boundary condition. In the developing region, a higher standard deviation is realized in the channel center as well as in the boundary layer with two lobes close to the walls. A significant variation in the modes (and thus the standard deviation) is realized up to 10-12 channel half-widths and further downstream, these modes become less significant. Figure 3 provides the axial velocity profile with confidence region (±2σ) at different locations in the downstream direction. The uncertainties tend to zero in the fully developed region, which is in accordance with the theory where the velocity is independent of the viscosity (see Eq. (16)). Figure 4 shows the estimated ratios of the modes of pressure gradient with respect to the mean pressure gradient along the channel centerline. Clearly, after the recirculating regions near the channel inlet, these ratios gradually reach their constant values further downstream. For x/δ > 20, the results are identical to analytical predictions (see Eq. (17)), characterizing the uncertainty in pressure due to the uncertainty in viscosity. Velocity mode strengths are shown in Fig. 5 . As evident, the results from the intrusive variant are in accordance with the non-intrusive counter-part, and due to the fast spectral convergence of the polynomial chaos representation, the magnitudes of the modes decrease as P increases.

A turbulent channel flow is a theoretical representation of a flow driven by a constant pressure gradient between two parallel planes extending infinitely. A 3D schematic is shown in Fig. 6 . Since the computational domain has to be finite, in addition to channel width h, we fix the stream-and span-wise truncation lengths, l x and l z , respectively. The values of h, l x and l z are adopted from [1] . These values ensures that the computational domain is large enough to accommodate the turbulent structures in the flow.

In order to maintain an equivalent flow, instead of the pressure gradient, the bulk velocity can also be prescribed, U b = 1 h h 0 u dy. The bulk Reynolds number is then defined as Re b = hU b /ν. In context of turbulent channel flows, another characteristic velocity called the friction velocity is usually introduced in terms of the wall shear stress τ w and the fluid density ρ, as u τ = τ w /ρ. The friction Reynolds number is then defined as Re τ = δu τ /ν, where δ = h/2 is the channel half-width. Then the pressure gradient and the wall shear stress relates as − dp dx = τw δ , wherep/ρ = p. Thus we have a choice between prescribing the bulk Reynolds number via bulk velocity and the friction Reynolds number via pressure gradient. Since we have to study the effect of the uncertain model parameter on the flow profile, we decide to fix the pressure-gradient and compute U b through (stochastic) simulation. The stochastic LES Smagorinsky model, as discussed in Sect. 3, is employed to solve for the turbulence. Since the model parameter C S may take a range of values, we assume it to be uncertain with a known PCE. We consider a Uniform random variable to model the parameter, for which the associated polynomial chaoses ψ i (ξ) are the Legendre polynomials. Table 1 summarizes the physical parameters used in the deterministic and stochastic simulations. 1D third-order Legendre polynomials are used for all the PCEs. The value of C S1 is set equal to the standard deviation, while the remaining mode strengths of C S are assigned a value of zero. Periodic boundary condition are applied in the stream-and span-wise directions, while no-slip boundary condition is used at the walls. The simulation results will be compared to the Direct Numerical Simulation (DNS) data from [9] at Re τ = 395. Based on the study of the effect of computational grid size from [1] , we use a reasonably fine mesh with the details in Table 1 . Note that Δx + = Δxu τ /ν, Δz + = Δzu τ /ν and y + = yu τ /ν, are calculated using value of u τ , corresponding to the value of Re τ in the DNS database. In order to capture the sharp gradients in the near-wall region, we specify a grading along y direction. We use the van Driest damping function to correct the behavior of the Smagorinsky model in the near-wall region [1] . Figure 7 presents the normalized time-averaged streamwise component of the velocity along with the mean and the confidence interval. The profile is compared with the DNS data, the deterministic solution (DET) at the mean value of C S and also with the results from non-intrusive polynomial chaos. The stochastic mean from IPC is found to be close to DET and mean of NIPC, and deviates slightly from DNS in the same manner as the DET solution. In contrast to IPC, the NIPC approach under predicts the variance. The uncertainty in the LES model parameter is reflected in the solution in the regions close to the wall and the channel center. In Fig. 8 the normalized square-root of the second order velocity moments and Reynolds shear stress are plotted together with their confidence regions. As evident, the stochastic mean of stresses are close to that of NIPC and DET solution. The deviation from DNS can mainly be attributed to the use of a relatively coarse mesh and the choice of LES model. Both the IPC and NIPC methods predicts high variance near the wall, with almost zero uncertainty in the channel center. This is expected as the Smagorinsky model parameter, when changed, usually affects significantly near the wall as compared to the channel center. While, the confidence region of intrusive method mostly overlaps with that of the non-intrusive counterpart, in some regions, NIPC still underestimates the uncertainty.

The IPC method, as it involves solving a lesser number of equations than NIPC, can be of great use when the deterministic simulation is already computationally expensive. To develop an IPC solver it is important to efficiently decouple the system of equations and ensure the overhead due to coupling is not significant.

In this work, intrusive polynomial chaos for CFD simulations using a popular finite-volume library OpenFOAM was presented. The aim was to use the existing deterministic solver for the incompressible Navier-Stokes equations, to develop a new stochastic solver for the quantification of the uncertainties involved and study its non-linear propagation.

To this end, we tested this solver for various standard CFD problems involving laminar and turbulent flows. The plane Poiseuille flow with uncertain kinematic viscosity was discussed first. Here, we realized a significant effect of the uncertainty in the re-circulation region of the flow. The results were also compared with the non-intrusive counterpart with same polynomial order. The results from IPC were found to be very close to NIPC, verifying its implementation using OpenFOAM. Next, we examined the turbulent channel flow with uncertain LES model parameter. The results were found to be in accordance with the non-intrusive gPC method. Deviations in the variance predicted by the two variants of gPC approaches can be attributed to the use of the two very different numerical methods to estimate the expansion coefficients in both variants.

Through this work, for UQ in CFD, we bring to light an alternate to the MC and the NIPC approaches. The promising results obtained from IPC method encourages to further pursue research in this direction. Pseudo-spectral methods along with early truncated expansions will be some next steps to reduce the computational cost. As a future work, Parallel-in-Time methods will be tested to overcome the saturation via space-only parallelization.

The potential of large eddy simulation for the modeling of wall bounded flows Eugene de Villiers

Chaospy: an open source tool for designing methods of uncertainty quantification

Computational Methods for Fluid Dynamics

Stochastic Finite Elements: A Spectral Approach

Non-intrusive polynomial chaos methods for stochastic CFD-theory and applications

Solution of the implicitly discretised fluid flow equations by operatorsplitting

A stochastic projection method for fluid flow. I. Basic formulation

Direct numerical simulation of turbulent channel flow up to Re 590

Subgrid-scale stresses and their modelling in a turbulent plane wake

Spectral Methods for Uncertainty Quantification with Applications to CFD

Large Eddy Simulation for Incompressible Flows

General circulation experiments with the primitive equations

Uncertainty Quantification: Theory, Implementation, and Applications

Uncertainty analysis for fluid mechanics with applications

The homogeneous chaos

Modeling arbitrary uncertainties using Gram-Schmidt polynomial chaos

Efficient collocational approach for parametric uncertainty analysis. Commun

The Wiener-Askey polynomial chaos for stochastic differential equations

Modeling uncertainty in flow simulations via generalized polynomial chaos

Acknowledgments. This work is supported by the Netherlands Organisation for Scientific Research (NWO) under the "Parallel-in-time methods for propagation of uncertainties in wind-farm simulations" project.