COMSOL Blog

FEM vs. FVM

by Ed Fontes

November 29, 2018

There has always been a debate among engineers working with fluid flow about the suitability of finite element methods for CFD. Some engineers have a firm opinion about the superiority of finite volume methods compared to finite element methods. Is there a scientific support for this opinion? No, not in general. Different methods may be suitable for different problems. Let us see why.

Science, Technology, and Tradition

Finite element methods are widely used by the numerical analysis community to study numerical methods for fluid flow. There is a vast amount of work described in scientific publications regarding CFD and finite element methods as well as more recent work for nodal discontinuous Galerkin (DG) methods, which are finite element methods with discontinuous basis functions.

Commercial packages for CFD are traditionally based on finite volume methods. This is due to the fact that basically all of the larger commercial packages for CFD have the same ancestors. There is a vast amount of work and technology invested in these methods in the industry. Different methods have been implemented so that they can efficiently and accurately compute and integrate the fluxes for both structured (for example, hexahedrons; i.e., bricks) and unstructured (for example, tetrahedrons) meshes.

That said, there is no theoretical or practical support for the hypothesis that finite volume methods are superior to finite element methods for fluid flow. First, there are many different finite volume methods and finite element methods and some of these methods overlap. Second, the implementation of a method has a very large impact on the practical use of a software. We can say that package X seems to do a better job than package Y for a certain family of problems. However, this is not because one uses FEM and the other FVM. Let us explain this below.

Finite Elements Versus Finite Volumes: What’s Best?

Mathematical Models and Numerical Models

Let us start with a generic mathematical model:

(1)

$P\left( {\frac{\partial }{{\partial {\mathbf{x}}}},\frac{\partial }{{\partial t}}} \right)u = F$

(2)

$\begin{gathered}
u = f{\text{ initial condition}} \hfill \\
Bu = g{\text{ boundary condition}} \hfill \\
\end{gathered} $

Here, P denotes a differential operator, u the dependent variable (the solution variable), F a source, f a function that describes the initial condition, B an operator, and g a function at the boundary. The space coordinate \mathbf{x} represents all three directions (x, y, and z) in this case.

The mathematical model may describe a physical phenomenon such as fluid flow. In that case, the model may represent the conservation of momentum and mass in space and time. Luckily, the mathematical model stemming from such conservation laws, and complemented with initial values and adequate boundary conditions, is often well behaved. This means that there is a unique solution and the solution behaves continuously on the data for the problem; for example, on initial values, source terms, or boundary conditions.

Despite the fact that a problem has a unique and stable solution, it may be difficult or almost impossible to find such a solution analytically. That is, it may be difficult to find an analytic expression for the solution with operations that can be easily calculated (+, -, x, ÷). Instead, we have to formulate a numerical model, which approximates the mathematical model. The numerical model equations can then be solved using a numerical method implemented in a computer program.

Finite element and finite volume methods are numerical methods based on discretization of space of the model equations. The time discretization is usually done with some type of time-stepping scheme for ordinary differential equations. The mathematical model defined above gives the following numerical model:

(3)

${P_h}{u_h} = {F_h}$

(4)

$\begin{gathered}
{u_h} = {f_h}{\text{ initial condition}} \hfill \\
{B_h}{u_h} = {g_h}{\text{ boundary condition}} \hfill \\
\end{gathered} $

where h denotes a discretization parameter; for example, the mesh element or cell size in a finite element or finite volume method.

Note that the building blocks in a mesh are referred to as elements in finite element methods and cells in finite volume methods.

There are several sources of error. The truncation error, $\tau $, tells us how well the numerical model approximates the mathematical model:

(5)

$\tau = \left( {P-{P_h}} \right)u$

The order of accuracy of the numerical model tells us how fast the truncation error decreases with decreasing h. This means that the smaller the element or the cell size, the smaller the difference between the numerical and the mathematical model should become. A numerical model is consistent if the truncation error decreases with h.

The discretization error in the solutions is given as the difference between the exact solution and the numerical solution to the model equations:

(6)

$e = \left\| {u-{u_h}} \right\|$

