Week 9 - FVM Literature Review

AIM

To describe the need for interpolation schemes and flux limiters in Finite Volume Method (FVM).

 

OBJECTIVE

To study and understand

1. What is Finite Volume Method(FVM)
2. Write down the major differences between FDM & FVM
3. Describe the need for interpolation schemes and flux limiters in FVM

 

INTRODUCTION

             We know that in Fluid Dynamics there are three numerical-method based solvers which are used to solve for the governing equation of Navier-Stokes: 

• Finite Difference Method
• Finite Element Method
• Finite Volume Method

 

THEORY

 

Finite Difference Method(FDM)

The finite difference method (FDM) is an approximate method for solving partial differential equations. Basically, it converts linear ordinary differential equations(ODEs) or non-linear partial differential equations(PDEs) into a system of equations that can be solved using matrix algebraic techniques. Here, the grid should be linear & must be aligned in such a way so that all the grid points lie in a straight line. Here, the governing equation's properties are said to be concentrated at each and every point along the grid. 

 

Finite Element Method(FEM)

            This method consists of the pictorial representation of the solutions in terms of defined basis functions. Here, the computational domain is partitioned into smaller domains of finite elements & the solution is built from the basis functions in each element. The actual equations are usually founded by restating the conservation equations in a weak form. Here, the field variables should be stated in terms of the basis functions, then the equation is multiplied by the appropriate test functions & the equation is then integrated over an element. Since the finite element method's solution is expressed in terms of specific basis functions hence far more information is available when compared to FDM or FVM. 

 

Finite Volume Method(FVM)

            FVM is a discretization method using which we can approximate a single PDE or even a system of PDEs that expresses the conservation of one or more quantities. These PDEs are often regarded as conservation laws that can be of various nature such as elliptic, parabolic, hyperbolic, etc. It is used to describe the relations between partial derivatives of unknown fields such as temperature, pressure, molar fraction, concentration, the density of electrons, etc with respect to the variables within the domain(time, space, etc) under consideration.

Here in FVM, a mesh is being generated that consists of a partition of the domain where the space variables reside. Elements of those meshes are termed as control volumes where the integration of those PDEs over each control volume results in a balanced equation. Then, those sets of balanced equations are then discretized with respect to the set of discrete unknowns. The issue here is that the discretization of those fluxes at the boundaries of each control volume as in order for the FVM to be efficient, the numerical fluxes are generally conservative and consistent. When the numerical fluxes are conservative then the flux entering a control volume from its neighbor should be the opposite of the one entering the neighbor from the control volume. When the numerical fluxes are consistent then the numerical flux of a regular function interpolation tends to the continuous flux as the mesh size vanishes. 

Sometimes it is possible to discretize the fluxes at the boundaries of the control volume by the FDM and then it is referred to as the conservative finite difference method. Specifically, in the case of FVM with respect to FDM, the discretization can be performed on the local balance equations, rather than on the PDEs, and hence the fluxes on the boundaries of the control volumes are discretized, rather than the continuous differential operator. 

 

Difference between FVM and FDM

