Finite volume Method(FVM) Literature Review

Aim:-

• Define Finite volume Method(FVM)
• Write down Difference Between FDM and FVM.
• Describe the need for interpolation schemes and also describe Flux limiters in FVM, use 1D linear heat conduction equation as an example.

Objective:-

• Define finite volume method (FVM). Also describe difference between FDM, FVM, and FEM.
• Consider example of 1D heat conduction equation to describe how to solve FVM problem.
• Explain about grid generation process in FVM
• Describe how to descretise the equation by using various interpolation method.
• Write down how flux limiter is important to reduce shocks and distortion is solution.
• Calculate the solutions for discretise equations.
• Lastly, compare calculated numerical results with actual results.

Finite Volume Method

• The Finite Volume method (FVM) is a discretization method for the approximation of a single or a system of partial differential equations expressing the conservation, or balance, of one or more quantities.
•  These partial differential equations (PDEs) are often called conservation laws. They describe the relationship between partial derivatives of unknown fields such as temperature, concentration, pressure, molar fraction, density of electrons or probability density function, with respect to variables within the domain (space, time,..) under consideration.
• Therefore, the finite volume method is used by various numerical techniques in which the solution domain is descretised into a number of control volumes over which integrals of conservation laws are directly discretised.

Difference between FDM, FVM and FEM

Finite difference method (FDM)

• Finite difference method is created from basic definition of differentiation that is `(df)/(dx)=(f(x+h)-f(x))/h`
• Easy method but not reliable for conservative differential equations and solutions having shocks.
• Though to implement in complex geometry where it needs complex mapping and mapping makes governing equation even tougher.
• Extending to higher order accuracy is very simple.
• In FDM discretization is based upon the differential form of the PDE to be solved. Each derivative is replaced with an approximate difference formula (that can generally be derived from a Taylor series expansion). The computational domain is usually divided into hexahral cells (the grid), and the solution will be obtained at each nodal point.
• The FDM is easiest to understand when the physical grid is Cartesian. We cannot implement this method, went the grids are not aligned along x, y, or z direction.
• The discretization results in a system of equation of the variables at nodal points, and once a solution is found, then we have a discrete representation of the solution.

Finite volume method (FVM)

• It is a numerical tool that is borrowed from calculus of variation. Here we assume some trial function and multiply that trial function with weighting function. Then this weighting function is multiplied with trial function then integrated over the control volume (weak from) and equated to zero(This procedure will differ for different types of FEM but theme is same). Then we get one set of algebraic equations. Solving that will give solution.
• Here we are working only in error and differential equation some times conservative law may be violated. This method is more accurate than PDEs. Here higher order accuracy is achieved by using higher order basis (i.e) shape functions. Extending to higher order accurate calculations are expensive in computation and Mathematical formulation especially for non-linear PDEs. Mostly suitable for Heat transfer, Structural mechanics, vibrational analysis etc.
• In finite volume method discritization is based upon an integral form of the PDE to be solved (e,g. conservation of mass, momentum, or energy). The PDE is written in a form which can be solved for a given finite volume (or cell). The computational domain is discretized into finite volume and then for every volume the governing equations are solved.
• The resulting system of equation usually involves fluxes of the conserved variable, and thus the calculation of fluxes is very important in FVM.
• The basic advantage of this method over FDM is it does not require the use of structured grid, and the effort to convert the given mesh into structured numerical grid internally is completely avoided. 
• As with FDM, the resulting approximate solution is a discrete, but the variables are typically placed at cell centers rather than at nodal points. This is not always true. As there are also face-centered finite volume methods. In any case, the values of field variables at non-storage locations (e.g. vertices) are obtained using interpolation.

Finite element method (FEM)

• This is similar to FDM, but, it didn’t kill the theme of differentiation because we are integrating the differential equation over a control volume and discretizing the domain.
• Since we have integrated the differential equation discetization is mathematically a valid one. It can be loosely viewed as FEM but weight here used is 1. Here fluxes are integrated and resultant is set to zero, so flux is conserved.
• It can handle almost any PDEs and complex domain. Interpolation is done from face to center will reduce the accuracy of this process. Here accuracy is based on order of polynomial used.
• In finite element method discretization is based upon a piecewise representation of the solution in terms of specified basis functions. The computational domain is divided up into smaller domains (finite elements) and the solution in each elements is constructed from the basis functions.
• The actual equation in weak form: the field variables are written in terms of the basis functions, the equation is multiplied b appropriate test functions, and then integrated over an element.
• A system of equations is obtained (usually for nodal values) that must be solved to obtain a solution.

