Application of the Double Potential Method to Simulate Incompressible Viscous Flows

Abstract

In this paper we discuss an application of the double potential method for modelling flow of incompressible fluid. This method allows us to avoid the known difficulties in calculating pressure and overcome the instability of numerical solution. Also, the double potential method enables us to simplify the problem of boundary conditions setting. It arises when computing the incompressible fluid flow by the Navier-Stokes equations in the vector potential - velocity rotor variables. In the approach given, the final system of equations is approximated through applying the finite volume method. In this case, an exponential transformation of the flow terms is applied. A parallel program was developed by means of using MPI and OpenMP technologies for the purpose of the numerical method computer implementation. We used two tasks to test. One of them deals with the classical calculation of the fluid flow establishment in a long round pipe. The other one is connected with the flow calculation in the pipe that in the output region contains a separation into two symmetrical parts. To perform numerical simulation, we take into consideration the steady flow with Reynolds numbers of 50 and 100. The numerical results obtained are consistent with computational results received through using the ANSYS CFD package.

1 Introduction

Modeling internal flow of viscous incompressible fluid is one of the most important and complex problems in continuum dynamics. This challenge implies the great potential to deal with applied tasks. For example, the problem of cleaning aquatic environment from impurities of the fine iron ions can be solved by applying the method of affecting an electromagnetic field on a collector considered in two-dimensional formulation in [1, 2]. To obtain the distribution of the impurities by the grid method we need to have a velocity field in the computational domain.

When simulating internal flow problems in 2D geometry, the Navier-Stokes model is widely used in the stream-vortex function formulation. It allows us to avoid the setting limitations in natural variables (speed-pressure). These limitations imply the discretization complexity of the equation to determine the pressure field, the high instability of the solution to this equation when using numerical simulation cellular methods and no performance guarantee of the mass conservation law [3]. This technique enables us to reduce the number of dependent variables. To receive the natural formulation, we use two coordinates of velocity and pressure. To get the current function – vortex formulation we have only two scalar variables.

Many sources suggest the lack of a generalization of such method for the three-dimensional problem [5,6,7,8,9,10,11,12,13,14,15,16,17,18], but the work [4] shows that it is possible. It is achieved through determining the vector potential and the vector vortex like the current function and the vortex for the two-dimensional variant. It should be noted that in this case, the formulation of the boundary conditions for the vector potential is nontrivial. The boundary conditions are extremely important for computational hydrodynamics and, as it is shown in [4], they have a determining influence on the resulting view of flow. To simplify the boundary conditions formulation for incompressible fluid modeling internal flow, the double potential method [19,20,21], was developed and successfully applied in [21]. The approach is discussed in this work.

2 Mathematical Model

To simulate the flow of viscous incompressible fluid, we use the Navier-Stokes system of equations and the incompressibility condition [3]. The dimensionless formulation of the equations has the following form:

$$ \frac{{\partial {\mathbf{u}}}}{\partial t} + \left( {{\mathbf{u}},\nabla } \right){\mathbf{u}} = \frac{1}{\text{Re}}\Delta {\mathbf{u}} - \nabla p $$

(1)

$$ \nabla {\mathbf{u}} = 0 $$

(2)

here u is speed vector, $\frac{\partial }{\partial t}$ is a derivative of a function with respect to time, $\text{Re} = \frac{{u_{0} \rho_{0} D_{0} }}{\eta }$ is the Reynolds number, Δ is the Laplace operator, ∇ is the Hamilton operator, p is pressure, u₀ is characteristic speed of flow, ρ₀ is medium density, D₀ is hydraulic diameter, η is coefficient of dynamic viscosity.

The system of Eqs. (1), (2) is completed with corresponding boundary and initial conditions. They will be presented further when obtaining the system of equations of the double potential method. Let’s consider the main calculations of this method.

From the Eq. (2) and the rotor property $\nabla \left( {rot\left[ {\mathbf{A}} \right]} \right) = 0$, we can write the formula for the vector A:

$$ {\mathbf{u}} = rot\left( {\mathbf{A}} \right) $$

(3)

When condition (3) is fulfilled, vector A is called the vector potential of flow. A characteristic feature of the viscous fluid motion is vorticity, which is defined as follows:

$$ {\varvec{\upomega}} = rot\left( {\mathbf{u}} \right) $$

(4)

Basing on the condition (3) we can write formula $rot\left( {rot\left[ {\mathbf{A}} \right]} \right) = {\varvec{\upomega}}$. Then, we use the rotor property: $rot\left( {rot\left[ {\mathbf{A}} \right]} \right) = \nabla \left( {\nabla {\mathbf{A}}} \right) - \Delta {\mathbf{A}}$, require that $\nabla {\mathbf{A}} = 0$ [4] and get equation to define the vector potential:

$$ \Delta {\mathbf{A}} = - {\varvec{\upomega}} $$

(5)

The equations for calculating the vortex are obtained by the action of the $rot$ operation on Eq. (1), then we take the final formula:

$$ \frac{{\partial {\varvec{\upomega}}}}{\partial t} + \left( {{\mathbf{u}},\nabla } \right){\varvec{\upomega}} - \left( {{\varvec{\upomega}},\nabla } \right){\mathbf{u}} = \frac{1}{\text{Re}}\Delta {\varvec{\upomega}} $$

(6)

So, Eqs. (3)–(6) with the boundary conditions allow us to calculate the velocity field in the computational geometry. This method is briefly described in [4] and it has applications for solving some problems of computational hydrodynamics. When using this method, there is a serious difficulty in setting the boundary condition on the vector potential. It can be seen from condition (3). In order to overcome this difficulty, a consequence of the theory of potentials was proposed in [19,20,21]. Write the speed in the form:

$$ {\mathbf{u}} = rot({\mathbf{A}}) - \nabla \varphi $$

(7)

here φ is scalar potential. To satisfy condition (2), it is necessary to require:

$$ \Delta \varphi = 0 $$

(8)

This condition is implemented as an equation to calculate the scalar potential.

The view of the velocity in the form of (7), in accordance with the theory of potentials, allows us to represent the boundary conditions on a simply connected computational domain for the vector potential in the form [19]:

$$ A_{\tau } = \frac{{\partial A_{n} }}{{\partial {\mathbf{n}}}} = 0,\quad {\text{for}}\;\delta \varOmega $$

(9)

here $A_{\tau }$ is the tangent component of the vector potential, $A_{n}$ is the normal component of the vector potential, n is normal to the border, δΩ is boundary of the computational area. The boundary conditions for the scalar potential are determined as follows:

$$ \frac{\partial \varphi }{{\partial {\mathbf{n}}}} = \left( {\left. {\mathbf{u}} \right|_{\delta \varOmega } ,{\mathbf{n}}} \right)\quad {\text{for}}\;\delta \varOmega $$

(10)

here $\left. {\mathbf{u}} \right|_{\delta \varOmega }$ is the flow velocity of medium across the boundary of the area. It is the boundary conditions for speed. The boundary and initial conditions for the vortex are determined from expression (4) as follows:

$$ {\varvec{\upomega}} = rot\left( {\left. {\mathbf{u}} \right|_{\delta \varOmega } } \right)\quad {\text{for}}\;\delta \varOmega , $$

(11)

$$ {\varvec{\upomega}} = rot\left( {\left. {\mathbf{u}} \right|_{t = 0} } \right)\quad {\text{for}}\;t = 0, $$

(12)

At last, the final formulation of the problem of flow modeling by the double potential method is represented by Eqs. (5)–(8) through using the boundary conditions (9)–(11) and the initial conditions (12).

3 Numerical Method for 3D Geometry

For 3D case, the grid consisting of triangular prisms is constructed in the computational domain (it is shown in Fig. 1). The finite volume method on the centers of the prismatic cells is used to approximate the equations by the double potential method (5)–(8) [22].

Fig. 1.

Computational domain (at the left) and grid consisting of triangle prisms (at the right).