SERIAL NO	FACTOR OF DIFFERENCE	FINATE DIFFERENCE METHOD (FDM)	FINATE VOLUME METHOD (FVM)
1	GRID GENERATION	In the case of FDM, the domain is divided into numerous points that are referred to as grid points each of which represents the respective variable of interest.	In the case of FVM, the domain is divided into numerous control volumes mostly referred to as cells which consist of variables of interest that are located at the centroid of each control volume or cell.
2	MODE OF DISCRETIZATION	In the case of FDM, the discretization is based on the derivatives since the differential form of governing equations is replaced by a system of algebraic equations at each grid point.	In the case of FDM, the discretization is based on the derivatives since the differential form of governing equations is replaced by a system of algebraic equations at each grid point.
3	APPLICABILITY	For FDM, it is almost mandatory to use structured grids as otherwise the conservation laws will be not be preserved.	The control volume can have grids in the shape of structured or unstructured grids. The conservation principles will be maintained in both cases regardless of the type of mesh. Conservation principles including mass, momentum, and energy are maintained by definition.
4	NUMERICAL SCHEMES	Numerical Schemes which are applicable in the case of FDM are Forward differencing, Central differencing & Backward differencing.	Numerical Schemes which are applicable in the case of FVM are Upwind differencing, Hybrid differencing, Power-law which is used for giving a more accurate approximation to the one-dimensional exact solution
5	ADVANTAGES	Implementation is easier, Higher-order schemes are available & it requires lower computational cost in general along with improved efficiency when compared to FVM while solving problems of similar nature as it uses simpler grid structures. If one can overcome the boundary-condition problem on curved boundaries, FDM gives very efficient and high-quality results.	FVM’s most significant advantage is that it only needs to do flux evaluation for the cell boundaries. Also, this is true for nonlinear issues as well, which, in turn, makes it an excellent choice for the handling of (nonlinear) conservation laws appearing in transport problems. The accuracy of the FVM, for example, close to a corner of interest, can be increased by refining the mesh around that corner.
6	DISADVANTAGES	We have to use the mapping functions in methods like the variable transformation technique where it is required when the points are not aligned along with the cartesian coordinates. The Finite difference method is more difficult to use for handling material discontinuities. It doesn't lend itself for local grid refinement or adaptive mesh refinement that is needed to resolve local rapid variations in solutions such as around a corner of a complex shape and hence it becomes increasingly challenging when dealing with more complex geometrical features.	In FVM, the functions that approximate the solution when using the finite volume method cannot be easily made of high-order as the presence of false diffusion exists since lower-order accuracy upwind schemes are preferred due to stability issues in higher-order accurate schemes like central differencing, and hence it is difficult to achieve higher-order schemes in comparison to FDM.
7	APPLICATIONS	We can deal with many physical problems using FDM such as predictions of quantities based on models developed to study weather along with its attributes, It has applications in the field of astrophysics and seismology as well. There aren't any known commercially available software using FDM although research codes are being developed to solve selective physical problems.	Using FVM, evaluation of external aerodynamics features of aircraft bodies, analysis of various kinds of heat transfer modes and mechanisms such as analysis of heat dissipation in case of a graphic processing unit(GPU) within an enclosed domain, assessment of wells in petroleum engineering, etc. COMMERCIAL SOFTWARES USING FVM: OpenFOAM, ANSYS CFX & FLUENT, PHOENICS, STAR-CD, STAR-CCM+

 

In FVM instead of writing the equation in a single point, equations are written in an entire volume. It means for a considered volume no assumptions are made over the shape of the volume. Inside the particular volume the flow is conserved, flow quantities are conserved, energy is conserved, mass is conserved, momentum is conserved and so on. This is called as 'Finite Volume Representation' which begins with 'Integral form of governing equation'.

Steps:-

1. once we have this integral form of equation.
2. Apply Gauss divergence theorem
3. We get descretized form which can be solved using computer ("Flux Terms").

consider a steady state heat conduction equation with heat generation:

`∂/(∂x)(α.(∂T)/(∂x))+S=0`

 

S is the source term

α is thermal diffusivity (function of temperature)

`(∂T)/(∂x)` is the temperature gradient along x direction

 

Considering finite volumes

Governing equation – `∂/(∂x) (α.(∂T)/(∂x))+S=0` is valid at every point in the domain, which means it is also valid in the entire control volume

If a small control volume (dv) is considered and multiplied with the governing equation and integrated over the entire volume, it gives the integral form of the governing equation

 

`∫(∂/(∂x) (α.(∂T)/(∂x))+S)dv=0`

`∫(∂/(∂x) (α.(∂T)/(∂x)))dv+∫Sdv=0`

`∵dv=A.dx`

`∫(∂/(∂x) (α.(∂T)/(∂x)))A.dx+∫Sdv=0`

`[αA((∂T)/(∂x))]_w^e+barS dv=0`

`bar S` is the average value of the source term (S), which when multiplied by dv is equal to`∫Sdv`

`α_e A (∂T)/(∂x_e )-α_w A (∂T)/(∂x_w )+barS dv=0`

 

Assuming α as a constant

`αA[(T_e-T_p)/(Δx)]-αA[(T_p-T_w)/(Δx)]+barS dv=0`

 

The terms inside the square braces are the heat flux terms. The source term needs to be linearized, if it is a function of temperature

`(barS underline("linearization" )→S_u+S_p T_p )`

`αA[(T_e-T_p)/(Δx)]-αA[(T_p-T_w)/(Δx)]+barS dv=0`

 

`[Heatflux]_"out" -[Heatflux]_"in" +Sourceterm=0`  The heat flow in a particular volume is conserved.

 

Different schemes are used to solve the equation

1. Gradient scheme
2. Interpolation schemes

 

1. Gradient scheme

When we calculate the temperature gradient in the east face, we use temperature from the neighbouring points. This is known as the gradient scheme.

 

2. Interpolation schemes

In FVM formulation, evaluation of surface and volume integrals may sometimes require the values of unknown variables at locations other than the computational nodes of the control volume. These values are obtained by employing certain interpolation schemes.

 

As an example, a simple 1D convection and diffusion is described.

`d/dx (ρuϕ)=d/dx (Γ (dϕ)/dx)`

 

