Scientific journal
European Journal of Natural History
ISSN 2073-4972


Mohammad Reza Akhavan Khaleghi 1
1 The Office of Counseling and Research Fluid Engineering and Aerodynamic
1. Zienkiewicz O.C. and Morgan K. Finite Elements and Approximation, John Wiley & Sons, 1983.
2. Zienkiewicz O.C. and Taylor R.L. The Finite Element Method: Volume 1, the Basis, Butterworth-Heinemann, 2001.
3. Laney C.B. Computational Gasdynamics. Publisher: Cambridge University Press. Pub. Date: June 28, 1998.
4. Hughes T.J.R, Cottrell J.A, Bazilevs Y. “Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement”. Computer Methods in Applied Mechanics and Engineering, 194(39–41), Р. 4135–4195, 2005.
5. Hoffmann K.A. and Chiang, S.T. Computational Fluid Dynamics, Vol. I and II, 3rd edition, Engineering Education System (1998).
6. Deng J., Chen F., Li X., Hu C., Tong W., Yang Z. and Feng Y. Polynomial splines over hierarchical T-meshes. Graph. Models 70(4):76–86, 2008.

First, I should mention that this is a basic package and is not limited to CFD, it can also be used for other problems.

Finite Element Method (FEM) is a powerful numerical method which has been used successfully for the solution of the existing problems in various scientific and engineering fields such as its application in CFD. Many algorithms have been expressed based on FEM, but none has been used in popular CFD software. In this section, full monopoly is according to Finite Volume Method (FVM) due to better efficiency and adaptability with the physics of problems in comparison with FEM. It doesn’t seem that FEM could compete with FVM unless it was fundamentally changed. In this paper, I am going to show those changes and its result will be a powerful method which has much better performance in all subjects in comparison with FVM and other computational method, I called it Reduced Finite Element Method (RFEM).

The general form of a linear differential equation with boundary conditions can be shown as follows [Zienkiewicz and Morgan (1983); Zienkiewicz and Taylor (2001)]:


mohamad02.wmf (1)

The operators mohamad03.wmf and mohamad04.wmf in equations (1) can be zero-order, odd-order, even-order or a combination of two or all three


mohamad06.wmf (2)

The weighted residual relationship for these relationships on any element is


mohamad08.wmf (3)


mohamad09.wmf (4)

By inserting (2) in (3) we have


mohamad11.wmf (5)



The approximate relation of the field also is equal to

mohamad14.wmf (6)

By inserting (6) in (5) we have



mohamad17.wmf (7)


Together all the elements a comprehensive system of equations is obtained which can be written quite generally as

mohamad19.wmf (8)


mohamad20.wmf (9)

Finally by using of equation (7) can write


mohamad22.wmf (10)


And this is the general form of approximation to a differential equation by the finite element method.

New Formulation for Finite Element Method (NFFEM)

Here I will introduce a new formulation for finite element method which its performance is much better than all conventional algorithms of finite element method in CFD. In NFFEM relationship (9) is written as follows:

mohamad24.wmf (11)

And each of the matrix coefficients are

mohamad25.wmf (12)

Where mohamad26.wmf become

mohamad27.wmf (13)

For mohamad28.wmf we have

mohamad29.wmf (14)

For mohamad30.wmf we have

mohamad31.wmf (15)

And for mohamad32.wmf we have

mohamad33.wmf (16)

These relations are used for boundary integrals the same way, these are relations completely of NFFEM.

In the equations (13) to (16), by choosing b = 1 full upwind difference scheme (FUDS) and by choosing b = 0 central difference scheme (CDS) is obtained (Flux Vector Splitting Methods (FVSM) can also use with NFFEM, in this form b = 1 is considered and equations are written separately for positive and negative-terms after that the sum of the equations is done). These relationships are also valid for non-linear operators (mohamad34.wmf), which leads to mohamad35.wmf. In the equations presented above, mohamad36.wmf is given by

mohamad37.wmf (17)

For the central difference scheme the following equation can also be used instead of mohamad38.wmf in equations (13) to (16)

mohamad39.wmf (18)

