ABSTRACT

Among the several methods used to solve the Navier-Stokes equations Hierarchical Expansion Method has demonstrated satisfactory results. This work aimed to apply the expansion of the variables in hierarchical functions for the solution of the Navier-Stokes equations for incompressible fluids in two dimensions in laminar flow. This method is based on the finite element method. The expansion functions in this study were based on Legendre polynomials, adjusted in the rectangular elements in such a way that corner, side and area functions were defined. The order of the expansion functions associated with the sides and with the area of the elements is adjusted to the necessary or desired degree. This method is denominated by Hierarchical Expansion Method. In order to validate the proposed numeric method three well-known problems of the literature in two dimensions were analyzed; however, for this paper only one problem was presented. The results demonstrated that method was able to provide precise results. From the results obtained in this paper it is possible to conclude that the hierarchical expansion method can be effective for the solution of fluid dynamic problems that involve incompressible fluids.

Keywords:

Finite Element, Petrov-Galerkin Formulation, Navier-Stokes Equations, Hierarchical Expansion Functions, Incompressible Fluid

1. Introduction

Among the several methods used to solve the Navier-Stokes equations Hierarchical Expansion Method has demonstrated satisfactory results. This work aimed to develop and apply the method of hierarchical expansion proposed by Zienkiewicz and Morgan [1] for the solution of the two-dimensional Navier-Stokes equations, for incompressible fluids in laminar flow. This method is based on the finite element method using the Petrov-Galerkin formulation [2] . All variables describing the fluid flow are expanded in hierarchical functions based on Legendre polynomials. This method has some advantages over other numerical methods as described in the following sections.

A recent literature review has demonstrated the efficiency of the method proposed in this study. Williams [3] reports that a particular space-time, hybridizable discontinuous Galerkin method is entropy stable for the compressible Navier-Stokes equations. In order to facilitate the proof, “entropy variables” are used to rewrite the compressible Navier-Stokes equations in a symmetric form. The resulting form of the equations is discretized with a hybridizable discontinuous finite element approach in space and a classical discontinuous finite element approach in time.

Chen et al. [4] demonstrates that a unit operator (first level) and an orthogonal project operator (second level) are constructed to act as the stability terms for Meshless Local Petrov-Galerkin (MLPG) method, which is called a two-level variational multiscale MLPG (VMS-MLPG) method. The VMS-MLPG method is applied to eliminate oscillation, overshoots and undershoots of MLPG method at large Pe. The results showed that the present VMS-MLPG method could guarantee the stable and reasonable solutions of convection-diffusion problems with large Pecle.

2. Conservation Equations (2D)

The equations governing the fluid dynamics are continuity, momentum and energy equations. In this work, the method proposed by Chorin [5] is used for the treatment of the coupling between pressure and velocity. In this case, a mathematical artifice is used to solve the equations of the momentum for each direction, which is divided into two equations. The first equation relates the fluid velocity components u and w, in terms of the components of a pseudo-velocity, called u* and w*. The second equation calculates the pressure gradient as a function of the pseudo-velocity components, u* and w*.

The following equations were used with this method:

$\rho \left(\frac{u^{\ast }-u^{t-\Delta t}}{\Delta t}\right)+\rho u\frac{\partial u^{\ast }}{\partial x}+\rho w\frac{\partial u^{\ast }}{\partial z}-\mu \left(\frac{{\partial }^2{u}^{\ast }}{\partial x^2}+\frac{{\partial }^2{u}^{\ast }}{\partial z^2}\right)=0$ (1)

$\rho u=\rho u^{\ast }-\Delta t\frac{\partial p}{\partial x};$ (2)

$\rho \left(\frac{w^{\ast }-w^{t-\Delta t}}{\Delta t}\right)+\rho u\frac{\partial w^{\ast }}{\partial x}+\rho w\frac{\partial w^{\ast }}{\partial z}-\mu \left(\frac{{\partial }^2{w}^{\ast }}{\partial x^2}+\frac{{\partial }^2{w}^{\ast }}{\partial z^2}\right)=0$ ; (3)

$\rho w=\rho w^{\ast }-\Delta t\frac{\partial p}{\partial z};$ (4)

${\nabla }^{2}p=\frac{\rho }{\Delta t}\left(\frac{\partial u^{\ast }}{\partial x}+\frac{\partial w^{\ast }}{\partial z}\right)$ ; (5)

where r is the fluid specific mass, t is the time, Dt is the time step, u and w are the velocity components in the x and z directions respectively, u* and w* are the intermediate values of the velocity components in the x and z directions respectively, p is the fluid pressure, μ is the fluid dynamic viscosity, and the superscript $t-\Delta t$ refers to the variables calculated at the previous time step. To make the nomenclature easier no superscript is used to denote the variables at the present time.

