THE FINITE ELEMENT METHOD (FEM): AN APPLICATION TO FLUID MECHANICS AND HEAT TRANSFER

In this paper the finite element method (FEM) is used to solve three problems that are of the paramount importance in Chemical Engineering. The first problem is related with the bidimensional flow of an ideal fluid around a cylindrical body, and the objective is to determine the velocity distribution of the flow. To model the flow, the potential formulation is used to obtain an analytical solution, and then, the approximated solution obtained by using FEM is compared with the analytical solution. From this comparison, it is deduced that both solutions have a good agreement. The second problem is the calculation of the temperature profile in a two-dimensional body with specified boundary conditions. This problem is modeled by the two-dimensional Laplace equation, and from the problem data and using variables separation, an analytical solution was obtained. Then, FEM was used to obtain an approximate solution and compared with analytical ones. Besides, from this comparison, it is concluded that both solutions agree. Finally, in the third problem the temperature distribution in a bidimensional body with internal heat generation is studied. This problem is modeled by Poisson equation in two dimensions, but due to the boundary conditions and the complications that arise by adding some heat sources in the final FEM discretization, the problem does not have an analytical solution. However, the analysis of FEM solution indicates that this solution is correct.


Introduction
For solving complex engineering problems, it is necessary to have methods that be computationally effective. Ideally, an effective computational method should have the following features: 1. It should have a sound mathematical as well as physical basis (i.e., yield convergent solutions and be applicable to practical problems).
2. It should not have limitations with regard to the geometry, the physical composition of the domain, or the nature of the 'loading'.
3. The formulative procedure should be independent of the shape of the domain and the specificity of the boundary conditions. 4. The method should be flexible enough to allow different degrees of approximation without reformulating the entire problem. 5. It should involve a systematic procedure that can be automated for use on a computer.
The finite element method is a technique in which a given domain is represented as a collection of simple domains called finite elements, so that it is possible to systematically construct the approximation functions needed in a variational or weighted-residual approximation of the solution of a problem in each element (Hughes [1], Lewis et al. [2], Zienkiewicz-Taylor [3,4], Ciarlet [5][6][7]). Thus, the finite element method differs from the traditional Rayleigh-Ritz, Galerkin, least squares, colocation, and other weighted-residual methods in the manner in which the approximation functions are constructed. But this difference is responsible for the following three basic characteristics of the finite element method: 1. Division of whole into parts, which allows representation of geometrically complex domains as collections of geometrically simple domains that enable a systematic derivation of the approximation functions.
2. Derivation of approximation functions over each element; the approximation functions are often algebraic polynomials that are derived using interpolation theory.
3. Assembly of elements, which is based on continuity of the solution and balance of internal fluxes; the assemblage of elements represents a discrete analog of the original domain, and the associated system of algebraic equations represents a numerical analog of the mathematical model of the problem being analyzed.
These three features, which constitute three major steps of the FEM formulation, are already closely related. The geometry of the elements used to represent the domain of a problem should be such that the approximation functions can be uniquely derived. The approximation functions depend not only on the geometry, but also on the number and location of points, called nodes, in the element and the quantities to be interpolated (e.g., solution, or solution and its derivatives). Once the approximation functions have been derived, the procedure to obtain algebraic relations among the unknown coefficients (which give the values of the solutions at the nodes of finite elements) is exactly the same as that used in the Rayleigh-Ritz and weightedresiduals (Gelfand-Fomin [8], Weinstock [9]). The FEM not only overcomes the shortcomings of the traditional variational methods, but it is also endowed with the features of an effective computational technique.

Mathematical Preliminaries
To understand how a differential equation set in a given bounded domain of the plane, is modeled by the use of the FEM, consider the following elliptic partial differential equation: where Ω is a bounded domain of the plane; and the unknown u(x, y) are functions defined on Ω. The specified boundary conditions are a combination of and its normal derivative on the boundary: , (2) where is called the residual. What is looked for is the best approximation of u in the class of continuous piecewise polynomials.
Therefore, the equation for is tested against all possible functions of that class. Testing formally means to multiply the residual against any function and integrate, i.e., determine such that: for all possible . The functions are usually called test functions. The Eq. (3) is integrated by using Green's formula (Evans [13]). Therefore, should satisfy: , where is the boundary of Ω, and ds is the arclenght differential on the boundary. Note that the integrals of this formulation are well-defined even if ? This gives a linear system KС = F, where the matrix K and the right-hand vector F contains integrals in terms of functions φ i , φ j and the coefficients defining the problem: c, a, f, q and g. The solution vector C contains the expansion coefficients of , which are also the values of at each node.