Where mohamad40.wmf is given by

mohamad41.wmf (19)

In this form, the number of participating elements in the equation of node i are limited to two elements in 1D, 2D and 3D. The mohamad43.wmf, mohamad44.wmf and mohamad45.wmf are equal to

mohamad46.wmf (20)



mohamad48.wmf (21)


And αi is the sign of the unknown variable or derivations coefficients (any term that changes the sign of the matrix coefficients)

mohamad50.wmf mohamad51.wmf mohamad52.wmf (22)

Equation (17) can be written with βi and – βi instead of the mohamad53.wmf that its result will be mohamad54.wmf, depending on the sign of βi that we have chosen, the mohamad55.wmf gives the forward or backward differencing scheme.

Superscript k, l and m represent the terms are located on lines k, l and m (The node i in three dimensions is located on three lines, line k at x direction, line l at y direction, and line m at z direction, see figures (1) and (2)).

For hierarchical shape functions, the shape functions dependent to the sides and element as shown in figures (3) and (4) attributable to the points (these points as mohamad56.wmf, mohamad57.wmf and mohamad58.wmf to calculate mohamad59.wmf also to determine lines k, l and m are used).

When wind direction is constant before and after discontinuity (shock), as an alternative (or as another option for all position), mohamad60.wmf can be used instead of mohamad61.wmf, where mohamad62.wmf is

mohamad63.wmf mohamad64.wmf

mohamad65.wmf. (23)



mohamad67.wmf (24)


Because βi is zero for internal nodes, there is no need to calculate mohamad69.wmf, however it can be calculated by the following equation


mohamad71.wmf (25)


Note that NFFEM only for quadrilateral and hexahedron elements can be used. Solving the system of equations of the NFFEM can be performed simply by element by element solving method and line by line sweeping method.

Reduced Finite Element Method (RFEM)

A new method for the finite element formulation presented in the previous section, In this section, I’m going to limit the participating nodes in equation for each node to the nodes that are located on lines k, l and m the node. For this purpose, I introduce reduced elements, characteristics of these elements is as follows:

1 – The different elements are used for different directions of operator, see figures (5) and (6).

2 – For each direction of operator, element only in the same direction has DOF (for one direction of shape function is used p – degree function and for other directions is used zero degree function, see figures (5), (6), (7) and (8)).

3 – For equation of node i, element only on the lines k, l and m the same node has the node, see figures (7) and (8).

In RFEM relationship (11) is written as follows:

mohamad74.wmf mohamad75.wmf (26)

And each of the matrix coefficients are

mohamad76.wmf (27)


mohamad77.wmf (28)

For mohamad78.wmf, mohamad79.wmf and mohamad80.wmf we have


mohamad82.wmf (29)


For mohamad84.wmf, mohamad85.wmf and mohamad86.wmf we have


mohamad88.wmf (30)


For mohamad90.wmf, mohamad91.wmf and mohamad92.wmf we have


mohamad94.wmf (31)


Note: for the mohamad96.wmf operator, unknown variable is written in all directions and each direction is no belongs to particular unknown variable. And for mohamad97.wmf we have

mohamad98.wmf (32)


Fig. 1. Lines k and l on two fourth-degree elements in two dimensions


Fig. 2. Lines k, l and m on a quadratic element in three dimensions


Fig. 3. Lines k and l also the places that β is calculated on for a hierarchical element in two dimensions


Fig. 4. Lines k, l and m also the places that β is calculated on for a fourth-degree hierarchical element in three dimensions


Fig. 5. A reduced standard quadratic element in two-dimensional, lines k and l also the participating nodes in the equation of node i for each of them


Fig. 6. A reduced hierarchical quadratic element in two-dimensional, lines k and l also the participating nodes in the equation of node i for each of them


Fig. 7. Nodes involved in the equation of node i for FEM and RFEM on several fourth-degree elements in two dimensions


Fig. 8. Nodes involved in the equation of node i for FEM and RFEM on two quadratic elements in three dimensions