To approximate the fluxes across the boundaries of the computational cell, we introduce the following notations (see Fig. 2): $V^{i}$ is the volume of the current i-th grid element, $P^{i}$ is the point of the center of the current i-th element, $P^{ij}$ is the point of the center of the j neighbor to the i-th element, $l^{ij}$ is the distance between the centers of neighboring cells, ${\mathbf{n}}^{ij} = \left\{ {n_{x}^{ij} ,n_{y}^{ij} ,n_{z}^{ij} } \right\}^{T}$ is the unit direction vector from the center of the current element to the center of the neighboring one, ${\mathbf{S}}^{ij} = \left\{ {S_{x}^{ij} ,S_{y}^{ij} ,S_{z}^{ij} } \right\}^{T}$ is the directional area of the j-th edge of the i-th element. The number of adjacent elements of a prismatic cell is five.

Fig. 2.

Characteristic grid values for prismatic elements.

To approximate the vector Eq. (6), we apply the exponential transformation similar to the two-dimensional case in [23] and write the resulting difference expressions:

$$ \bar{W}_{kp}^{ij} = \frac{{\bar{G}_{k}^{ij} \omega_{p}^{ij} - \omega_{p}^{i} }}{{0.5\left( {\bar{G}_{k}^{ij} - 1} \right)l^{ij} }}n_{k}^{ij} ,\quad \bar{G}_{k}^{ij} = \exp \left( { - \text{Re} l^{ij} n_{k}^{ij} \frac{{u_{k}^{ij} + u_{k}^{i} }}{2}} \right) $$

(13)

$$ \begin{array}{*{20}l} {\hat{\omega }_{p}^{i} = \omega_{p}^{i} + \frac{\tau }{{V^{i} }}\left[ {\sum\limits_{j = 1}^{5} {\left( {{\bar{\mathbf{W}}}_{p}^{ij} ,{\mathbf{S}}^{ij} } \right)} + \sum\limits_{j = 1}^{5} {\left( {\overline{{\varvec{\upomega}}}^{ij} ,{\mathbf{S}}^{ij} } \right)\bar{u}_{p}^{ij} } } \right],} \hfill \\ {\overline{{\varvec{\upomega}}}^{ij} = 0.5\left( {{\varvec{\upomega}}^{ij} + {\varvec{\upomega}}^{i} } \right),} \hfill \\ {\bar{u}_{p}^{ij} = 0.5\left( {u_{p}^{ij} + u_{p}^{i} } \right)} \hfill \\ \end{array} $$

(14)

here $k,p = x,y,z,\bar{W}_{kp}^{ij}$ is approximation of the special operator that acts on the vortex on the j-th face in the i-th cell and it is a (3 × 3) matrix, $\hat{\omega }_{p}^{i}$ is a difference analog of the p-th component of the vortex in the center of the current cell at the next time layer, $\omega_{p}^{ij}$ is the difference analog of the p-th component of the vortex in the center of the current cell at the current time layer, $\omega_{p}^{i}$ is the difference analog of the p-th component of the vortex in the center of the current cell on the current time layer, ${\bar{\mathbf{\omega }}}^{ij}$ is linear interpolation of the vortex in the middle of the face, $\bar{u}_{p}^{ij}$ is linear interpolation of the p-th coordinate of velocity in the middle of the face. The boundary conditions for the vortex are realized by discretization of the speed rotor with the first order of accuracy.

Further, we write the difference analogues for the equation of the scalar potential (8) and the vector potential (5):

$$ \hat{\varphi }^{i} = \sum\limits_{j = 1}^{5} {\frac{{\varphi^{ij} }}{{l^{ij} }}\left( {{\mathbf{n}}^{ij} ,{\mathbf{S}}^{ij} } \right)} \left[ {\sum\limits_{j = 1}^{5} {\frac{{\left( {{\mathbf{n}}^{ij} ,{\mathbf{S}}^{ij} } \right)}}{{l^{ij} }}} } \right]^{ - 1} $$

(15)

$$ \hat{A}_{k}^{i} = \left[ {\sum\limits_{j = 1}^{5} {\frac{{A_{k}^{ij} }}{{l^{ij} }}\left( {{\mathbf{n}}^{ij} ,{\mathbf{S}}^{ij} } \right)} + \hat{\omega }_{k}^{i} V^{i} } \right]\left[ {\sum\limits_{j = 1}^{5} {\frac{{\left( {{\mathbf{n}}^{ij} ,{\mathbf{S}}^{ij} } \right)}}{{l^{ij} }}} } \right]^{ - 1} $$