To understand FVM more effectively consider following problem

Problem statement

As below figure shows a large plate of thickness L = 2 cm with constant thermal conductivity k = 0.5 W/m.K and uniform heat generation q = 1000 kW/m3 . The faces A and B are at temperatures of 100°C and 200°C respectively.

Governing equation for 1D linear heat conduction equation is

                              `d/(dx)(k(dT)/(dx))+q=0`                                                                                                                                      `…(1)`

Where,

k is the thermal conductivity

T is temperature

dx is the control volume

q is heat generation in the system

                                                                                                      

Write all the given condition in program.

%% Given input for problem
% plate thickness in mm
L = 0.02;
% Thermal conductivity in W/m.K
k = 0.5;
% Heat generation in kW/m^3
q = 1000e3;

% Initial temperature at point A
TA = 100;
% Initial temperature at point B
TB = 200;

%% Consider Area 
A = 1;

>> Grid generation

The fist step in the finite volume method is to divide the domain into descrete control volumes. Lest us place a number of nodal points in the space between `T_A` and `T_B`. The boundaries (or faces) of control volume are positioned mid-way between adjacent nodes. Thus each nodes is surrounded by a control volume or cell. It is common practice to set up control volumes near the edge of the domain in such a way that the physical boundaries coincide with the control volume boundaries.

                                          

In above figure, a nodal point is identified by p and its neighbours in a one-dimensional geometry, the nodes to the west and east, are identified by W and E. The west side face of the control volume is referred to by e and the east side control volume face by e. and length of the control volume is referred as `Delta x`.

Defining grid in program

%% Grid generation
% divided by 5 volumes
nx = 5;
dx = L/nx;
% specify center node of each volume
x = dx/2:dx:L-(dx/2);

>> Discretisation

The key step of the finite volume method is the integration of the governing equation (or equation) over a control volume to yield a discretised equation at its nodal point P. For the control volume defined above gives.

                              `int_(Delta V) d/(dx)(k (dT)/(dx))dV + int_qdV=0`                                                                                                      `…2`

                              `[(kA(dT)/(dx))_e-(kA(dT)/(dx))_w]+q Delta V = 0`                                                                                       `…3`

Here A is the cross-sectional are of the control volume face, `Delta V` is the volume

It is a very attractive feature of the finite volume method that the discretised equation has a clear physical interpretation.

Above equation states that the temperature leaving the east face minus the temperature entering the west face is equal to the generation of T, i.e. it constitutes a balance equation for T over the control volume.

In practical situations, the heat generating source may be a function of the dependent variable. In such cases the finite volume method approximates the heat generation term by means of a linear form:

                              `bar q DeltaV = qAdelx`  

In order to derive useful forms of the discretised equations, the interface thermal conductivity k and the gradient `(dT)/(dx)` at east(e) and west(w) are required. To calculate gradients temperature at the control volume faces an approximate distribution of properties between nodal points is used. Linear approximations seems to be the obvious and simplest way of calculating inner face value and the gradients. To do this we first need to understand interpolation scheme.

Finite Volume Interpolation schemes

In the field of numerical analysis, interpolation is a type of estimation, a method of constructing new data points within the range of a discrete set of known data points

The approximation of surface and volume integrals may require values of the variables at locations other than the computational nodes of the control volume I.e. to evaluate temperature values and temperature gradient we are going to use interpolation technique through the nodal value.

Now, these interpolation schemes are not so different from Taylor series expansion in the finite difference methods, nevertheless, in case of FVM, differentials need to be evaluate on the overall control volume rather than nodal points.

In equation (3) term `kA(dt)/(dx)` depicts the flux in the defined region from west to east of control volume and equation itself define how flux is being conserved inside the control volume. The heat is generated in control volume as the flux enter in the control volume form west and leave from east as shown in below figure.

                                                                                   

Following are the various interpolation scheme for different orders.

[1] First-order Upwind differencescheme (UD)

As its name suggest it is the first order accurate simplest numerical scheme, where value at the interface assume same as the value at the center of the control volume depending upon the flow direction.

This interpolation scheme is equivalent to using a backward or forward difference approximation (depending upon flow direction) which gives stable calculation but it is very diffusive.

                              `((del phi)/(del x))_P = (phi_E - phi_P)/(Delta x)`

Since it uses values at point E and P (where `x_E>x_P`) to evaluate the gradient `((del phi)/(del x))` at P, the formula is called forward differencing

Similarly for backward difference  `((del phi)/(del x))` at P will be

                              `((del phi)/(del x))_P = (phi_P - phi_W)/(Delta x)`