The numerical method is said to converge if the numerical solution approaches the exact solution as h decreases:

(7)

$\left\| {u-{u_h}} \right\| \to 0{\text{ as }}h \to 0$

The order of accuracy of the discretization, p, tells us how fast the numerical solution converges to the exact solution with decreasing h.

(8)

$e\left( h \right) \leqslant {C_1}{h^p}$

The larger p is, the faster the approximation converges.

So, is there any inherent difference in the order of accuracy of finite element methods and the finite volume methods? By increasing the order of the base functions, we can theoretically achieve any order of accuracy with finite element methods (in practice, there are other limitations). The most common finite element methods are second- to third-order accurate, and finite volume methods are first- to second-order accurate.

Where Are the Similarities and Where Are the Differences?

Let us look at a flux balance equation, which forms the base in mathematical models for fluid flow:

(9)

\frac{\partial u}{\partial t} + \nabla \cdot \Gamma = F{\text{ in }}\Omega

In this equation, u denotes the conserved physical quantity, such as momentum or mass, and \Gamma denotes the flux of that quantity; for example, momentum flowing across a control surface per unit area and per unit time.

Finite element methods start by formulating an integral equation where the equations are weighted with test functions, \varphi, and averaging is done by integrating over the model domain:

(10)

$\int\limits_\Omega {\frac{{\partial u}}{{\partial t}}\varphi dV} + \int\limits_\Omega {\left( {\nabla \cdot \Gamma } \right)\varphi dV} = \int\limits_\Omega {F\varphi dV} $

However, before we continue, let us apply the divergence theorem to \[{\Gamma \varphi }\]. This gives the following expression:

(11)

\[\int\limits_\Omega {\nabla \cdot \left( {\Gamma \varphi } \right)dV = } \int\limits_{\partial \Omega } {\left( {\Gamma \varphi } \right) \cdot {\mathbf{n}}dS} \]

Here, \partial \Omega denotes the boundary of the domain \Omega and \mathbf{n} denotes the normal vector to the domain boundary. Integration by parts of the left-hand side in the equation above gives:

(12)

\[\int\limits_\Omega {\nabla \cdot \left( {\Gamma \varphi } \right)dV = } \int\limits_\Omega {\left( {\nabla \cdot \Gamma } \right)\varphi dV + } \int\limits_\Omega {\Gamma \cdot \nabla \varphi dV} \]

which means that from Eq. 11, we get:

(13)

\[\int\limits_\Omega {\left( {\nabla \cdot \Gamma } \right)\varphi dV + } \int\limits_\Omega {\Gamma \cdot \nabla \varphi dV} = \int\limits_{\partial \Omega } {\left( {\Gamma \varphi } \right) \cdot {\mathbf{n}}dS} \]

We now have the following expression:

(14)

\[\int\limits_\Omega {\left( {\nabla \cdot \Gamma } \right)\varphi dV} = -\int\limits_\Omega {\Gamma \cdot \nabla \varphi dV} + \int\limits_{\partial \Omega } {\left( {\Gamma \varphi } \right) \cdot {\mathbf{n}}dS} \]

We can use Eq. 14 above for the second term in Eq. 10. This is done in order to naturally include the flux boundary conditions in the integral equations. Down the road, in the numerical implementation, this also leads to the advantage that the flux vector does not have to be differentiable. It results in the equations used as the starting point in finite element methods:

(15)

$\int\limits_\Omega {\frac{{\partial u}}{{\partial t}}\varphi dV}-\int\limits_\Omega {\Gamma \cdot \nabla \varphi dV} {\text{ + }}\int\limits_{\partial \Omega } {\Gamma \cdot {\mathbf{n}}{\text{ }}\varphi dS} = \int\limits_\Omega {F\varphi dV} $

This is the so-called weak formulation. The third term on the left-hand side integrates the flux of u, \Gamma, over the boundary of the domain, \partial \Omega.

The weak formulation is only related to the physics model if it holds for a large variation of test functions \varphi. A common choice is to use polynomials, but they can also be other types of functions. One special case is the case of a constant test function; for example, that \varphi = 1. The last equation above, Eq. 15, then yields:

(16)

$\int\limits_\Omega {\frac{{\partial u}}{{\partial t}}dV} {\text{ + }}\int\limits_{\partial \Omega } {\Gamma \cdot {\mathbf{n}}{\text{ }}dS} = \int\limits_\Omega {F{\text{ }}dV} $

This relation is used as the starting point for finite volume methods.

So far, there is no difference between the finite element and finite volume methods. As we can see above, the formulation for finite volume methods, Eq. 16, is just a special case of the generic weak formulation used in finite element methods, Eq. 15.

The difference is found in the discretization of Eq. 15 and Eq. 16. The finite element method is obtained from picking a finite number of test functions \varphi = \varphi_h and requiring Eq. 15 to hold for all of them. The finite volume method is obtained by picking a finite number of control volumes \Omega = \Omega_h and requiring Eq. 16 to hold for all of them. If we use a triangulation as a basis for both these methods, Figures 1 and 2 show possible discretized forms of finite element and finite volume formulations, respectively.

If we look at the most common finite element methods, then the test functions are only nonzero in the vicinity of so-called nodes (local supported function). This means that the integrals only need to be computed over elements (here, triangles) in this vicinity. See, for example, the highlighted domain node and the elements in its vicinity colored gray in Figure 1 below. The contribution from the boundary fluxes, the third term on the left-hand side in Eq. 15, only needs to be included for elements that have a face (3D) or an edge (2D) at the boundary, since interelemental boundary contributions cancel out for continuous base functions. Figure 1 below shows how the equation may be formulated for internal elements in the domain (white, gray, and green) and for elements that have a face (3D) or edge (2D) on the boundary (light blue). For the highlighted node in the domain in Figure 1, only the surrounding gray elements give a contribution to the domain integral (over \Omega). For the highlighted node at the boundary, only the two adjacent light blue elements give a boundary flux contribution to the integral over \partial \Omega, while both these two light blue elements and the light green element give contributions to the domain integral (over \Omega).

A schematic of the domain elements and equations for an FEM model in COMSOL Multiphysics.
Figure 1. The domain contribution for internal elements (white and gray) and for elements that have a face (3D) or edge (2D) on the boundary. The basis function for the node in the middle of the gray hexagon has support in all surrounding elements; i.e., in all the gray elements. The basis functions for the nodes at the boundaries have support in the light blue elements but also in any element that has a node at the boundary; for example, the light green element. The integral for the flux only gets contributions from the elements that have an edge on the boundary (light blue elements).

The domain expression for the discretized flux, \Gamma_h, is obtained from the constitutive relations for the flux; for example, the convection-diffusion equation, (the diffusion term in the fluid flow equations is the viscous term for momentum transport). The boundary expressions for the flux in Eq. 15 are obtained from the boundary conditions for the specific model. These expressions are then included as a contribution for the boundary elements (light blue) in Figure 1 above.

If we instead look at a common finite volume method, the cell-centered method, each cell (triangle) is treated as an individual domain. The boundary term for the flux, the second term on the left-hand side in Eq. 16, is integrated for all cells — both internal cells and cells that have a face (3D) or edge (2D) on the boundary. The constitutive relation for the flux is used for the domain faces or edges, while the boundary conditions are used for faces or edges on the boundary; see Figure 2 below.

A schematic of the cell faces, cells, and equations for an FVM model.
Figure 2. The flux is integrated over all cell faces (3D) or edges (2D), both for internal cells and cells that have a face or edge on a boundary.

How do we then express u and \Gamma in the two different methods?

In finite element methods, often the same basis functions are used to approximate the solution as the test functions. As long as the approximation of the solution has a higher polynomial degree than zero and the first-order derivatives can be approximated, nothing special needs to be done to handle a flux given by convection and diffusion. The flux vector is also a local polynomial function.

In finite volume discretizations, on the other hand, the solution on the boundary is not well defined. The method only defines the value of the solution for each cell normally interpreted as the value at the center of the cell. A finite volume method therefore needs to be completed with some sort of reconstruction method to be useful. Typically, a local interpolation method is used where neighboring cell values are taken into account; see Figure 3 for an example. To obtain a higher-order interpolation of the solution and the flux, more cell values need to be taken into account. This is not only complicated but also leads to a less local method.