(16)

here k = x, y, z, $\hat{A}_{k}^{i}$ is difference analog of the k-th component of the vector potential in the center of the current cell on a new time layer, $A_{k}^{ij}$ is difference analog of the k-th component of the vector potential in the center of a neighboring cell on the previous time layer.

So, the difference analogue of the three-dimensional velocity vector is written as follows:

$$ \begin{aligned} & \hat{u}_{x}^{i} = \frac{0.5}{{V^{i} }}\left[ {\sum\limits_{j = 1}^{5} {\left( {\hat{A}_{z}^{i} + \hat{A}_{z}^{ij} } \right)}\; S_{y}^{ij} - \sum\limits_{j = 1}^{5} {\left( {\hat{A}_{y}^{i} + \hat{A}_{y}^{ij} } \right)} \; S_{z}^{ij} - \sum\limits_{j = 1}^{5} {\left( {\hat{\varphi }^{i} + \hat{\varphi }^{ij} } \right)} \; S_{x}^{ij} } \right], \\ & \hat{u}_{y}^{i} = \frac{0.5}{{V^{i} }}\left[ {\sum\limits_{j = 1}^{5} {\left( {\hat{A}_{x}^{i} + \hat{A}_{x}^{ij} } \right)}\; S_{z}^{ij} - \sum\limits_{j = 1}^{5} {\left( {\hat{A}_{z}^{i} + \hat{A}_{z}^{ij} } \right)} \; S_{x}^{ij} - \sum\limits_{j = 1}^{5} {\left( {\hat{\varphi }^{i} + \hat{\varphi }^{ij} } \right)} \; S_{y}^{ij} } \right], \\ & \hat{u}_{z}^{i} = \frac{0.5}{{V^{i} }}\left[ {\sum\limits_{j = 1}^{5} {\left( {\hat{A}_{y}^{i} + \hat{A}_{y}^{ij} } \right)} \; S_{x}^{ij} - \sum\limits_{j = 1}^{5} {\left( {\hat{A}_{x}^{i} + \hat{A}_{x}^{ij} } \right)} \; S_{y}^{ij} - \sum\limits_{j = 1}^{5} {\left( {\hat{\varphi }^{i} + \hat{\varphi }^{ij} } \right)} \; S_{z}^{ij} } \right] \\ \end{aligned} $$

(17)

here ${\hat{\mathbf{u}}}^{i} = \left\{ {\hat{u}_{x}^{i} ,\hat{u}_{y}^{i} ,\hat{u}_{z}^{i} } \right\}$ is difference analogue of the velocity vector in the center of the current cell.

The final algorithm to calculate the velocity field on 3D prismatic grid is the sequential calculations of difference expressions (13)–(17).

4 Parallel Realization

To parallelize the numerical method, we use two-level domain decomposition technique. At the first level, we divide the computational area into domains in accordance with number of the supercomputer system nodes. At the second level, we apply splitting of domains into subdomains in accordance with threads implemented in the cores of the central processors of each node.

The resulting program code was developed in the C++ programming language using MPI [24] and OpenMP [25] parallel technologies.

MVS-10P Supercomputer of JSCC RAS (see www.jscc.ru) was chosen for testing parallel software. The system has parameters:

• 207 nodes, each includes:

• 2 x CPU Intel Xeon 2,7 GHz microprocessors with 6 cores

• 2 x VPU Intel Xeon Phi 7110X microprocessors with 61 cores

• Peak performance is 523.8 TFlops

• Transfer is FDR InfiniBand

The results of parallel computations on the prismatic grid with 2500000 cells for simple geometry (round tube) are shown in Fig. 3. They show that constructed parallel software is quite effective (the cylindrical computational domain is shown in Fig. 1).

Fig. 3.

Comparison of the acceleration and ideal on the prismatic grid with 2500000 cells (at the left) and the efficiency and ideal on the prismatic grid with 2500000 cells (at the right).

5 Test and Results

Two tasks were used to test our approach. One of them deals with the classical calculation of the fluid flow establishment in a long round pipe (see Fig. 1). The other one is connected with the flow calculation in the pipe that in the output region contains a separation into two symmetrical parts (see Fig. 4).