By using first order up wind scheme, we will get

                              `(kA(dT)/(dx))_e = k_eA((T_E - T_P)/(delx))`

                              `(kA(dT)/(dx))_w = k_wA((T_P-T_W)/(delx))`

After put it in the equation, we will get

                              `[k_eA((T_E - T_P)/(delx))-k_wA((T_P-T_W)/(delx)] + qAdelx = 0`                                                                       `…(4)`

[2] Central difference scheme (CD)

It is second order accurate approximation scheme in which value of the flux is defined by  linear interpolation amid the center of the control volume. 

                              `((del phi)/(del x))_P = (phi_E - phi_W)/(2 Delta x)`

Above equation uses values at E and W to evaluate the gradient at the mid-point P, and is called a central difference formula. The quadratic dependence of the error on grid spacing means that after grid refinement the error reduces more quickly in a second-order accurate different scheme than in a first-order accurate.

By using central difference scheme, we will get

                              `(kA(dT)/(dx))_e = k_eA((T_E - T_P)/(2delx))`

                              `(kA(dT)/(dx))_w = k_wA((T_P-T_W)/(2delx))`

After put it in the equation, we will get

                              `[k_eA((T_E - T_P)/(2delx))-k_wA((T_P-T_W)/(2delx))] + qAdelx = 0`                                                                  `…(5)`

[3] Second order upwind or linear upwind differencing (LUD) scheme

The spatial accuracy of the first-order scheme can be improves by including 3 data points instead of just 2, which offers a more accurate finite difference stencil for the approximation of spatial derivative

                              `((del phi)/(del x))_P = (3phi_E-phi_W)/(2 Delta x)`

By using linear upwind differencing scheme, we will get 

                              `(kA(dT)/(dx))_e = k_eA((3 T_E - T_P)/(2delx))`

                              `(kA(dT)/(dx))_w = k_wA((3T_P-T_W)/(2delx))`

This scheme is less diffusive compare to the first-order accurate scheme and is called linear upwind differencing (LUD) scheme.

After put it in the equation, we will get

                              `[k_eA((3T_E - T_P)/(2delx))-k_wA((3T_P-T_W)/(2delx))] + qAdelx = 0`                                                               `…(6)`

[4] Quadratic upwind interpolation (QUICK)

A quadratic upwind interpolation (QUICK) scheme can be derived using polynomial fitting. It is relatively straightforward to demonstrate the third-order accuracy of the QUICK differencing scheme for the convective flux at a midpoint cell face in a uniform grid of spacing Delta x. The QUICK scheme calculates the value phi_e at the east cell face of a general node as

                              `phi_e=3/8phi_E+6/8phi_P-1/8phi_W`

                              `(kA(dT)/(dx))_e = k_eA((3/8T_E+6/8T_P-1/8T_W)/(delx))`

                              `(kA(dT)/(dx))_w = k_wA((3/8T_E+6/8T_P-1/8T_W)/((delx)/2))`

After put it in the equation, we will get

                              `[k_eA((3/8T_E+6/8T_P-1/8T_W)/(delx))-k_wA((3/8T_E+6/8T_P-1/8T_W)/((delx)/2))] + qAdelx = 0`                   `…(7)`

 [5] Total variation diminishing (TVD) scheme

Schemes of third-order and above have been developed for the dicretisation of convective terms with varying degrees of success. Implementation of boundary conditions can be problematic with such higher-order schemes.  The fact that the QUICK scheme and other higher-order schemes can give

undershoots and overshoots has led to the development of second-order schemes that avoid these problems. The class of TVD (total variation diminishing) schemes has been specially formulated to achieve oscillation-free solutions and has proved to be useful in CFD calculations. TVD is a property used in the discretisation of equations governing time-dependent gas dynamics problems. More recently, schemes with this property have also become popular in general-purpose CFD solvers.

The general form of the east face value phi_e within a discretisation scheme my be written as

                              `phi_e = phi_p+1/2 psi( r)(phi_E-phi_P)`

The function of psi for various discretisation scheme

For the UD scheme `psi(r)` = 0

For the CD scheme `psi(r)` = 1

For the LUD scheme `psi(r)` = r

For the QUICK scheme `psi(r)` = `(3+r)/4`

                                                                    

Above figure shows the relationship between `psi(r)` VS r for above four scheme. This diagram is known as the `r-psi` diagram. All the above expression assume that the flow direction is positive (i.e. from west to east). it can be shown that similar expression exist for negative flow direction and r will still be the ratio of upwind-side gradient.