where ϕ is the solution variable (scalar), ρ is the density, u is the velocity and Γ, the diffusion coefficient.

 

The above figure shows the three neighbouring control volumes, used for discretization.

 

Integration of the governing differential equation over the control volume around grid point P yields

`(ρuϕ)_e-(ρuϕ)_w=(Γ (dϕ)/dx)_e-(Γ (dϕ)/dx)_w`

 

The evaluation of the LHS of the above equation requires the value of ϕ over the faces

E and w. These values are estimated by employing certain interpolation schemes.

The approximations of surface and volume integrals require values of the variable at location other than the computational nodes of the control volume. Values at these locations are obtained using interpolation formulae.

Interpolation methods:

            Upwind interpolation

            Linear interpolation

            Quadratic upwind interpolation

            Hybrid, TVD and ENO interpolation schemes

 

1. Upwind interpolation (UDS)

In this scheme, values from the upwind/upstream of the cell centers should be defined. It is used in convection dominated problems. The upwind interpolation for approximating the value of a variable f at the east face of a control volume is given by

`f_e={(f_P,if (v,n)_e>0) , (f_E,if (v,n)_e<0):}`

 

This interpolation scheme is equivalent to the forward differencing method or backward differencing method depending on the flow direction.

 the flow is along the -ve x-direction, then the values at the interface are given by

for v<0;

`u_(i-1/2)≈u_(i-1);u_(i-1/2)≈u_(i-1)`

`I_c=v (u_i-u_(i-1))/(Δx)` which is equivalent to Backward differencing scheme

 

If the flow is along the +ve x-direction, then the values at the interface are given by

for v>0;

`u_(i-1/2)≈u_(i-1); u_(i+1/2)≈u_i`

`I_c=v (u_(i-1)-u_i)/(Δx)`  which is equivalent to Forward differencing scheme

 

It is said to be equivalent because all three methods have an order of accuracy 1.

Advantage:

• UDS is the only approximation that satisfies the boundedness criterion unconditionally irrespective of velocity and grid spacing i.e. never yield oscillatory solutions.

Disadvantages:

• UDS is numerically diffusive
• Numerical diffusion is magnified in multi-dimensional problems if the flow is oblique to the grid
• Peaks or rapid variations in the variables will be smeared out
• Only first-order accurate, hence very fine grids are required to obtain accurate solutions.

 

 

2. Linear interpolation (CDS)

 In this scheme, approximate values of the variable at the face centre of the control volume are obtained by the linear interpolation of the values from the two nearest computational nodes (one upstream and one downstream).

`f_e=f_E λ_e+f_P (1-λ_e )`

 

Where `λ_e=(x_e-x_p)/(x_E-x_P )`

It is more simple and widely used.

The linear interpolation is equivalent to the use of the central differencing formula of the first-order derivative and hence this scheme is also termed as Central Differencing Scheme (CDS)

Average interface values are given by

`u_(i+1/2)=(u_i+u_(i+1))/2`

`u_(i-1/2)=(u_(i-1)+u_i)/2`

`I_c=v (u_(i+1)-u_(i-1))/(2Δx)`

which is equivalent to CDS.

It is said to be equivalent because both these methods have an order of accuracy 2.

Advantages:

• Approximation of gradients has second-order of accuracy on uniform grids
• It's formally of first-order on non-uniform grids

Disadvantage:

• It may produce oscillatory solutions in case of large grid space

 

 

3. Hybrid interpolation scheme

It is a blend of UDS and CDS. The hybrid differencing scheme exploits the favourable properties of the upwind and central differencing schemes. It switches to upwind differencing when central differencing produces inaccurate results at high Pe numbers.

The scheme is fully conservative and since the coefficients are always positive it is unconditionally bounded. It satisfies the transportive requirement by using an upwind formulation for large values of Peclet number.

`f_e=γf_(CDS)+(1-γ) f_(UDS)`

Advantage:

• Produces physically realistic solutions and highly stable compared to QUICK

Disadvantage:

• Accuracy is of first-order

 

4. Quadratic Upwind Interpolation for Convective Kinematics (QUICK)

It approximates the values of the variable at the face center of a CV by quadratic interpolation of the values at three nearest computational nodes (1 downstream node D & 2 upstream nodes U, UU)

Quadratic - As it uses 3 nodes

Upwind - As it uses two upstreams

Convective Kinematics - As it is used for convective domain problems

for `u_w>0` ; `ϕ_w` is evaluated based on the values at WW, W, and P

for `u_e>0` ; `ϕ_e` is evaluated based on the values at W, P, and E

