Abstract
A two-dimensional stochastic solver for the incompressible Navier-Stokes equations is developed. The vorticity-stream function formulation is considered. The polynomial chaos expansion was integrated with an unstructured node-centered finite-volume solver. A second-order upwind scheme is used in the convection term for numerical stability and higher-order discretization. The resulting sparse linear system is solved efficiently by a direct parallel solver. The mean and variance simulations of the cavity flow are done for random variation of the viscosity and the lid velocity. The solver was tested and compared with the Monte-Carlo simulations and with previous research works. The developed solver is proved to be efficient in simulating the stochastic two-dimensional incompressible flows.
1. Introduction
In engineering fields, most models are represented as partial differential equations (PDEs), assuming all input data is perfectly known. Unfortunately, geometry and material characteristics for instance will present uncertainties. Under those conditions, the output data become also uncertain. To deal with propagation of the input data uncertainties to the output data, probabilistic models are more appropriate than deterministic ones. Several methods of solution are developed to assess the response due to the uncertainties. This response depends on two main factors: the first factor is the geometric domain discretization and the second is the discretization involved random process [1, 2]. The solution methods may be classified according to the first factor to meshless methods [3], stochastic finite difference methods [4], and stochastic finite element methods [5]. On the other hand, according to the second factor the methods of solution may be classified into Monte-Carlo (MC) simulations [6], perturbations [7], and Neumann series expansion method [2]. Perturbation and Neumann series expansion methods are limited to small perturbations and some simple cases. MC methods are simple and can be easily used in many fields. But, usually, they are only used as the last resort since they are inefficient and normally need lots of samples which consume great computational resources. MC simulations should be of order 104 or more for reliable comparisons [8].
In recent years, some analysis methods based on spectral expansion were developed. Polynomial chaos expansion (PCE) method is one of them. The basic theory of this method was proposed by Wiener in 1938, where it is called homogenous chaos. The pioneering development work of this method was done by Ghanem and Spanos in the solid mechanics field [9, 10]. They integrated the PCE with the finite-element method which is known as the spectral stochastic finite element methods (SSFEM). Recently, the SSFEM is one of the most widely used methods [7]. In recent years, the applications of SSFEM in computational fluid dynamics (CFD) are expanded. In [11] a stochastic finite-volume upwind technique is used to solve viscous Burgers’ equation with stochastic viscosity over a wide range of the mean viscosity. The stochastic developed solver has higher performance compared with the MC simulations. In [12], the PCE was implemented for the 2D incompressible stochastic Navier-Stokes equations to simulate the stochastic behavior of a lid-driven cavity flow under the influence of uncertainties. The driven velocity and fluid viscosity were supposed, respectively, to be the uncertain variables which have Gaussian probability distributions. It was found there that the influence of the driven velocity uncertainty is larger than that of viscosity. In [13, 14], the cavity flow with stochastic variation of the viscosity is also simulated. More applications in CFD are given in [15] and the references herein.
In 2D, the incompressible laminar flow can be modeled using the Navier-Stokes equations in different formulations. The most famous formulations are the primitive variables (velocity and pressure) formulation and the vorticity-stream function formulation [16]. There are advantages in using the vorticity-stream function formulation of the incompressible Navier-Stokes equations to compute 2D flows: the continuity equation is automatically satisfied, only one (vorticity equation) transport equation has to be solved, the streamlines of the flow are given by level curves of the stream function, the vorticity is a conserved quantity [17], and there is no difficulty in matching a pressure to the velocity field [18]. The basic method and the treatment of the vorticity boundary conditions are reviewed in [19]. There is no much research work in the literature that uses the vorticity-stream function formulation. This can be due to difficulties in extending this formulation to 3D and also due to some difficulties in the vorticity boundary treatment. Some research efforts using the vorticity-stream function formulation are given in [16, 17, 20]. In [21], the vorticity-stream function formulation is discretized on unstructured grids with the upwind finite-volume cell-centered technique, and it was applied successfully in simulating the 2D incompressible flow in the cavity.
Solution of the stochastic Navier-Stokes equations in vorticity-stream function formulation has only little attention in the literature. In [14], the Wiener-chaos expansion (WCE) is used to treat the vorticity-stream function formulation driven by Gaussian Brownian motion. They introduced a compression technique similar to the sparse tensor products to handle the constant flux of new random variable due to the Brownian motion. There is no consideration in the literature for vorticity-stream function formulation with random variation of the parameters and/or boundary values. The current work fills this gap. As we shall see in this work, the vorticity-stream function formulation when extended in the stochastic dimension has advantages in the storage and computational complexity over the stochastic extension of the primitive variables formulation. In the current work, a node-centered finite-volume solver is developed and then extended in the stochastic dimension by using PCE to simulate the 2D incompressible flows with random variations in the viscosity and/or the boundary conditions.
A rigorous theory of the stochastic Navier-Stokes equation has been a subject of several papers [22]. Most papers rely on martingale-type methods and a direct theory of strong solutions providing the stochastic analog of the well-known Lions and Prodi [23] solvability theorem for the deterministic Navier-Stokes equation. Proofs of the existence and uniqueness of a strong local solution of stochastic Navier-Stokes equations are given in [24]. Also, existence and uniqueness of a global strong solution in 2D are established in [25]. The existence of a unique invariant measure is established and the properties of this measure are described in [26].
This paper is organized as follows: the governing equations with the boundary conditions are outlined in Section 2. The deterministic solver is derived in Section 3. The PCE and the extension to the stochastic dimension are explained in Section 4. The test cases, results, and discussions are given in Section 5. Summary of the work and conclusions are outlined in Section 6.
2. Governing Equations
The stream-function-vorticity formulation of the incompressible Navier-Stokes equations in 2D Cartesian coordinates, fixed boundaries, and neglected source terms are: where is the stream function, the components of are the horizontal and vertical velocities, respectively, of the flow field, and is the scalar vorticity. The constant is the kinematic viscosity. On the boundary of rigid bodies, the stream function must satisfy two boundary conditions: one of Dirichlet type corresponding to a no-normal-flow condition and one of Neumann type corresponding to a no-slip condition. No boundary conditions are explicitly given for the vorticity [17]. These equations are easily derived by taking the curl of the Navier-Stokes equations in primitive variables formulation and eliminating any dependence in the third () component.
3. Derivation of the Deterministic Solver
A 2D finite-volume node-centered unstructured solver is considered. The control volume is taken around the nodes by connecting the cell centroids surrounding each node . Figure 1(a) shows a part (shaded) from the control volume. In case of the near-boundary nodes, the boundary faces are considered as cells of zero area. The area of the control volume around node will be denoted as . Another quad control volume, Figure 1(b), is used in obtaining the gradient vector of any scalar variable . The quad control volume is consisting of the endpoints, and , of the face and the centroids of the two neighbor cells, nodes 1 and 2, as shown in Figure 1(b). The area of the quad control volume will be denoted as .
(a)
(b)
The components of a gradient vector of a scalar quantity are needed at the face centroids . For any scalar variable , a second-order midpoint approximation of the horizontal component is obtained by integrating over the quad control volume, Figure 1(b), and after applying Green’s theorem. The same procedure can be done for the vertical component of the gradient vector. This will give where is the vector connecting the face endpoints, and . The vector is the vector connecting the centroids (nodes 1 and 2) of the two neighbor cells surrounding the face. The gradient vector will be needed also at the nodes. This can be done by integrating over any control volume around the node, for example, the control volume in Figure 1(a). At the cell centroids, the gradient vector can be approximated by averaging it from the cell nodes.
The integration over the control volume of the Laplacian operator of any scalar variable around node is obtained, after applying Green’s theorem and using (4), as where is the value at node , is the number of neighbor nodes around , is the value at each neighbor node , is the number of neighbor cells around , and is the value of the neighbor cell centroid . The coefficients , , and are
The nonboundary nodal values will be treated implicitly (from the current time level) and the cell-centered values will be treated explicitly (from the previous time level) by averaging them from the nodal values from the previous time level.
The convective term in the vorticity transport equation is integrated over the control volume . The discretization of this integration will be equal to the summation of the exit fluxes from all surrounding faces of the control volume ; see Figure 1(a). The variable is the value of at the centroid, , of the face and is the mass flow rate across the face and in the direction of the face outer normal vector. The exit flux across each surrounding face can be approximated using a second-order upwind scheme as where and . The gradient vectors and are computed at nodes and , and they will be obtained from the previous time level. The vectors and are connecting the face centroid with the nodes and , respectively.
Integration of the time derivative term over the control volume around node can be computed with the first-order forward difference after assuming it constant over the control volume. This will give
Assembly of the stream function (2) at each node yields where is the average vorticity over the control volume . Assembly of the vorticity transport equation at each node yields
The superscript denotes the current time level while denotes the previous time level. Equation (9) is solved for the stream function and then the velocities are obtained from the gradient components of the Stream function, (3). The outlet velocities will be updated before they are used with the other boundary velocities to get the vorticity at all boundaries. The boundary vorticities are then used with the vorticity gradients from the previous time level to solve (10) and get the vorticity at the nonboundary nodes. Time iterations are done until the required convergence to the steady state is satisfied.
4. Random Discretization Using PCE
Polynomial chaos expansion (PCE) provides orthogonal basis for the space of second-order random variables, where the stochastic quantities are uniquely identified in terms of their coordinates with respect to the basis. Ghanem and Spanos [9] evaluate the stochastic system response as a summation of nonlinear functional of a set of centered, normalized, mutually orthogonal, and Gaussian random variables multiplied by deterministic coefficients. PCE has advantages in evaluating both statistical moments of any order and the probability distribution function (p.d.f) of the stochastic system response.
There are two significant numbers to define in the PCE; so-called order () and dimension (). The order depicts the maximum power used in the basis and the dimension describes the number of the used basis in the expansion. Thus, PCE of order and dimension consists of all polynomials of order , involving all possible combinations of the random variables .
The system variable in terms of PCE is written, after truncation at order and random variables , in the form: where are deterministic kernels and is the number of orthogonal polynomials of the PCE of order and dimension . This number is computed as [9]:
The polynomial of order and dimension is defined as [9]:
Stochastic variation of any variable, parameter, and/or boundary conditions will cause all other variables to be stochastic and hence are expanded as in (11) using PCE.
Assuming that the current system ((1)–(3)) has one or more parameters and/or boundary conditions that have stochastic variations, this will lead to a stochastic system. Following the procedure given in [9, 11, 15] by expanding all variables and parameters using PCE, (11), substituting in (1), (2), and (3) and taking the ensemble average after multiplying by ; to get the following stochastic system: where the stochastic multiplication tensor is the ratio between the ensemble average of the three polynomials, , and and the variance of [15], that is,
As it is clear in (14), the stochastic stream function equation is decomposed into independent linear systems with the same coefficient matrix resulting from the deterministic solver. This coefficient matrix is fixed (only dependent on the geometry) and diagonally dominant symmetric positive definite that can be solved efficiently using a conjugate gradient solver or by the SSOR. This means that we shall have only one linear system with multiple right hand sides. The storage cost will be and the computational cost will be , where is the number of degrees of freedom (number of nonboundary nodes in our case) and and are the order and dimension, respectively.
The stochastic vorticity transport equation should be handled carefully as we have the velocities , and the vorticity gradients are from different stochastic levels (levels and ). Note that we have stochastic levels. In the current work, the developed deterministic solver, described in Section 3, is modified to handle variables from different stochastic levels. The deterministic solver will be called with two indices, and , that represent two stochastic levels. A local matrix is then returned and will be multiplied by the multiplication tensor element before it will be inserted in the global matrix. The sparsity of the global matrix is preserved by mapping only nonzeros of the local matrices to the global matrix. This is important for saving storage memory and in enhancing the performance of the algebraic solver.
As in [13], when using the PCE with a 2D transport equation, the storage cost grows polynomially as and the computational cost grows as . In the vorticity-stream function formulation, only the vorticity transport equation has these expensive storage and computational costs. When using the primitive variables formulation, we shall have two expensive transport equations. This is an advantage in using the vorticity-stream function formulation over using the primitive variables formulation when considering 2D incompressible flow with random variations.
Figure 2 shows the sparsity patterns of two global matrices using different order and dimension using a unit cavity grid with 500 triangular elements. As we can notice, the sparsity ratio (ratio of nonzero entries) decreases with the increase of the stochastic order and/or dimension .
(a)
(b)
The linear system resulting from the discretized stochastic vorticity transport equation will be a large sparse system. A suitable efficient solver should be chosen. The PARDISO solver [28, 29] is used in the current work. PARDISO is a direct sparse solver that runs in parallel by automatically utilizing the multicore architecture of the used workstation. Other conjugate gradient solvers, such as BiCGSTAB, are also tested but PARDISO was found to be robust and more efficient in our work. More comparisons between different algebraic sparse solvers for the vorticity transport equation are given in [16] for the deterministic solver. The same results are also held for the stochastic solver.
The MC simulations are used in comparing and validating the stochastic solver. MC simulations are done also in parallel using the multithreading to exploit the available multicore architecture of the recent workstations.
5. Results and Discussion
The stochastic solver, described in Section 4, was tested and applied to the cavity flow. A unit lid velocity is considered with zero horizontal and vertical velocities on other boundaries. The stream function will be zero for all boundaries. The boundary vorticity is computed from the curl of the velocity. The stochastic variation of the viscosity is chosen to be in the form: where is the mean viscosity. This means that, for first order and first dimension , the viscosity standard deviation is 10% of the mean value . Unstructured triangular grid of 500 elements is used in the simulations. More finer unstructured triangular grids of 1900 and 3000 elements are also used in the simulations and in testing the developed solver. Figure 3 shows the convergence history of the stochastic solver on different grid sizes (500 and 1900 elements) with different values of the mean viscosity . As it is expected, the number of iterations to reach a certain error is increased as the mean viscosity decreases. Also, more iterations are required to reach a certain error as the grid size increases.
(a)
(b)
The mean of the horizontal velocity component is tested and compared with the known benchmark solution in [27]. Good agreements are obtained when comparing the horizontal velocity component at the vertical midline for mean viscosity of 0.01 and 0.001 as shown in Figure 4.
(a)
(b)
The stochastic solver is tested and validated against the MC simulations. For a fair comparison between the two approaches, the deterministic solver used in the MC simulations is taken to be the same stochastic solver but with zero order and zero dimension, that is, . The standard deviation relative error between the stochastic solver and the MC simulations for the horizontal and vertical velocity components was measured over the whole solution domain as shown in Figure 5. At a mean viscosity of 0.01, the maximum relative error in the mean is 0.035% for the horizontal velocity component and is 0.014% for the vertical velocity component. The figure shows also that the relative errors are larger in the regions with small values of the velocity components.
(a)
(b)
The p.d.f of the velocity components, vorticity, and the stream function are simulated and compared between the stochastic solver (with and ) and the MC simulations as shown in Figures 6 and 7. The p.d.f comparisons of the vorticity and stream function are done at and . As it is shown in the figures, the distributions are near from the standard Gaussian distribution. Good agreements are obtained for both the mean and variance even with 100 MC simulations. Increasing the MC simulations enhances the agreement but it will be very expensive both in time and storage. In the current work, the MC simulations are done in parallel (multithreading) to exploit the current available multicore architecture.
(a)
(b)
(a)
(b)
Simulations of the vorticity and stream function mean values for a mean viscosity of 0.01 and 0.001 are shown in Figure 8. Figure 9 shows the mean values of the two velocity components and for mean viscosity of 0.01 and 0.001, respectively. The standard deviations of the vorticity and the stream function using the current stochastic solver and for mean viscosity of 0.01 and 0.001 are shown in Figure 10. Figure 11 shows the standard deviations of the two velocity components for mean viscosity of 0.01 and 0.001, respectively. Qualitatively, there is a good agreement with the previous work in [12]. The quantitative comparison for the mean horizontal velocity was shown in Figure 4. Comparisons of the standard deviations are given in Tables 1, 2, and 3 and Figure 12.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
It is worth to note that Figures 8–11 are identical in case of MC simulations. The quantitative comparisons between using PCE and MC computations are shown Figures 5–7.
Table 1 shows good agreement in the standard deviation of the horizontal velocity component specially with the results in [14]. Table 2 is for the standard deviation of the vertical velocity component. Table 3 compares the standard deviation of the horizontal velocity component in case of and . As we can notice only small changes occurred.
Also, random variation of the lid velocity is considered and compared with [12]. The standard deviation of the lid velocity was taken to be 10% of the mean lid velocity. Figure 13 shows good agreements with the solution in [12]. The mixed random variation of both the viscosity and the lid velocity is evaluated. Comparison between the standard deviation of the horizontal velocity component and 10% standard deviation in the viscosity, 10% standard deviation in the lid velocity, and mixed effect of 10% standard deviation of both viscosity and the lid velocity is shown in Figure 14. As it is clear, the mixed effect is less than the unique variation in the lid velocity but greater than the unique effect of the random viscosity variation. We may notice that for the mixed effect is the superposition of both unique effects, while, for , the mixed effect is the difference between the two unique effects. Also, the effect of different values of the standard deviation of the lid velocity is tested; 2.5%, 5%, and 10% of the mean lid velocity, as shown in Figure 15. The maximum value of the standard deviation of the horizontal velocity component was 0.0076, 0.0152, and 0.0304, respectively, and occurs at in all cases. Figure 16 shows a comparison of the standard deviation of the vertical velocity component for different values of the viscosity standard deviation; 2.5%, 5%, and 10% of a mean viscosity of 0.01. We can notice that the stochastic response at the midline only scales up and down but does not change its profile.
6. Summary and Conclusions
The vorticity-stream function formulation of the 2D incompressible Navier-Stokes equations with random variation in the viscosity and/or the boundary conditions is considered. A node-centered finite-volume stochastic solver is developed and tested on the cavity flow. The results are compared with the MC simulations and with previous research works with good agreements. The effect of different stochastic variations of flow viscosity and/or lid velocity is quantified. The developed solver was shown to be efficient in simulating the stochastic 2D incompressible flow.
Acknowledgments
The authors thank Prof. M. A. El-Tawil (passed away) for the continuous support and motivation for this work. Thanks are due also to Prof. E. F. Abdel-Gawad and Dr. O. H. Galal for the continuing support.