Computational aspects of FEM
From previous exposition, it is deduced that for the FEM applications, computer programs that help in the many calculations involved in such method are needed. The following are the main steps that are necessary for solving through FEM, physical problems that are governed by partial differential equations (case of the present work): 1. Preprocessing. In this step the following is considered: • Meshing: the bidimensional region Ω is divided in triangles. It is necessary to have a computer program for this task.
• Geometrical and physical data related with the problem at hand.
• The assemble of triangular element equations for obtaining the final system of algebraic equations. For linear problems, this system is linear.
2. Processing. This consists in the solution of equations system for obtaining the nodal values of the scalar quantity that is being approximated. To solve the problems of this paper, a Matlab® based program was written that implements the main steps of FEM.

Example 1 4.1.1. Statement
Determine the velocity distribution for the irrotational flow of an ideal fluid around a cylindrical body, 40 mm in diameter, centered between two parallel walls which are 80 mm apart. The fluid has a uniform velocity of 40 mm/s at a location far removed from the cylinder (see Figure 1).

Mathematical model
For a dimensional and incompressible flow, the continuity equation (Cengel-Cimbala [14], Reddy [15]) is written as: , (6) where and are the vertical and horizontal components of the velocity V of the fluid. The flow is called irrotational, if =0, which implies: .
The irrotational flow of an ideal fluid may be formulated in terms of a stream function or in terms of a potential velocity. The potential formulation is used here.
The normal component to fixed boundary is: In polar coordinates, the function ϕ(x, y) satisfies the following BVP (Edwards [16]): (14) Using the separation variables method, the solution of BVP (14) (Asmar [17]) is: (15) It is possible, from (15), to derive expressions for u, ν, u r , u θ , r, θ being the polar coordinates, R the cylinder radius, and U ∞ is the velocity of approximation of the fluid to cylinder.

Solution
The triangular mesh that was used for the solution of the problem is shown in Figure 2. The mesh has 445 nodes and 64 triangular elements. The boundary conditions of the problem are shown in Figure 3. The symmetry of the geometrical configuration makes possible to analyze a cylinder quadrant.
A form to validate the solution obtained from FEM is to compare its solution with the value in the left border of Figure 2   These values are close to the value . Therefore, the FEM solution is a good approximation to the solution of the considered problem.

Example 2 4.2.1. Statement
Use the FEM to calculate the temperature distribution in the dimensional body shown in Figure 4 and compare its results with the analytical solution.

Mathematical model
The mathematical model of the problem is the bidimensional Laplace equation , where , is the temperature. The boundary conditions consist of specified temperature in inferior border and insulation in superior border; the left border is subjected to a heat flux, and the right, to convection with specified surrounding temperature.

Solution
By using the separation of variable method (Asmar [17]), the analytical solution of the problem is: (16) Now, the solution by FEM is presented. Figure 5 shows the triangular mesh of region in Figure 4. The mesh has 121 nodes and 200 triangular elements. In Table 1 the analytical and FEM solution for 18 nodes are shown. This table shows the good agreement between FEM and analytical solution. In particular, the FEM solution reproduces the temperatures specified in the nodes.

Example 3 4.3.1. Statement
Use the FEM to calculate the temperature distribution in the dimensional body shown in Figure 6. The dimensions of the body are 5×12 m.

Mathematical model
The partial differential equation that models the problem is the Poisson equation , where T = T (x, y), and Q is the heat generation in units of E/m 3 . The boundary conditions are: specified heat flux in the straight sides of left side and specified temperature in the curve part of the left side; the superior side is insulated; the inferior side has a specified temperature; and the right side has convection with surroundings. The value of Q will be specified later (see Solution).

Solution
The figure 7 shows the triangular mesh for the bidimensional region of Figure 6. It has 36 nodes and 48 triangular elements. After discretization, nodal heat sources of Q= 4000 W/m 3 are added to nodes 6, 18 and 31; and to elements 15 and 34. This problem does not have an analytical solution.
In Table 2 the results obtained by FEM are presented. As it is expected, the FEM solution reproduces the correct The maximum temperature is located at node 31 with value 266.5 °C that corresponds to a node with a specified heat source. The above considerations permit to conclude that the FEM solution is correct.

Conclusions
The finite element method is a numerical procedure for approximating the solution of complex problems governed by differential equations that are partial differential equations in the presented examples of this paper. Three important problems in chemical engineering were presented and solved by means of FEM.
Two problems solved by FEM have analytical solutions, and both solutions were compared. The FEM solution has a good agreement with the analytical solution. The Example 3, a complex problem of heat transfer, does not have an analytical solution, but based on the analysis of the solution, it may be concluded that the FEM solution is correct.