Fig. 9. A reduced quadratic element in three-dimensional with two DOF for approximation the mixed derivatives ∂2ϕ/∂x∂y, ∂2ϕ/∂x∂z and ∂2ϕ/∂y∂z


Fig. 12. 1D advective equation ut + ux = 0


Fig. 13. Euler equations 1D based on sod’s shock tube problem (t = 0.2)


Fig. 14. Grid for flow around a cylinder


Fig. 15. Flow around a cylinder (steady-state Euler equations)


Fig. 16. Mach contours for flow around a cylinder (steady-state Euler equations)


Fig. 17. Grid for flow over a NACA 0012 airfoil


Fig. 18. Surface skin friction coefficient distributions on NACA 0012 airfoil (steady-state Navier-Stokes equations)


Figure 19. Surface pressure distribution for subcritical NACA 0012 airfoil (steady-state Euler equations)


Fig. 20. Solution domain and lines grid for flat plate flow


Fig. 21. Skin friction distributions for flat plate flow (steady-state Navier-Stokes equations)


Fig. 22. 1D convection equation with source term, auх = сos(x) on [0, π] and a = 1

For approximating the mixed derivatives we will need to more DOF for example; if in three dimensions, the mixed derivative be in two-direction, like following derivative

mohamad99.wmf (33)

We use the shape function with two DOF in x and y directions, see figure (9). This elements are used only for mixed operator.

When we use the hierarchical shape functions, for some equations have two unknown variable (ϕi and ai, see figure (6)), the additional equation for additional unknown can be written as follows:


RFEM can be used to isogeometric analysis as well (see the following sections and examples).

A new hybrid difference scheme (HDS) for liner shape functions

In this section, I will introduce a technique to eliminate oscillations for using the FEM and RFEM on liner shape functions, in this case b is written as follows:

mohamad102.wmf (35)

Where mohamad103.wmf is “Upwind Parameter” and is chosen in the range mohamad104.wmf. It is clear that there are many choices that can be used for. I used the following equation for it

mohamad105.wmf (36)


mohamad106.wmf (37)



mohamad108.wmf (38)


The performance of this scheme is much better than other second-order schemes that in the FEM are used for CFD, see figure (11).

A new shape functions for general problems

Equation (35) can be used only for the linear shape functions, in this section another technique is provided for eliminating oscillations that can be used for higher-order shape functions in any degree. This technique works directly on the shape functions (this technique is equivalent to creating new shape functions), for this purpose, the shape functions are written as following

mohamad110.wmf (39)

Where δi is Element Degree Limiter, Field Variable Limiter (EDL, FVL), and N(p) can be any function of p ≥ 2 (for example, the functions of Lagrangian, Serendipity and Bezier for standard shape functions and the functions of Legendre, Chebyshev, Fourier, etc. for hierarchical shape functions). mohamad112.wmf, mohamad113.wmf and mohamad114.wmf in equation (39) are reference functions, original function for them become

mohamad115.wmf mohamad116.wmf mohamad117.wmf (40)

Where mohamad118.wmf and mohamad119.wmf are liner shape functions.

The δi in equation (39) has two range:

1 – Using the δi as EDL by choosing it in the range mohamad120.wmf.

2 – Using the δi as FVL by choosing it in the range mohamad121.wmf.

Note: if the reference function is greater than the shape function the sign of EDL and FVL will change.

The EDL can be used only in the directions of the shape function (directions of operator) that βi ≠ 0 is for, and the FVL is used only in the directions of the shape function that βi = 0 is for or inverse, or the FVL and EDL can be added on all nodes with identical or non-identical values, and their relationship is written as follows:

mohamad122.wmf (41)

The value of the mohamad123.wmf and mohamad124.wmf depending on the application is, see examples. Note that δi is applied to shape functions as one-dimensional (separately for each direction of shape functions).

The following functions can also be used as reference function on the Lagrangian, Serendipity and Bezier functions to make the shape functions of non-oscillatory

mohamad125.wmf mohamad126.wmf

mohamad127.wmf (42)

And for boundary elements become

mohamad128.wmf mohamad129.wmf

mohamad130.wmf mohamad131.wmf