for `u_w<0`; `ϕ_w` is evaluated based on the values at W, P, and E

for `u_e<0` ; `ϕ_e` is evaluated based on the values at P, E, and EE

 

`f_e=f_E+g_1 (f_D-f_U )+g_2 (f_(UU)-f_U )`

 

where,

`g_1=((x_e-x_U)(x_e-x_(UU)))/((x_(UU)-x_U)(x_(UU)-x_D))`

`g_2=((x_e-x_U )(x_e-x_D ))/((x_(UU)-x_U )(x_(UU)-x_D ))`

 

Advantages:

• It is the better version of UDS and CDS because it is more accurate than UDS and more damping in nature than CDS as it uses 3 nodes to solve
• It has an order of accuracy 3

Disadvantages:

• The use of QUICK can lead to subtle problems caused by unbounded results
• Implementation of boundary conditions could be problematic with such higher-order schemes

 

5. Total Variation Diminishing scheme (TVD)

The fact that the QUICK scheme and other higher-order schemes can give overshoots and undershoots has led to the development of second-order schemes that avoid these problems. The class of Total Variation Diminishing (TVD) scheme has been specifically formulated to achieve oscillation-free solutions and has proved to be useful in CFD calculations.

TVD scheme is designed to address the undesirable oscillatory behavior of higher-order schemes. In the TVD scheme, the tendency towards oscillation is counteracted by adding an artificial diffusion fragment or by adding a weighting parameter towards upstream contribution. TCD has an order of accuracy 2.

`f_e=f_P+1/2 ψ(r)(f_E-f_P)`

where,

ψ(r) is called as flux limiter function

where,

`r=(u_i-u_(i-1))/(u_(i+1)-u_i )`

 

FLUX LIMITER:

For convective fluid flow, it is observed that the lower-order schemes are usually stable but quite dissipative in nature around the points of discontinuity or shocks while the higher-order schemes are unstable in nature and show oscillations in the vicinity of discontinuities or shocks. Highly accurate and oscillation-free schemes are knowns as high resolution schemes.

They operate functionally as a switch and tune the numerical fluxes obtained from lower order schemes and higher order schemes such that a lower order (stable) scheme is employed near the discontinuity and a higher order (accurate) scheme away from it.

In order to make the scheme TVD, we must constrain or limit the range of possible values of the additional convective flux `f_e=f_P+1/2 ψ(r)(f_E-f_P)` 

which was originally introduced to make the scheme higher-order. Hence, the function ψ(r) is called a flux limiter function.

Flux limiter functions are used to fine tune the higher order and lower order schemes in such a way that the resulting scheme gives a higher order accuracy in the smooth region of the flow and maintains first order accuracy in the vicinity of shocks and discontinuities.

Consider 1D semi discrete scheme

`(du_i)/dt+1/(Δx_i ) [F(u_(i+1/2) )-F(u_(i-1/2) )]=0`

`F(u_(i+1/2) ) & F(u_(i-1/2) )`  are edge fluxes for the "i"th cell

`F(u_(i+1/2) )=f_(i+1/2)^(low)-ϕ(r_i)(f_(i+1/2)^(low)-f_(i+1/2)^(high) )`

`F(u_(i-1/2) )=f_(i-1/2)^(low)-ϕ(r_i)(f_(i-1/2)^(low)-f_(i-1/2)^(high) )`

 

f^high - low precision, high resolution flux

f^low - high precision, low resolution flux

ϕ(r) - flux limiter function

`r_i=(u_i-u_(i-1))/(u_(i+1)-u_i )`

 

r represents the ratio of successive gradients on the solution mesh.

`ϕ(r)≥0`

This is the limiter function case

If the limiter function = 0, there will be a low resolution and sharp gradient.

If the limiter function > 0, there will be a high resolution and positive gradient.

 

NEED FOR FLUX LIMITERS IN FVM:

•  Interpolation schemes give good results.
• In case of shock flow problem its discontinuities, for that we need flux limiter.
• The limiter function is constrained to be greater than pr equal to zero.
• Therefore when the limiter is equal to zero the flux is represented by a low resolution scheme.
• similarly when the limiter is equal to  1 it is represented by a high resolution scheme. 

 

CONCLUSION

The Finite Volume method (FVM) uses a cell cantered approach to account for flow quantities. And there is need to discretize or more precisely interpolate the flow quantities with some dedicated interpolation scheme. For every physical characteristic, there is necessity to satiate the dependencies and flow characteristics within the domain in an appropriate manner. Most fundamental interpolation schemes were studied in detailed approach. There are some advanced schemes used in CFD codes which are combination of linear, upwind and linear upwind / hybrid differencing schemes.