All Courses
All Courses
Courses by Software
Courses by Semester
Courses by Domain
Tool-focused Courses
Machine learning
POPULAR COURSES
Success Stories
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…
Syed Saquib
updated on 10 May 2023
Aim:-
Objective:-
Finite Volume Method
Difference between FDM, FVM and FEM
Finite difference method (FDM)
Finite volume method (FVM)
Finite element method (FEM)
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
ddx(kdTdx)+q=0���(�����)+�=0 …(1)…(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 TA�� and TB��. 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 Δ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.
∫ΔVddx(kdTdx)dV+∫qdV=0∫Δ����(�����)��+∫���=0 …2…2
[(kAdTdx)e−(kAdTdx)w]+qΔV=0[(������)�-(������)�]+�Δ�=0 …3…3
Here A is the cross-sectional are of the control volume face, Δ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:
¯qΔV=qA∂x�¯Δ�=��∂�
In order to derive useful forms of the discretised equations, the interface thermal conductivity k and the gradient dTdx���� 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 kAdtdx������ 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.
(∂ϕ∂x)P=ϕE−ϕPΔx(∂�∂�)�=��-��Δ�
Since it uses values at point E and P (where xE>xP��>��) to evaluate the gradient (∂ϕ∂x)(∂�∂�) at P, the formula is called forward differencing
Similarly for backward difference (∂ϕ∂x)(∂�∂�) at P will be
(∂ϕ∂x)P=ϕP−ϕWΔx(∂�∂�)�=��-��Δ�
By using first order up wind scheme, we will get
(kAdTdx)e=keA(TE−TP∂x)(������)�=���(��-��∂�)
(kAdTdx)w=kwA(TP−TW∂x)(������)�=���(��-��∂�)
After put it in the equation, we will get
[keA(TE−TP∂x)−kwA(TP−TW∂x]+qA∂x=0[���(��-��∂�)-���(��-��∂�]+��∂�=0 …(4)…(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.
(∂ϕ∂x)P=ϕE−ϕW2Δx(∂�∂�)�=��-��2Δ�
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
(kAdTdx)e=keA(TE−TP2∂x)(������)�=���(��-��2∂�)
(kAdTdx)w=kwA(TP−TW2∂x)(������)�=���(��-��2∂�)
After put it in the equation, we will get
[keA(TE−TP2∂x)−kwA(TP−TW2∂x)]+qA∂x=0[���(��-��2∂�)-���(��-��2∂�)]+��∂�=0 …(5)…(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
(∂ϕ∂x)P=3ϕE−ϕW2Δx(∂�∂�)�=3��-��2Δ�
By using linear upwind differencing scheme, we will get
(kAdTdx)e=keA(3TE−TP2∂x)(������)�=���(3��-��2∂�)
(kAdTdx)w=kwA(3TP−TW2∂x)(������)�=���(3��-��2∂�)
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
[keA(3TE−TP2∂x)−kwA(3TP−TW2∂x)]+qA∂x=0[���(3��-��2∂�)-���(3��-��2∂�)]+��∂�=0 …(6)…(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
ϕe=38ϕE+68ϕP−18ϕW��=38��+68��-18��
(kAdTdx)e=keA(38TE+68TP−18TW∂x)(������)�=���(38��+68��-18��∂�)
(kAdTdx)w=kwA(38TE+68TP−18TW∂x2)(������)�=���(38��+68��-18��∂�2)
After put it in the equation, we will get
[keA(38TE+68TP−18TW∂x)−kwA(38TE+68TP−18TW∂x2)]+qA∂x=0[���(38��+68��-18��∂�)-���(38��+68��-18��∂�2)]+��∂�=0 …(7)…(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
ϕe=ϕp+12ψ(r)(ϕE−ϕP)��=��+12�(�)(��-��)
The function of psi for various discretisation scheme
For the UD scheme ψ(r)�(�) = 0
For the CD scheme ψ(r)�(�) = 1
For the LUD scheme ψ(r)�(�) = r
For the QUICK scheme ψ(r)�(�) = 3+r43+�4
Above figure shows the relationship between ψ(r)�(�) VS r for above four scheme. This diagram is known as the r−ψ�-� 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.
Below figure shows the shaded TVD region in a r−ψ�-� diagram along with the r−ψ�-� relationships for all the finite difference schemes that we have discussed above.
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−ψ�-� 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 Feψ(r)ϕE−ϕP2���(�)��-��2, which was originally introduce to make the scheme higher order. Hence, the function ϕ(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
(keA∂x+kwA∂x)Tp−(kwA∂x)Tw−(ke∂x)TE=qA∂x(���∂�+���∂�)��-(���∂�)��-(��∂�)��=��∂�
Thermal conductivity is constant all over so ke=kw=k��=��=�
(2kA∂x)Tp−(kA∂x)Tw−(kA∂x)TE=qA∂x(2��∂�)��-(��∂�)��-(��∂�)��=��∂�
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
[(kAdTdx)e−(kAdTdx)w]+qΔV=0[(������)�-(������)�]+�Δ�=0
Introduction of the linear approximation for temperatures between A and P yields
[keA(TE−TP∂x)−kAA(TP−TA∂x2)]+qA∂x=0[���(��-��∂�)-���(��-��∂�2)]+��∂�=0
The above equation can be rearranged, using ke=kA=k��=��=�, to yield the discretised equation for boundary 1 :
(3kA∂x)TP−(kA∂x)TE=(2kA∂x)TA+qA∂x(3��∂�)��-(��∂�)��=(2��∂�)��+��∂�
% 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
[(kAdTdx)e−(kAdTdx)w]+qΔV=0[(������)�-(������)�]+�Δ�=0
[kBA(TB−TP∂x2)−kwA(TP−TW∂x)]+qA∂x=0[���(��-��∂�2)-���(��-��∂�)]+��∂�=0
The above equation can be rearranged, noting that k_B = k_w = k, to give the discretised equation for boundary node 5 :
(3kA∂x)TP−(kA∂x)TW=(2kA∂x)TB+qA∂x(3��∂�)��-(��∂�)��=(2��∂�)��+��∂�
% 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=[TB−TAL+q2k(L−x)]x+TA�=[��-���+�2�(�-�)]�+��
%% 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));
Leave a comment
Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.
Other comments...
Week - 4 - 2D meshing for Plastic components
14 Feb 2024 04:24 PM IST
Week 3 - 2D meshing for Sheet metal
14 Feb 2024 04:10 PM IST
Project
AIM: To carry out a system-level simulation of an All-Terrain Vehicle (ATV). OBJECTIVES : To carry out a Simulation of ATV. To prepare a technical report explaining the model properties & comments on the results. THEORY : All-Terrain Vehicle (ATV) An All-Terrain Vehicle (ATV), also known as a light utility…
03 Jan 2024 10:45 AM IST
Project 1
Aim : Develop a double-acting actuator model using Simscape Multibody and Simscape components. Objective : The mechanical system of the cylinder needs to be built using Simscape Multibody library components/blocks, and the hydraulic system needs to be modeled using Simscape library physical components. Theory : The…
16 Oct 2023 03:59 PM IST
Related Courses
Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.
© 2025 Skill-Lync Inc. All Rights Reserved.