The equations of fluid mechanics – the Navier-Stokes equations, are solvable analytically for only a limited number of flows under certain assumptions. The known solutions are extremely useful in helping to understand fluid flow but rarely can they be used directly in engineering analysis or design. Therefore it is obligatory to use other approaches in order to solve practical problems.
In experimental analysis, the problem is that many flows require several dimensionless parameters for their specification and it may be impossible to set up an experiment which correctly scales the actual flow. Examples are flows around aircraft or ships. In order to achieve the same Reynolds number with smaller models, fluid velocity has to be increased. For aircraft, this may give too high a Mach number if the same fluid (air) is used; one tries to find a fluid which allows matching of both parameters. For ships, the issue is to match both the Reynolds and Froude numbers, which is nearly impossible.
Today, due to increased computational power of advance computers, interest in numerical techniques increased dramatically. Solution of the equations of fluid mechanics on computers has become so important that it now occupies the attention of perhaps a third of all researchers in fluid mechanics and the proportion is still increasing. This field is known as computational fluid dynamics (CFD). Contained within it are many subspecialties. We shall discuss only a small subset of methods for solving the equations describing fluid flow and related phenomena.
What is CFD?
Flows and related phenomena can be described by partial differential equations, which cannot be solved analytically except in special cases. To obtain an approximate solution numerically, we have to use a discretization method which approximates the differential equations by a system of algebraic equations, which can then be solved on a computer. The approximations are applied to small domains in space and/or time so the numerical solution provides results at discrete locations in space and time. Much as the accuracy of experimental data depends on the quality of the tools used, the accuracy of numerical solutions is dependent on the quality of discretizations used. Contained within the broad field of computational fluid dynamics are activities that cover the range from the automation of well-established engineering design methods to the use of detailed solutions of the Navier-Stokes equations as substitutes for experimental research into the nature of complex flows. At one end, one can purchase design packages for pipe systems that solve problems in a few seconds or minutes on personal computers or workstations. On the other, there are codes that may require hundreds of hours on the largest super-computers.
Finite Difference Method:
This is the oldest method for numerical solution of PDE’s, believed to have been introduced by Euler in the 18th century. It is also the easiest method to use for simple geometries. The starting point is the conservation equation in differential form. The solution domain is covered by a grid. At each grid point, the differential equation is approximated by replacing the partial derivatives by approximations in terms of the nodal values of the functions. The result is one algebraic equation per grid node, in which the variable value at that and a certain number of neighbor nodes appear as unknowns.
In principle, the FD method can be applied to any grid type. However, in all applications of the FD method known to the authors, it has been applied to structured grids. The grid lines serve as local coordinate lines. Taylor series expansion or polynomial fitting is used to obtain approximations to the first and second derivatives of the variables with respect to the coordinates. When necessary, these methods are also used to obtain variable values at locations other than grid nodes (interpolation).
On structured grids, the FD method is very simple and effective. It is especially easy to obtain higher-order schemes on regular grids. The disadvantage of FD methods is that the conservation is not enforced unless special care is taken. Also, the restriction to simple geometries is a significant disadvantage in complex flows.
Basic Steps in Finite Difference Method:
- In order to solve a given PDE by numerical methods, the partial differentials of the dependent variable in PDEs must be approximated by finite difference relations (algebraic equations).
- Solution of a steady PDE in a two-dimensional rectangular domain with initial and boundary conditions.
- Division of domain in uniform/non-uniform mesh.
- Formulation of algebraic equations for each grid point/control volume (matrix system Ax = b).
- Methods for finite difference approximations.
- Taylor series expansions.
- Finite Difference by Polynomials.
2. Finite Volume Method:
The FV method uses the integral form of the conservation equations as its starting point. The solution domain is subdivided into a finite number of contiguous control volumes (CVs), and the conservation equations are applied to each CV. At the centroid of each CV lies a computational node at which the variable values are to be calculated. Interpolation is used to express variable values at the CV surface in terms of the nodal (CV-center) values. Surface and volume integrals are approximated using suitable quadrature formulae. As a result, one obtains an algebraic equation for each CV, in which a number of neighbor nodal values appear.
The FV method can accommodate any type of grid, so it is suitable for complex geometries. The grid defines only the control volume boundaries and need not be related to a coordinate system. The method is conservative by construction, so long as surface integrals (which represent convective and diffusive fluxes) are the same for the CVs sharing the boundary.
The FV approach is perhaps the simplest to understand and to program. All terms that need be approximated have physical meaning which is why it is popular with engineers. The disadvantage of FV methods compared to FD schemes is that methods of order higher than second are more difficult to develop in 3D. This is due to the fact that the FV approach requires three levels of approximation: interpolation, differentiation, and integration.
Basic Steps in Finite Volume Method:
- Divide the continuous domain into a number of discrete subdomains (control volumes) by a grid. The grid defines the boundaries of a control volume, while the computational node lies at the center of each control volume.
- For each sub-domain, derive governing algebraic equations from the governing differential equations.
- Obtain a system of algebraic equations from above.
- Solve the above system of algebraic equations to obtain values of the dependent variables at identified discrete points (computational nodes).
3. Finite Element Method:
The FE method is similar to the FV method in many ways. The domain is broken into a set of discrete volumes or finite elements that are generally unstructured; in 2D, they are usually triangles or quadrilaterals, while in 3D, tetrahedra or hexahedra are most often used. The distinguishing feature of FE methods is that the equations are multiplied by a weight function before they are integrated over the entire domain. In the simplest FE methods, the solution is approximated by a linear shape function within each element in a way that guarantees continuity of the solution across element boundaries. Such a function can be constructed from its values at the corners of the elements. The weight function is usually of the same form.
This approximation is then substituted into the weighted integral of the conservation law and the equations to be solved are derived by requiring the derivative of the integral with respect to each nodal value to be zero; this corresponds to selecting the best solution within the set of allowed functions (the one with minimum residual). The result is a set of non-linear algebraic equations.
An important advantage of finite element methods is the ability to deal with arbitrary geometries; there is an extensive literature devoted to the construction of grids for finite element methods. The grids are easily refined; each element is simply subdivided. Finite element methods are relatively easy to analyze mathematically and can be shown to have optimality properties for certain types of equations. The principal drawback, which is shared by any method that uses unstructured grids, is that the matrices of the linearized equations are not as well structured as those for regular grids making it more difficult to find efficient solution methods.
Basic Steps in Finite Element Method:
The analysis of a structure by the Finite Element Method can be divided into several distinctive steps. These steps are to a large extent similar to the steps defined for the matrix method. Here we give a theoretical approach to the method, and its different steps.
- Discretization is the process of dividing your problem into several small elements, connected with nodes. All elements and nodes must be numbered so that we can set up a matrix of connectivity. The picture to the right shows discretization of a transverse frame into beam elements and discretization of a plane stress problem into quadrilateral elements. It is important to remember that the order the nodes and elements are numbered greatly affects the computing time. This is because we get a symmetrical, banded stiffness matrix, which bandwidth is dependent on the difference in the node numbers for each element, and this bandwidth is directly connected width the number of calculations the computer has to do. Computer FEM-programs have internal numbering that optimizes this bandwidth to a minimum by doing some internal renumbering of nodes if they are not optimal.
- The element analysis have two key components; Expressing the displacements within the elements, and maintaining equilibrium of the elements. In addition, stress-strain relationships are needed to maintain compatibilty. The final result is the element stiffness relationship: S = kv. For beam elements this relationship was obtained using the exact relationships between forces and moments and the corresponding displacements.These results could therefore be interpreted as being obtained by the governing differential equation and boundary condition of the beam elements. For e.g. a plane stress problem it is not possible to use an exact solution. The displacements within the elements are expressed in terms of shape functions scaled by the node displacements. Hence, by assuming expressions for the shape functions, the displacements in an arbitrary point within the element is determined by the nodal point displacement. The section of the structure that the element is representing is kept in place by the stresses along the edges. In the finite element analysis it is convenient to work with nodal point forces. The edge stresses may in the general case be replaced by equivalent nodal point forces by demanding the element to be in an integrated equilibrium using work or energy considerations. This technique is often referred to as to “lump” the edge forces to nodal forces.
This requirement result in a relationship between the nodal point displacements and forces to be given as:
S = kv + S0, where:
S– generalized nodal point forces
k– element stiffness matrix
v– nodal point displacements
S0– nodal point forces for external loads
Computer programs usually have many options for types of elements to choose among. Here the most usual elements:
- In System Analysis a relationship between the load and the nodal point displacements is established by demanding equilibrium for all nodal points in the structure:
The stiffness matrix is established by directly adding the contributions from the element stiffness matrices. Similarly the load vector R is obtained from the known nodal forces.
- Boundary conditions are introduced by setting nodal displacements to known values or spring stifnesses are added.
- The global displacements are found by solving the linear set of equations stated above:
r = K-1(R – R0)
- The stresses are determined from the strains by Hooke’s law. Strains are derived from the displacement functions within the element combined with Hooke’s law. They may be expressed generally by:
v = a r
D – Hooke’s law on matrix form
B – Derived from u(x, y, z)
Output interpretation programs, called post-processors, help the user sort out the output and display it in graphical form.
- Computational Methods for Fluid Dynamics by J.H. Ferziger and M. Peric.
- Lecture notes on Computational Fluid Dynamics by Dr. Tariq Talha, College of EME, NUST, Islamabad, Pakistan.
- Finite volume method basics, principles and examples by Aniket Batra, IIT Guwahati.