Fig. 4.

Longitudinal section of the pipe for z = 0 (at the left) and the pipe cross section for $3 \le x \le 6$ (at the right).

Let’s consider the first of our tasks. In this case, the well-known stationary Poiseuille flow is realized in the three-dimensional cylindrical geometry. The boundary and initial conditions are written as follows:

$$ \begin{aligned} \begin{array}{c} \left. {\mathbf{u}} \right|_{x = 0} = \left\{ {1 - y^{2} - z^{2} ,0,0} \right\}^{T} ,\quad \left. {\mathbf{u}} \right|_{{|y^{2} + z^{2} | = 1}} \equiv 0,\quad \left. {\frac{{\partial {\mathbf{u}}}}{\partial x}} \right|_{x = 6} \equiv 0 \\ \left. {\mathbf{u}} \right|_{t = 0} \equiv 0. \\ \end{array} \end{aligned} $$

The test calculation was carried out on the prismatic computational grid consisting of 133632 elements - 87 steps along the X axis, 1536 triangles at the base in the YZ plane. The Reynold number used in the calculation is Re = 100. The Figs. 5, 6, 7 and 8 show the results of calculations. As can be seen from Figs. 5 and 6, we have got steady flow of Poiseuille throughout the volume of the cylinder.

Fig. 5.

The distribution of the scalar potential (at the left) and velocity modulus (at the right) in sections: XY and XZ.

Fig. 6.

The distribution of z-component of the vector potential in XY section (at the left) and the distribution of the velocity module in the array of YZ sections (at the right).

Fig. 7.

The distributions of the velocity modulus in the longitudinal section z = 0 were calculated by SIMPLEC method of the ANSYS CFD software tools for Re = 50 (at the top) and with the help of double potential method for Re = 50 (in the middle section) and Re = 100 (below).

Fig. 8.

The distribution of the velocity modulus in the longitudinal section z = 0. Distributions of the velocity modulus in cross sections YZ for Re = 50, calculated in ANSYS CFD using the SIMPLEC method (at the top) and calculated by the double potential method (at the bottom).

Let’s consider the solution to the second problem. Here, the numerical simulation of the stationary flow was performed for the Reynolds numbers 50 and 100. The obtained numerical results are shown in Figs. 7 and 8. For comparison, we present at the same figures the calculated data obtained through using the ANSYS CFD package [26] and based on the traditional SIMPLEC method [4]. Analysis of the results shows that qualitatively and quantitatively, the calculations performed by the traditional SIMPLEC method and by the double potential method are very close. We could obtain more accurate estimates if setting parameters of the SIMPLEC method were known.

6 Conclusion

This paper examines the double potential method to calculate the viscous incompressible fluid internal flows. The main advantage of the method is based on the formulation preservation of the Navier-Stokes equations in the vector potential - vector vortex variables. In addition to that, it enables us to avoid setting complex boundary conditions on the vector potential.

To implement the double potential method, the original difference scheme was developed. The scheme is based on the finite volume method on cell centers through using the flux vector exponential transformation. This scheme realizes both the double potential method in the three-dimensional case and its consequence in the two-dimensional formulation.

The numerical algorithm obtained was parallelized and the software was implemented in C++ be means of MPI and OpenMP technologies. The two test problems were examined by the generated program code. One of them referring to the classical problem deals with the establishing the Poiseuille flow in the long circular pipe. The other task deals with calculating the flow in a square tube divided into two symmetrical parts in the outlet region.

The calculation results demonstrate the great potential to apply the double potential method to model viscous incompressible fluid internal flows in technical systems of complex three-dimensional geometry. The comparison of the results given with the ANSYS CFD package calculated data showed that the double potential method presented allows us to calculate the internal flow of fluid with great accuracy. The parallel implementation of the method in complex geometry enables us to solve tasks of this class much faster than be means of standard methods.

The work was supported partially by Russian Foundation for Basic Research, projects No. 18-07-01292-a, 18-07-00841-a. This work was partially supported by the Academic Excellence Project of the NRNU MEPhI under contract with the Ministry of Education and Science of the Russian Federation No. 02.A03.21.0005.