mohamad132.wmf mohamad133.wmf (43)

Where e is a constant coefficient and N0, Nj and NL can be the Lagrangian, Serendipity or Bezier functions and degree them can be chosen 1 or p (for 1 – degree and p – degree the equations and results are different). mohamad134.wmf and mohamad135.wmf are a constant coefficients and are given by


mohamad137.wmf (44)

Where ξj is the location of the nodes (control points) and ω is the weight coefficient (optional). For p ≥ 3 the hierarchical shape functions with mohamad138.wmf can be used. e is chosen as, e = 1 – p for ω = ∞, e = – 1 for ω = 0 and e = 1/1 – p for ω = – ∞.

Another reference function is



mohamad141.wmf (45)

When αi = 1 equation (45) give backward difference approximation and when αi = – 1 give forward difference approximation. A reference function with high performance around discontinuities in derivative form is as follows:








mohamad149.wmf (46)





The difference between equation (46) and equations (40), (42) and (45) on a fourth-degree element around the discontinuity in figure (10).


Fig. 10. The difference between equation (46) and equations (40), (42) and (45) on three fourth-degree elements around the discontinuity with FUDS


Fig. 11. Liner functions for quadratic and cubic IGA functions

The r for higher-order shape functions is written as re (the same for all nodes lines k, l and m)


mohamad155.wmf (47)


Where ri is given by equation (38).

New isogeometric analysis and new functions

In this section, I’m going by reference functions presented in the previous section, I introduce a new isogeometric analysis. The isogeometric analysis [Hughes et al. (2005)] is done by the B-Spline and NURBS functions, although the reference functions presented above are applicable with B-Spline and NURBS functions, here I introduce a new Isogeometric Analysis Functions that are closer to the standard finite element method and more comfortable for the use in CFD and engineering, I made the following algorithm to make these functions on the Bezier functions

mohamad157.wmf mohamad158.wmf


mohamad160.wmf mohamad161.wmf (48)

Note that mohamad162.wmf are local and can add as hierarchical as well, for p = 2 these functions are equivalent to the B-Spline functions and for p = 3 are similar the PHT-spline functions work of [Deng et al. (2008)] and for more than a third degree are completely new. The equation (48) for first and last boundary elements become

mohamad163.wmf mohamad164.wmf


mohamad166.wmf mohamad167.wmf (49)

mohamad168.wmf mohamad169.wmf mohamad170.wmf

mohamad171.wmf mohamad172.wmf (50)

For IGA functions, the liner functions to use in equations (40), (42), (45) and (46) are a liner functions between two control points, see figure (10), while the p-order functions for using in equations (42) and (45) are the same IGA functions.

Solutions and Examples

In this section, some examples that have been solved by standard shape functions that are including functions of the Lagrangian and Bezier and hierarchical shape functions that are including functions of the Legendre, Chebyshev (first and second kind) and Fourier Sine Series are presented. It explained that the volume examples were solved to test RFEM is very high and here are just a few of them select and presented. Examples include 1D advective equation, 1D and 2D Euler equations, 2D Navier – Stokes equations and1D convection equation. Here details of the solution of these equations due to the length of paper and be less important are not presented and only the results are shown. Note: All examples are RFEM (the results for the NFFEM is almost identical). I used the infinite elements for non-solid boundaries in all examples were solved.


The best criterion for evaluate a numerical method is the results of the method and as can be seen in the examples were solved, the result of RFEM comparable to the best results were obtained by other methods is (for these examples, because of the small size of the elements, the difference between the results from different methods can not be seen) and given that the system of equations resulting from this approach similar (in terms of density) Finite Difference Method (FDM) is its efficiency can be compared with finite difference methods (although the use of the Legendre shape functions gives an efficiency much higher than FDM) so do not think other methods that are used in CFD be able to compete with it. Also, due to the similarity of relations many of the techniques used in the FDM can be used to RFEM [Hoffmann, and Chiang (1998)].

The work is submitted to the International Scientific Conference “Engineering science and modern manufacture”, France, October, 18–25, 2015, came to the editorial office оn 13.08.2015.