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  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  . 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  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.  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  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:
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 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:
where T is the temperature, k is the thermal conductivity and cp 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, Pm, is based on the Petrov-Galerkin formulation and it is given, according to Brooks and Hughes  , by:
where and are respectively the average velocity in the x and z directions in node and Nm is the mth expansion function for node .
For each node the variables u*, w*, u, w, p and T are expanded in a series as follows:
where the parameters , , um, wm, pm and Tm are the coefficients of the variables associated with the mth expansion function for that node.
For example, this process applied to the pressure equation, Equation (5), results in the following expression for node :
where Dxi,j and Dzi,j are the length of node 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:
where is the pressure equation matrix for node , which contains all the right hand side terms of Equation (9), p is the vector of the pressure expansion coefficients, and is a vector for node , 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  ).
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
Figure 1. One-dimensional elements and expansion functions (a) linear (b) quadratic e (c) cubic, according to Zinkiewics and Morgan (1983).
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 are the corner parameters associated with the corner (linear) expansion functions, and are the side parameters associated with the side expansion functions, and 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.
Figure 2. Four rectangular grid nodes and its associated parameters.
Figure 3. Backward-facing step geometry.
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  .
For the accomplishment of this analysis the experimental results obtained by Denham and Patrick  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. Profile of velocity along a flow through a step (“backward-facing step”) third order expansion (degree = 3).
Figure 5. Calculated results (third order expansion) and experimental results from Denhan and Patrick experiment with Reynolds = 73 .
Table 1. Mean error between the calculated and experimental results for the velocity.
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.
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.
 Brooks, N.A. and Hughes, T.J.R. (1982) Streamline Upwind Petrov-Galerkin Formulation for Convection-Dominated Flows with Particular Emphasis on the Compressible Navier-Stokes Equations. Computer Methods in Applied Mechanics and Engineering, 32, 199-259.
 Chen, Z.J., Li, Z.Y., Xie, W.L. and Wu, X.H. (2017) A Two-Level Variational Multiscale Meshless Local Petrov-Galerkin (VMS-MLPG) Method for Convection-Diffusion Problems with Large Peclet Number. Computers & Fluids, 1, 1-10.
 Sabundjian, G. (1999) Hierarchical Expansion Method in the Solution of the Navier-Stokes Equations for Incompressible Fluids in Laminar Two-Dimension Flow. Ph.D. Thesis, Mechanical Engineer Department, University of São Paulo, Brazil (in Portu-guese).