A schematic of a cell-centered finite volume method used to compare FEM vs. FVM.
Figure 3. In a cell-centered finite volume method, the flux vector is constructed by interpolation between points centered in the cell.

Depending on the basis functions used in a finite element method and the type of construction of the flux used in a finite volume method, different accuracies can be achieved. A coarse mesh with a second-order accuracy method can yield a more accurate solution than a finer mesh with a first-order accuracy method.

Linear test functions and basis functions for the finite element method normally lead to methods that are second-order accurate. Finite elements have large flexibility when it comes to the discretization. Using quadratic basis functions, for example, is rather straightforward. There is no need for reconstruction or interpolation of the solution either. The method is very symmetric in nature and flux boundary conditions on boundaries can be imposed in a natural, direct way.

A drawback with FEM is that for continuous test and basis functions, there is no local conservation property; only global conservation is guaranteed. In other words, only the net flux over the domain boundaries is guaranteed to be in balance. Another drawback is that there is no control of the local fluxes, which means that stabilizing the discretization for convection-dominated flow is not straightforward. Stabilizing means, in this case, removing unphysical oscillations that are artifacts of the discretization. Both local conservation and stabilization for convection-dominated flow can be addressed by modifying the weak formulation, either directly or indirectly, by modifying the test functions. However, these methods can be expensive in terms of computation.

As mentioned earlier, finite volume methods correspond to piecewise constant finite element base functions, possibly with a higher-order interpolation scheme for the fluxes. This leads to methods that are first- or second-order accurate. The local statement for the finite volume method leads to local conservation, which is an attractive feature of the method. It means that the net flux for each cell is guaranteed to be in balance. It also leads to natural and direct ways to stabilize the discretization for flow problems that are convection dominated. So-called upwind stabilization and other stabilizations are naturally obtained by modifying the fluxes on the intercell boundaries. Upwinding gives nonsymmetry in the discretization in the direction of the convective flux.

The finite element method has the benefit of being able to formulate methods for basis functions of different orders. Higher orders for the basis functions give higher-order, accurate methods, which have the important benefit of being able to improve the accuracy for a given mesh. The finite volume uses zeroth-order base functions but can use higher-order interpolation schemes for the flux instead — also yielding improved accuracy. When using higher-order methods, the resulting system becomes larger and the solution time increases for the same mesh. However, this also comes with a higher order of accuracy. So, when we compare performance, we have to compare it for a given accuracy. Measuring the CPU time and the memory required to solve a fluid flow problem to the same accuracy, with different methods, is the correct way of comparing the performance for these different methods, not the number of cells or elements.

Future Finite Element Methods

At COMSOL, we work mainly with finite element methods for CFD, since this is where we have our expertise. During the last 15 years, major achievements have been made by the research community in the development of finite element methods with discontinuous test and basis functions. These are the DG methods mentioned in the introduction. The test functions in these methods are local to each element and the weak equations hold for each element. Local conservation is obtained automatically and higher-order discretizations are straightforward. There is no need for reconstruction of the solution either, except for the zeroth-order method, where it is equivalent to the finite volume method. In addition, the local flux on the element boundary is a natural part of the formulation, which means that stabilization is straightforward. A drawback with DG methods is the relatively large amount of extra degrees of freedoms introduced. This has lead to the study of so-called hybrid DG methods, where the discretization is leaner and yields fewer degrees of freedom.

Concluding Thoughts on FEM vs. FVM

As we have shown here, there are pros and cons in both finite volume and finite element methods. There are many more aspects that are important for efficient computing of fluid flow problems on a large scale: efficient and stable handling of time-stepping; efficient implementation of both implicit and, to some extent, explicit time-stepping; solving large systems of linear equations; and so on. There is still a lot of work to be done and there are still many possibilities to advance different methods and technologies!

We at COMSOL are constantly working to provide the best and latest technology using finite element methods. We are also committed to developing the best methods in general — not only for multiphysics, but also for “single” physics applications such as CFD.