Assuming constant fluid properties, the energy equation written in terms of the temperature for Cartesian coordinates in two dimensionis the following:

$\rho \frac{\partial T}{\partial t}+\rho u\frac{\partial T}{\partial x}+\rho w\frac{\partial T}{\partial z}=\frac{k}{c_p}{\nabla }^{2}T+\frac{\mu }{c_p}{\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)}^2+\frac{2\mu }{c_p}\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial w}{\partial z}\right)}^2\right]$ , (6)

where T is the temperature, k is the thermal conductivity and c_p is the specific heat at constant pressure.

3. Mathematical Development of the Conservation Equations for a Single Element

A rectangular structure mesh was used to divide the flow domain into nodes. Although this is the simplest form of grid, it is enough to show the capabilities of the proposed method.

For each node, Equations (1)-(5), and Equation (6) are multiplied by a weighting function and integrated. The Green’s Theorem was used to simplify the diffusive terms in Equations (1)-(5), and Equation (6).

The weighting function, P_m, is based on the Petrov-Galerkin formulation and it is given, according to Brooks and Hughes [2] , by:

$P_m=N_m+\frac{\Delta t}{2}{\bar{u}}_{i,j}\frac{\partial N_m}{\partial x}+\frac{\Delta t}{2}{\bar{w}}_{i,j}\frac{\partial N_m}{\partial z}$ , (7)

where ${\bar{u}}_{i,j}$ and ${\bar{w}}_{i,j}$ are respectively the average velocity in the x and z directions in node $i,j$ and N_m is the m^th expansion function for node $i,j$ .

For each node the variables u*, w*, u, w, p and T are expanded in a series as follows:

$\begin{array}{l}{u}^{\ast }={\displaystyle \underset{m=1}{\overset{M}{\sum }}{u}_m^{\ast }}{N}_m;w^{\ast }={\displaystyle \underset{m=1}{\overset{M}{\sum }}{w}_m^{\ast }}{N}_m;u={\displaystyle \underset{m=1}{\overset{M}{\sum }}{u}_m}{N}_m;\\ w={\displaystyle \underset{m=1}{\overset{M}{\sum }}{w}_m}{N}_m;p={\displaystyle \underset{m=1}{\overset{M}{\sum }}{p}_m}{N}_m;T={\displaystyle \underset{m=1}{\overset{M}{\sum }}{T}_m}{N}_m;\end{array}$ (8)

where the parameters $u_m^{\ast }$ , $w_m^{\ast }$ , u_m, w_m, p_m and T_m are the coefficients of the variables associated with the m^th expansion function for that node.

For example, this process applied to the pressure equation, Equation (5), results in the following expression for node $i,j$ :

$\begin{array}{l}{\displaystyle \underset{n=1}{\overset{M}{\sum }}{p}_n\{\frac{\Delta z_{i,j}}{\Delta x_{i,j}}}[{\displaystyle \underset{-1}{\overset{1}{\int }}{\left(N_m\frac{\partial N_n}{\partial \xi }\right)}_{\xi =-1}^{\xi =1}\text{d}\eta }-{\displaystyle \underset{-1}{\overset{1}{\int }}{\displaystyle \underset{-1}{\overset{1}{\int }}\frac{\partial N_m}{\partial \xi }}}\frac{\partial N_n}{\partial \xi }\text{d}\xi \text{d}\eta ]\\ -\frac{\Delta x_{i,j}}{\Delta z_{i,j}}\left[{\displaystyle \underset{-1}{\overset{1}{\int }}{\left(N_m\frac{\partial N_n}{\partial \eta }\right)}_{\eta =-1}^{\eta =1}\text{d}\xi }-{\displaystyle \underset{-1}{\overset{1}{\int }}{\displaystyle \underset{-1}{\overset{1}{\int }}\frac{\partial N_m}{\partial \eta }}}\frac{\partial N_n}{\partial \eta }\text{d}\xi \text{d}\eta \right]\}\\ =\frac{\rho \Delta z_{i,j}}{2\Delta t}{\displaystyle \underset{n=1}{\overset{M}{\sum }}{u}_n^{\ast }}{\displaystyle \underset{-1}{\overset{1}{\int }}{\displaystyle \underset{-1}{\overset{1}{\int }}{P}_m\frac{\partial N_n}{\partial \xi }\text{d}\xi \text{d}\eta }}-\frac{\rho \Delta x_{i,j}}{2\Delta t}{\displaystyle \underset{n=1}{\overset{M}{\sum }}{w}_n^{\ast }}{\displaystyle \underset{-1}{\overset{1}{\int }}{\displaystyle \underset{-1}{\overset{1}{\int }}{P}_m}}\frac{\partial N_n}{\partial \xi }\text{d}\xi \text{d}\eta \text{for}m=1,\cdots ,M\end{array}$ (9)