Flux Limiter:-

Before we understand flux limiter we first have to know criteria for TVD schemes

        Criteria for TVD scheme

There are necessary and sufficient conditions for a scheme to be TVD in terms of r-psi.

• If 0
• If r>=1 the upper limit is `psi(r )`=2, so for TVD schemes `psi(r )`<=2

Below figure shows the shaded TVD region in a `r-psi` diagram along with the `r-psi` relationships for all the finite difference schemes that we have discussed above.

                                                                     

• the UD scheme is TVD
• the LUD scheme is not TVD for r > 2.0
• the CD scheme is not TVD for r < 0.5
• the QUICK scheme is not TVD for r < 3/7 and r > 5

Except for UD, all the above schemes are outside the TVD region for certain values of r. The idea of designing a TVD scheme is to introduce a modification to the above schemes so as to force the `r-psi` relationship to remain within the shaded region for all possible values of r. This would imply that, in order to make the scheme TVD, we must constrain or limit the range of possible values of the additional convective flux `F_e psi(r)(phi_E-phi_P)/2`, which was originally introduce to make the scheme higher order. Hence, the function `phi(r)` is called a flux limiter function.

As we learn about TVD, flux limiters are used in high-resolution schemes to avoid spurious oscillation that would otherwise occur in high order discretization schemes due to shocks, discontinuities or sharp changes in the solution domain.

Use of flux limiters, together with an appropriate high resolution scheme, make the solutions total variation diminishing (TVD)

The flux limiter limit the spatial derivatives to physically realizable and meaningful values. It limits the solution gradient near shock or discontinuities.

>> Solution of equations

Discretised equations must be set up at each of the nodal points in order to solve a problem. For control volumes that are adjacent to the domain boundary the general discretised equation is modified to incorporate boundary conditions. The resulting system of linear algebraic equations is then solved to obtained the distribution of the temperature at nodal points.

We are going to solve this problem using First-order Upwind difference scheme (UD)

So equation(3) can be rearranged as

                              `((k_eA)/(delx)+(k_wA)/(delx))T_p - ((k_wA)/(delx))T_w-((k_e)/(delx))T_E=qAdelx`           

Thermal conductivity is constant all over so  `k_e=k_w=k`

                              `(2(kA)/(delx))T_p - ((kA)/(delx))T_w-((kA)/(delx))T_E=qAdelx`                

Equation(9) is valid for control volumes at nodal points 2, 3, and 4.

Now, lest define above equation in coding

%% Defining matix size according to nx
LHS = zeros(nx,nx);
RHS = zeros(nx,1);
 
%% apply First-order upwind scheme (points between A and B)
for i=2:nx-1
    LHS(i,i) = 2*(k*A)/dx;
    LHS(i,i-1) = -(k*A)/dx;
    LHS(i,i+1) = -(k*A)/dx;
    RHS(i) = q*A*dx;
end

To incorporate the boundary conditions at nodal 1 and 5 we apply the linear approximation for temperature between a boundary point and adjacent nodal point. At node 1 the temperature at the west boundary is known. Integration of equation at the control volume surrounding node 1 give

                              `[(kA(dT)/(dx))_e-(kA(dT)/(dx))_w]+q Delta V = 0`

Introduction of the linear approximation for temperatures between A and P yields

                              `[k_eA((T_E - T_P)/(delx))-k_A A((T_P-T_A)/((delx)/2))] + qAdelx = 0`              

The above equation can be rearranged, using `k_e=k_A=k`, to yield the discretised equation for boundary 1 :

                              `(3(kA)/(delx))T_P - ((kA)/(delx))T_E=(2(kA)/(delx))T_A+qAdelx`                    

% At A point
LHS(1,1) = 3*(k*A)/dx;
LHS(1,2) = -(k*A)/dx;
RHS(1) = 2*(k*A*TA)/dx + q*A*dx;

At the nodal point 5, the temperature on the east face of the control volume is known. The node is treated in a similar way to boundary node 1. At boundary point 5 we have

                              `[(kA(dT)/(dx))_e-(kA(dT)/(dx))_w]+q Delta V = 0`

                              `[k_BA((T_B - T_P)/((delx)/2))-k_w A((T_P-T_W)/(delx))] + qAdelx = 0`   

The above equation can be rearranged, noting that k_B = k_w = k, to give the discretised equation for boundary node 5 :

                              `(3(kA)/(delx))T_P - ((kA)/(delx))T_W=(2(kA)/(delx))T_B+qAdelx`              

