Control volume finite element method (CVFEM)

1.1 Introduction

Fluid flow has several applications in engineering and nature. Mathematically, real flows are governed by a set of nonlinear partial differential equations in complex geometry. So, suitable solutions can be obtained through numerical techniques such as the finite difference method, the finite volume method (FVM), and the finite element method (FEM). In the past decade the FEM has been developed for use in the area of computational fluid dynamics; this method has now become a powerful method to simulate complex geometry. However, the FVM is applied most in calculating fluid flows. The control volume finite element method (CVFEM) combines interesting characteristics from both the FVM and FEM. The CVFEM was presented by Baliga and Patankar [1,2] using linear triangular finite elements and by Raw and Schneider [3] using linear quadrilateral elements. Several authors have improved the CVFEM from then to now. Raw et al. [4] applied a nine-nodded element to solve heat conduction problems. Banaszek [5] compared the Galerkin and CVFEM methods in diffusion problems using six-nodded and nine-nodded elements. Campos Silva et al. [6] developed a computational program using nine-nodded finite elements based on a control volume formulation to simulate two-dimensional transient, incompressible, viscous fluid flows. Campos Silva and Moura [7] and Campos Silva [8] presented results for fluid flow problems. The CVFEM combines the flexibility of FEMs to discretize complex geometry with the conservative formulation of the FVMs, in which the variables can be easily interpreted physically in terms of fluxes, forces, and sources. Saabas and Baliga [9,10] referenced a list of several works related to FVMs and CVFEMs. Voller [11] presented the application of CVFEM for fluids and solids. Sheikholeslami et al. [12] studied the problem of natural convection between a circular enclosure and a sinusoidal cylinder. They concluded that streamlines, isotherms, and the number, size, and formation of cells inside the enclosure strongly depend on the Rayleigh number, values of amplitude, and the number of undulations of the enclosure.

1.2 Discretization: Grid, Mesh, and Cloud

In general there are three ways to place node points into a domain [11].

1.2.1 Grid

A basic approach assigns the location of nodes using a structured grid where, in a two-dimensional domain, the location of a node is uniquely specified by a row and a column index (Fig. 1.1a). Although such a structured approach can lead to convenient and efficient discrete equations, it lacks flexibility in accommodating complex geometries or allowing for the local concentration of nodes in solution regions of particular interest.

1.2.2 Mesh

Geometric flexibility, usually at the expense of solution efficiency, can be added by using an unstructured mesh. Figure 1.1b shows an unstructured mesh of triangular elements. In two-dimensional domains triangular meshes are good selections because they can tessellate any planar surface. Note, however, that other choices of elements can be used in place of or in addition to triangular elements. The mesh can be used to determine the placement of the nodes. A common choice is to place the nodes at the vertices of the elements. In the case of triangles this allows for the approximation of a dependent variable over the element, by linear interpolation between the vertex nodes. Higher-order approximations can be arrived by using more nodes (e.g., placing nodes at midpoints) or alternative elements (e.g., quadrilaterals). When considering an unstructured mesh recognizing the following is important:

1. The quality of the numerical solution obtained is critically dependent on the mesh. For example, avoiding highly acute angles is a key quality requirement for a mesh of triangular elements. The generation of appropriate meshes for a given domain is a complex topic worthy of a monograph in its own right. Fortunately, for two-dimensional problems in particular, there is a significant range of commercial and free software that can be used to generate quality meshes.

2. The term unstructured is used to indicate a lack of a global structure that relates the position of all the nodes in the domain. In practice, however, a local structure—the region of support—listing the nodes connected to a given node i is required. Establishing, storing, and using this local data structure is one of the critical ingredients in using an unstructured mesh.

1.2.3 Cloud

The most flexible discretization is to simply populate the domain with node points that have no formal background grid or mesh connecting the nodes. Solution approaches based on this “meshless” form of discretization create local and structures, usually based on a “cloud” of neighboring nodes that fall within a given length scale of a given node i [13] (Fig. 1.1c).

1.3 Element and Interpolation Shape Functions

A building block of discretization is the triangular element (Fig. 1.2). For linear triangular elements the node points are placed at the vertices. In Fig. 1.2, the nodes, moving in a counterclockwise direction, are labeled 1, 2, and 3. Values of the dependent variable ϕ are calculated and stored at these node points.

In this way, values at an arbitrary point (x, y) within the element can be approximated with linear interpolation

≈ax+by+c,

where the constant coefficients a, b, and c satisfy the nodal relationships

i=axi+byi+c,i=1,2,3.

Equation (1.1) can be more conveniently written in terms of the shape function N1, N2, and N3, where

ixy=1Atnodei0Atallpointsonsideoppositenodei

i=13Nixy=1Ateverypointintheelement

such that, over the element the continuous unknown field can be expressed as the linear combination of the values at nodes i = 1,2,3:

xy≈∑i=13Nixyϕi.

With linear triangular elements a straightforward geometric derivation for the shape functions can be obtained. With reference to Fig. 1.2, observe that the area of the element is given by

123=121x1y11x1y11x1y1=12x2y3−x3y2−x1y3−y2+y1x3−x2

and the area of the subelements with vertices at points (p, 2, 3), (p, 3, 1), and (p, 1, 2), where p is an arbitrary and variable point in the element, are given by

p23=x2y3−x3y2−xpy3−y2+ypx3−x2Ap31=x3y1−x1y3−xpy1−y3+ypx1−x3Ap12=x1y2−x2y1−xpy2−y1+ypx2−x1.

Based on these definitions, it follows that the shape functions are given by

1=AP23/A123,N2=AP31/A123,N3=AP12/A123.

Note that, when point p coincides with node i (1, 2, or 3), the shape function Ni = 1, and when point p is anywhere on the element side opposite node i, the associated subelement area is zero, and, through Eq. (1.8), the shape function Ni = 0. Hence the shape functions defined by Eq. (1.8) satisfy the required condition in Eq. (1.3). Further, note that at any point p, the sum of the areas:

P23+AP31+AP12=A123

is such that the shape functions at (xp, yp) sum to unity. Hence the shape functions defined by Eq. (1.8) also satisfy the conditions in Eq. (1.4)....