where Dx_i,j and Dz_i,j are the length of node $i,j$ in the x and z direction respectively, x and h are the node locally space variables in the x and z directions respectively. This expression represents a system of M equations, where M is the total number of expansion functions used.

In matricial form this system of equations is written as follows:

$A_p^{i,j}p=b_p^{i,j}$ , (10)

where $A_p^{i,j}$ is the pressure equation matrix for node $i,j$ , which contains all the right hand side terms of Equation (9), p is the vector of the pressure expansion coefficients, and $b_p^{i,j}$ is a vector for node $i,j$ , which contains the left hand side terms of Equation (9). For Equation (1), Equation (3), and Equation (6) the coefficient matrices and the left hand side vectors are functions of the velocity expansion coefficients.

After applying this process for all equations at every node of the mesh, a system of equations, which must be solved interactively at each time step to calculate the expansion coefficients, was obtained.

4. Classic Expansion Functions

In the classical finite element method the coefficients of expansion of the variables are identified in specific points of the mesh. Figure 1 shows a typical one-dimensional element of dimension Dx with linear, quadratic and cubic expansion functions.

5. Hierarchical Functions

In the classic finite element method the expansion coefficients are identified with the variables at specified locations. This identification has been widely followed in the finite element literature and has the merit of assigning a “physical” meaning to these coefficients. There is, however, a great disadvantage in this practice. If it is desired to change the order of the expansion it is necessary to restart the problem due to the complete change of all expansion functions involved (Zienkiewicz and Morgan [1] ).

In the case of the hierarchical expansion functions the expansion coefficients are not identified with the physical variables at specific points of the mesh. In this case, the coefficients are associated with the expansion functions, which are adjusted in the rectangular nodes, defining corner, side and area functions. This

association allows starting the solution of a problem with a linear expansion and, if necessary, during the solution process, adding new shape functions to increase the order of expansion, and obtaining a more accurate solution. In this case, the expansion functions do not change by adding or deleting new expansion functions to change the order of the expansion. Thus, it is not necessary to recalculate the matrices. Changing the order of expansion without the need to restarting the problem turns out to be the great advantage of this method.

The hierarchical expansion functions are based on Legendre polynomials, which form a set of functions with orthogonality properties. Due to the association of two functions in each direction to form the two-dimensional expansion functions, complete orthogonality is impossible to achieve, but the main diagonal terms in the matrices are dominant.

Figure 2 shows four nodes and the parameters associated with their corners, sides, and areas. In this Figure ${\phi }_{i,j}^c$ are the corner parameters associated with the corner (linear) expansion functions, ${\phi }_{i,j,g}^x$ and ${\phi }_{i,j,g}^z$ are the side parameters associated with the side expansion functions, and ${\phi }_{i,j,g}^A$ are the area parameters associated with the area expansion functions. The order of the approximation for the side and area expansion functions is up to desired or necessary degree.

Each corner parameter is connected to four elements through four different expansion functions and each side parameters is connected to two elements by two different functions. The area parameters belong to a single element only. The association of the corners and side expansion coefficients with adjacent elements guarantees the singularity and the continuity of the solution at all nodes boundaries.

It should be observed that it is very easy to increase the order of expansion locally to achieve a refinement in the region where the parameters vary most rapidly and where the approximation is therefore prone to the largest errors.

6. Results

In order to validate the numerical method proposed in this work, some well-known problems of the literature were simulated. This paper presents the results of a backward-facing step problem, and additional results can be found in Sabundjian [6] .

For the accomplishment of this analysis the experimental results obtained by Denham and Patrick [7] were used. Figure 3 shows the geometry of the problem. A course mesh with 25 × 12 elements and expansion orders 1, 2 and 3 were used.

Figure 4 and Figure 5 present: The profile of velocity along a flow through a step and the graphical results for the velocity for the case of third order expansion, respectively.

Table 1 presents the mean errors between the calculated and the experimental results. As the expansion order increases, the calculated results are closer to the experimental results. However, the differences between the numerical and experimental results increase with the distance from the inlet.

7. Conclusion

Based on the results obtained in the solution of the proposed problems, it is possible to conclude that the method of hierarchical expansion is suitable for solution of incompressible fluid problems in two dimensions. The calculated results are in good agreement with the experimental results. The great advantage of the hierarchical expansion method is the capacity to adapt the expansion order to the necessary or desired value during the flow calculation, without the necessity to restart the problem, as it happens in the conventional finite element method, or in the method of mesh refinement.

Cite this paper

Sabundjian, G., das Neves Conti, T. and Cabral, E.L.L. (2018) Hierarchical Expansion Method in the Solution of the Navier-Stokes Equations for Incompressible Fluids in Laminar Two-Dimensional Flow. Energy and Power Engineering, 10, 1-9. https://doi.org/10.4236/epe.2018.101001

References