% At B point
LHS(end,end) = 3*(k*A)/dx;
LHS(end,end-1) = -(k*A)/dx;
RHS(end) = 2*(k*A*TB)/dx + q*A*dx;

By putting above value in codding, we will get matrix for RHS and LHS and we can arrange matrix like this `[RHS][T]=[LHS]`

By solving above matrices we can get all the value for T

%% Rearranging matrices to calculate temperature matrix 
T = LHS^-1*RHS

Analytical solution for comparison

The analytical solution to this problem may be obtained by integrating equation (1) twice with respect to x and by subsequent application of the boundary conditions. This gives

                              `T=[(T_B-T_A)/L+q/(2k)(L-x)]x+T_A`                                                          

%% Analytical solution 
x_an = linspace(0,L,100); % Break L into 100 parts
T_an = ((TB-TA)/L+(q/(2*k))*(L-x_an)).*x_an+TA

Using above calculated value of T by analytical and numerical method we can plot its graph and compare the results

%% Plotting
plot(x_an*100,T_an,'linewidth',1.5) %x*100 to convert unit from mm to cm
hold on
plot(x*100,T,':o','markerfacecolor','g','color','r')
grid on 
axis([0 2 90 280])
title({'Comparison of the Numerical result';'with the Analytical solution'},'fontsize',14)
xlabel('Length in cm','fontweight','bold')
ylabel('Temperature in ^oC','fontweight','bold')
legend('Analytical solution','Numerical solution','location','southeast')

Output for above codes

                                                           

Results

Above figure exhibits the temperature value calculated for plate which have thickness 2 cm with both analytical and numerical method. We can see, even with a coarse grid of five nodes, the agreement is very good. We can, obviously, increase the accuracy by increasing the number of grid within length L.

Now, lets found out the maximum error we are getting from our numerical solution.

Error

To find out error just find out the difference between temperature calculated by analytical method and temperature calculated by numerical method.

%% Error
T_an_pt = ((TB-TA)/L+(q/(2*k))*(L-x)).*x+TA;
error = max(abs(T'-T_an_pt));

                                                                                           

So, error for our numerical solution for 5 control volume is 4. we can reduce it by making grid fine.

Reference

Book - An Introduction to Computational Fluid Dynamics - THE FINITE VOLUME METHOD -  H K Versteeg and W Malalasekera

Complete Program

clear all
close all
clc


%% Given input for problem
% plate thickness in mm
L = 0.02;
% Thermal conductivity in W/m.K
k = 0.5;
% Heat generation in kW/m^3
q = 1000e3;

% Initial temperature at point A
TA = 100;
% Initial temperature at point B
TB = 200;

%% Consider Area 
A = 1;


%% Grid generation
% divided by 5 volumes
nx = 5;
dx = L/nx;
% specify center node of each volume
x = dx/2:dx:L-(dx/2);



%% Defining matix size according to nx
LHS = zeros(nx,nx);
RHS = zeros(nx,1);


%% apply First-order upwind scheme (points between A and B)
for i=2:nx-1
    LHS(i,i) = 2*(k*A)/dx;
    LHS(i,i-1) = -(k*A)/dx;
    LHS(i,i+1) = -(k*A)/dx;
    RHS(i) = q*A*dx;
end

% At A point
LHS(1,1) = 3*(k*A)/dx;
LHS(1,2) = -(k*A)/dx;
RHS(1) = 2*(k*A*TA)/dx + q*A*dx;

% At B point
LHS(end,end) = 3*(k*A)/dx;
LHS(end,end-1) = -(k*A)/dx;
RHS(end) = 2*(k*A*TB)/dx + q*A*dx;

%% Rearranging matrices to calculate temperature matrix 
T = LHS^-1*RHS


%% Analytical solution 
x_an = linspace(0,L,100); % Break L into 100 parts
T_an = ((TB-TA)/L+(q/(2*k))*(L-x_an)).*x_an+TA

%% Plotting
plot(x_an*100,T_an,'linewidth',1.5) %x*100 to convert unit from mm to cm
hold on
plot(x*100,T,':o','markerfacecolor','g','color','r')
grid on 
axis([0 2 90 280])
title({'Comparison of the Numerical result';'with the Analytical solution'},'fontsize',14)
xlabel('Length in cm','fontweight','bold')
ylabel('Temperature in ^oC','fontweight','bold')
legend('Analytical solution','Numerical solution','location','southeast')

%% Error
T_an_pt = ((TB-TA)/L+(q/(2*k))*(L-x)).*x+TA;
error = max(abs(T'-T_an_pt));