Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective :
In this challenge, we will simulate an axi-symmetric flow through a cylinderical pipe.Simulating flow through a pipe will be computationally costly ,so as the flow through a pipe is axi-symmetric in nature we will split the domain and simulate only for the sector of the domain considering a wedge portion .Here we are simulating the Hagen poiseuille's flow which is an incompressible laminar flow occuring through a through a constant cross-section smooth cylinderical pipe and we consider an inlet Reynold's number of as 2100.We will do transient state simulation and validate our results with the Hagen- Poiseuille's equations for the fully developed flow.Here we will write a code in Matlab to automate the generation of blockMeshDict file and use it in OpenFoam. The following are the results we will show in this project and validate with analytical solution i.e, by Hagen poiseuille's flow equations,
- We will validate Hydro-dynamic length from Hagen poiseuille's flow with the numerical result.
- Then validate the fully developed flow velocity profile with its analytical profile by Hagen poiseuille's flow equations.
- We then validate maximum velocity and pressured drop for fully developed flow.
- Post-process Shear stress and validate wall shear stress for fully developed flow.
Theory :
The Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, is a physics law that gives the pressure drop in an incompressible and Newtonian fluid inlaminar flow flowing through a long cylindrical pipe of constant cross section.
Consider a fluid entering a circular pipe at a uniform velocity. Because of the no-slip condition, the fluid particles in the layer in contact with the surface of the pipe come to a complete stop. This layer also causes the fluidparticles in the adjacent layers to slow down gradually as a result of friction.To make up for this velocity reduction, the velocity of the fluid at the mid-section of the pipe has to increase to keep the mass flow rate through thepipe constant. As a result, a velocity gradient develops along the pipe.The region of the flow in which the effects of the viscous shearing forces caused by fluid viscosity are felt is called the boundary layer.
The hypothetical boundary surface divides the flow in a pipe into two regions: the boundary layer region, in which the viscous effects and the velocity changes are significant, and the irrotational (core) flow region, in which the frictional effects are negligible and the velocity remains essentially constant in the radial direction. The thickness of this boundary layer increases in the flow direction until the boundary layer reaches the pipe center and thus fills the entire pipe, as shown in below figure. The region from the pipe inlet to the point at which the boundary layer merges at the centerline is called the hydrodynamic entrance region, and the length of this region is called the hydrodynamic entry length L_h.

Flow in the entrance region is called hydrodynamically developing flow, since this is the region where the velocity profile develops. The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged is called the hydrodynamically fully developed region.
Operating conditions and Analytical solution :
Reynolds Number :The transition from laminar to turbulent flow depends on the geometry, surface roughness, flow velocity, surface temperature, and type of fluid, among other things.The flow regime depends mainly on the ratio of inertial forces to viscous forces in the fluid. This ratio is called the Reynolds number and is expressed for internal flow in a circular pipe as

where Vavg = average flow velocity (m/s),rho =density, D = characteristic length of thegeometry (diameter in this case, in m), and nu = mu/rho = kinematic viscosity of the fluid (m2/s). Note that the Reynolds number is a dimensionless quantity. Also, kinematic viscosity has the unit m2/s, and can beviewed as viscous diffusivity or diffusivity for momentum.
Hydrodynamic entry length :
since the flow is laminar ,The hydrodynamic entry length is given as, L_h = 0.05*Re*D.

Velocity profile and Average velocity :
The velocity profile in fully developed laminar flow in a pipe is parabolic with a maximum at the centerline and minimum (zero) at the pipe wall.

This is a convenient form for the velocity profile since V_avg can be determined easily from the flow rate information. The maximum velocity occurs at the centerline and is determined from
above equation by substituting r = 0,
V_max =2*V_avg
Average velocity,V_avg is defined as the average speed through a cross section. For fully developed laminar pipe flow, V_avg is half of maximum velocity.In fluid flow, it is convenient to work with an average velocity V_avg, which remains constant in incompressible flow when the cross-sectional area of the pipe is constant.Also from the definition of renold's number we can get the average velocity.
Pressure drop(ΔP) : The Pressure gradient shows the difference in the pressure between the two ends of the tube , determined by the fact that any fluid will always flow from high pressure to low pressure region and the flow rate is determined by the pressure gradient.Here as the Hagen poiseuille's equation is valid only for fully developed flow we will consider the pressure difference between the hydrodynamically fully developed region.

R= Radius of tube .The liquid flow varies directly with the radius.
Viscosity( μ) : The flow of the fluid varies inversely with the viscosity of the fluid and as the viscosity of the fluid increases, the flow decreases vice versa.
Length of the Tube (L): The liquid flow is inversely proportional to the length of the tube, therefore longer the tube, greater is the resistance to the flow.
Solution procedure :
- First we start by calculating the Analytical Solution for flow variables by Hagen–Poiseuille equations from the above mentioned equations .These are useful to validate our numerical solution and also we give initial conditions based on these values .
- Then comming to Numerical solution we start with writing a code in Matlab to automate the generation of blockMeshDict file to use it in OpenFoam.
- Then as the flow is laminar and incompreesible we use the solver IcoFoam which solves laminar incompressible navier stokes equations.Then we select a tutorial under IcoFoam i.e,cavity then we replace the existing BlockMeshDict file which was present , with the BlockMeshDict file which we created in matlab.
- Then we give the initial conditions for pressure and velocity by editing the p and U file in 0 folder.
- As the working fluid is water we change the kinematic viscosity which is already present with kinematic viscosity of water which is present in transport properties under constant folder.
- Then we change the time step and time in the controlDict file such that the simulation is stable and perfectly executed.
- Then we run the solver IcoFoam and and plot the results and validate them with analytical solution i.e, by Hagen poiseuille's flow equations.
Analytical Solution for flow variables by Hagen–Poiseuille equation :
Hence to find the analytical solution we need to have the properties of working fluid i.e, here water ,and the pipe dimensions and flow conditions.Therefore lets consider ,
- Diameter of pipe, D =0.02m
- radius of pipe, r= D/2= 0.01m
- Reynolds's number , Re =2100
- Wedge angle, theta = 4(should be less than 5 degrees)
- Density of water, rho = 997 kg/m^3
- Dynamic viscosity of water ,mu = 8.9*1e-4 Pa.s
- kinematic viscosity, nu = mu/rho= 8.9268e-07 m^2/s
Now from the Hagen–Poiseuille flow conditions we get,
- Hydrodynamic entry length,L_he = 0.05*Re*D
L_he =0.05*2100*0.02 =2.1m
- Length of pipe, L_p =L_he+1
L_p =2.1+1 =3.1m
- Average velocity, V_avg = Re*nu/D
V_avg = 2100*8.9268e-07/0.02 = 0.0937 m/sec
- Maximum velocity ,V_max =2*V_avg
V_max = 2*0.0937 = 0.1875 m/sec
- Pressure drop,ΔP = (8*V_avg*mu*L_p)/(r^2)
del_P(ΔP) = (8*0.0937*8.9*1e-4*1)/(0.01)^2 = 6.6737 Pa
(Here as the Hagen poiseuille's equation is valid only for fully developed flow we will consider the pressure difference between the hydrodynamically fully developed region. i.e,L_p =1m)
- Since here OpenFoam deals with kinematic pressure ,k_P =del_P/rho
k_P = 6.6737/997 = 0.0067 m^2/sec^2
Numerical Solution :
We start with writing a code in Matlab which automatically creates the vertices by the wedge angle that we specify.Here the geometry we create is wedge which is made of single block.This code then automate the generation of blockMeshDict file to use it in the OpenFoam.This code also gives the analytical solution and initial conditions for the simulation by taking above initial values and operating conditions by solving Hagen poiseuille's flow conditions.
Matlab code :
clear all
close all
clc
Re =2100 ;% reynolds's number
D =0.02 ; %diameter of pipe in metres
r= D/2; % radius of pipe
theta = 4; %Wedge angle (should be less than 5 degrees)
L_he = 0.05*Re*D; %Hydrodynamic entry length
L_p =L_he+1; % length of pipe
rho = 997; %density of water kg/m^3
mu = 8.9*1e-4; %dynamic viscosity of water in Pa.s
nu = mu/rho; %kinematic viscosity in m^2/s
L =1 ; %Length of hydrodynamically fully developed region.
V_avg = Re*nu/D; %Average flow velocity in m/sec
V_max =2*V_avg; %maximum velocity in m/sec
del_P = (8*V_avg*mu*L)/(r^2); %pressure drop
k_P =del_P/rho; %kinematic pressure drop m^2/sec^2
% vertices
V_0 =[0 0 0];
V_1 =[0 r*cosd(theta/2) r*sind(theta/2)];
V_2 =[0 r*cosd(theta/2) -r*sind(theta/2)];
V_3 =[L_p 0 0];
V_4 =[L_p r*cosd(theta/2) r*sind(theta/2)];
V_5 =[L_p r*cosd(theta/2) -r*sind(theta/2)];
% creating blockMeshDict file for openFoam
header{1}='/*--------------------------------*- C++ -*----------------------------------*\';
header{2}=' ========= | |';
header{3}=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
header{4}=' \\ / O peration | Website: https://openfoam.org';
header{5}=' \\ / A nd | Version: 8';
header{6}=' \\/ M anipulation |';
header{7}='\*---------------------------------------------------------------------------*/';
f1 = fopen('test.txt','wt');
for i =1 :length(header)
fprintf(f1,'%s\n ',header{i});
end
fprintf(f1,'\n');
fprintf(f1,'%s\n','FoamFile');
fprintf(f1,'%s\n','{');
fprintf(f1,'%16s% 10s\n','version' ,' 2.0;');
fprintf(f1,'%16s %12s\n',' format' ,' ascii;');
fprintf(f1,'%16s %17s\n','class ','dictionary;');
fprintf(f1,'%16s %20s\n','object', 'blockMeshDict;');
fprintf(f1,'%s\n\n',' }');
fprintf(f1,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%s\n','convertToMeters 1;');
fprintf(f1,'\n');
% vertices
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%2s\n','(');
fprintf(f1,'\t(%d %d %d)\n',V_0);
fprintf(f1,'\t(%d %d %d)\n',V_1);
fprintf(f1,'\t(%d %d %d)\n',V_2);
fprintf(f1,'\t(%d %d %d)\n',V_3);
fprintf(f1,'\t(%d %d %d)\n',V_4);
fprintf(f1,'\t(%d %d %d)\n',V_5);
fprintf(f1,'%2s\n\n',');');
% defining the block
z=[300 50 1]; %no of cells in X,Y and Z direction
d=[0.2 1 1 ]; %simple grading
fprintf(f1,'blocks \n (\n');
fprintf(f1,'\t%s (%d %d %d) %s (%d %d %d)','hex (0 3 5 2 0 3 4 1)',z,'simpleGrading',d);
fprintf(f1,'\n\n%2s',');');
fprintf(f1,'\n%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n','arc 1 2 (0 0.01 0)');
fprintf(f1,'%s\n','arc 4 5 (3.100000e+00 0.01 0)');
fprintf(f1,'%s\n',');');
% defining boundaries
fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
%inlet
fprintf(f1,'\t%5s\n','inlet','{');
fprintf(f1,'\n%17s %s\n', 'type', 'patch;');
fprintf(f1,'%18s \n\t %s\n','faces','(');
fprintf(f1,'%24s \n%11s \n\n%15s\n\n','(0 1 2 0)',');','}');
% outlet
fprintf(f1,'\t%5s\n','outlet','{');
fprintf(f1,'\n%17s %s\n', 'type', 'patch;');
fprintf(f1,'%18s \n\t %s\n','faces','(');
fprintf(f1,'%24s \n%11s \n\n%15s\n\n','(3 5 4 3)',');','}');
% top wall
fprintf(f1,'\t%5s\n','topWall','{');
fprintf(f1,'\n%17s %s\n', 'type', 'wall;');
fprintf(f1,'%18s \n\t %s\n','faces','(');
fprintf(f1,'%24s \n%11s \n\n%15s\n\n','(1 4 5 2)',');','}');
%bottom axis
fprintf(f1,'\t%5s\n','bottomAxis','{');
fprintf(f1,'\n%17s %s\n', 'type', 'empty;');
fprintf(f1,'%18s \n\t %s\n','faces','(');
fprintf(f1,'%24s \n%11s \n\n%15s\n\n','(0 3 3 0)',');','}');
% wedge front
fprintf(f1,'\t%5s\n','wedgeFront','{');
fprintf(f1,'\n%17s %s\n', 'type', 'wedge;');
fprintf(f1,'%18s \n\t %s\n','faces','(');
fprintf(f1,'%24s \n%11s \n\n%15s\n\n','(3 4 1 0)',');','}');
% wedge back
fprintf(f1,'\t%5s\n','wedgeBack','{');
fprintf(f1,'\n%17s %s\n', 'type', 'wedge;');
fprintf(f1,'%18s \n\t %s\n','faces','(');
fprintf(f1,'%24s \n%11s \n\n%15s\n\n','(2 5 3 0)',');','}');
fprintf(f1,');\n\n');
% mergepatchpairs
fprintf(f1,'mergePatchPairs\n');
fprintf(f1,'(\n');
fprintf(f1,');\n\n\n');
fprintf(f1,'// ************************************************************************* //');
fclose(f1);
The code below is the blockMeshDict text file which was generated by Matlab.Then we replace the existing BlockMeshDict file which was present , with the BlockMeshDict file which we created in matlab in the OpenFoam.
BlockMeshDict file :
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.993908e-03 3.489950e-04)
(0 9.993908e-03 -3.489950e-04)
(3.100000e+00 0 0)
(3.100000e+00 9.993908e-03 3.489950e-04)
(3.100000e+00 9.993908e-03 -3.489950e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (100 100 1) simpleGrading (0.1 0.1 1)
);
edges
(
arc 1 2 (0 0.01 0)
arc 4 5 (3.100000e+00 0.01 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
topWall
{
type wall;
faces
(
(1 4 5 2)
);
}
bottomAxis
{
type empty;
faces
(
(0 3 3 0)
);
}
wedgeFront
{
type wedge;
faces
(
(3 4 1 0)
);
}
wedgeBack
{
type wedge;
faces
(
(5 3 0 2)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Then we will give the initial and boundary conditions for pressure and velocity by editing the p and U file in 0 folder.
Velocity file :
Here we will give the average velocity calculated above as initial velocity at time t=0.
*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.0937 0 0);
}
outlet
{
type zeroGradient;
}
topWall
{
type noSlip;
}
bottomAxis
{
type empty;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
}
// ************************************************************************* //
Pressure file :
//*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
topWall
{
type zeroGradient;
}
bottomAxis
{
type empty;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //
TransportProperties file :
Here we change the existing viscosity with the kinematic viscosity of water.
transportProperties file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9268e-07;
// ************************************************************************* //
ControlDict file :
Here we give the time step in such a way that the solution should be stable and also it should satisfy the courant number condition(CFL<1).
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 40;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Now after setting up the problem by giving the initial and boundary conditions ,the next thing to do is meshing by using blockMesh and after that we solve the problem.Here we use icoFoam solver, which solves the incompressible laminar Navier-Stokes equations using the PISO algorithm. The code is inherently transient.
Post processing and plotting the results :After meshing by using blockMesh and running the case with icoFoam solver the next thing to do is to post process the results by using paraFoam in paraview and plotting the results by using plotover line command.Let's have a look at velocity distribution along the length of the pipe.
Velocity distribution at the inlet of the pipe :

Velocity distribution at the outlet of the pipe :
Here in the below plot we can see the velocity distribution and boundary layer region in the fully delveloped flow region.

Validation of Results :
- Validation of Hydro-dynamic length with the numerical result :
Now let us validate Hydro-dynamic length from Hagen poiseuille's flow with the numerical result.This is done by verifying the velocity profile,which we calculate at different points along the flow direction .
Velocity profile at X= 0.1m
It is the velocity profile when the fluid enters and cross the inlet of the pipe.

Velocity profile at X= 1m

Velocity profile at X= 2.1m
This is the velocity profile at the hydrodynamic length i.e, fully developed.

Velocity profile at X= 3.1m
It is the velocity profile at the end of the pipe i.e, at the outlet.

From the above plots we can observe that the velocity profile changes along the the flow direction continuously as the fluid progress through the pipe.As the fluid approaches the hydrodynamic length the maximum velocity goes on increasing and you can also observe that the velocity profile also becomes progressively more parabolic.When once the fluid reaches the hydrodynamic length i.e,at X =2.1, the flow becomes fully developed.And once it cross the hydodynamic length the flow becomes fully parabolic and we can notice that the velocity profile is almost constant till it reaches the outlet.You can also observe from the plots that the maximum velocity is at the outlet i.e, 0.1867 m/sec.
- Validation of the fully developed flow velocity profile and maximum velocity with its analytical profile :
We validate the velocity profile of the fully developed flow with its analytical profile by comparing both the profiles, where we take the velocity profile data at the outlet from OpenFoam and take it in a spead sheet and transfer it to matlab and compare this profile with the analytical profile which is generated from the below formula in the matlab.

Matlab code to validate velocity profile :
clear all
close all
clc
r = linspace(0 ,0.01,1001); % radial distance along Y direction
v_avg =0.093731193580742; %average velocity
R= 0.01; %Maximum raidus of the pipe.
for i =1:length(r)
U_analytical(i) = 2*v_avg*(1-(r(i)^2/R.^2)); %Analytical profile
end
velocity_o = xlsread('result.xlsx'); %Numerical data imported from OpenFoam
U_numerical =velocity_o(1:end,1); % seperated velocity data from the whole data
%plotting analytical velocity profile vs numerical velocity profile
plot(r,U_analytical,'r',r,U_numerical,'g--')
legend('U Analytical','U Numerical')
title('Analytical velocity profile vs Numerical velocity profile for fully developed flow')
xlabel('radial distance along Y direction')
ylabel('Velocity magnitude in (m/sec)')
Result :

From the above plot it is clear that our numerical values is in good agreement with the analytical values.Let's compare the maximum velocity of numerical solution with analytical solution.
Numerical value Analytical value %Error
0.1867 m/sec 0.1875 m/sec 0.3%
Hence the difference is in the order of 1e-3 the results that obtained in our numerical simulation are in good agreement with the analytical solution and are acceptable.
- Validation of pressure drop for fully developed flow :
Now lets us validate the pressure drop of the fully developed flow by comparing the pressure drop in the numerical simulation after the flow is fully developed with the analytical value calculated from the Hagen poiseuille's flow equation.

From the above plot we can observe that the numerical value of maximum kinematic pressure drop after the flow is fully developed is 0.0067 m^2/sec^2.Now let's compare this with the analyitcal value of kinematic pressure drop.Which is k_P =del_P/rho= 6.6737/997 = 0.0067 m^2/sec^2.As the numerical and analytical pressure drop are almost same,and hence we can conclude that numerical simulation are in close agreement with the analytical solution.
- Post-process of Shear stress and validation of wall shear stress for fully developed flow :

The above plot is the shear stress plot calculated at X=3.1 along radial direction .From the plot we can observe that the maximum shear stress is near the pipe walls i.e, 0.0332 pa
The Wall shear stress in a pipe for a fully developed flow is calculated analytically by form the formula,
Here τau_w is the maximum shear stress for a fully developed flow i.e, near the wall, R is the radius of the pipe.
therefore tau_w =(8*8.9*1e-4*0.0937)/(2*0.01) = 0.0334 Pa
Therefore we can finally conclude that numerical simulation value of shear stress are in close agreement with the analytical solution. Hence shear stress for fully developed flow also validated.
Conclusion :
- The Hagen poiseuille's flow was succesfully simulated numerically using OpenFoam and the results are found to be in good agreement with the analytical solution calculated from Hagen poiseuille's equations.
- The Validation of the numerical results also succesfully done,where we found that once flow cross the hydodynamic length the flow becomes fully parabolic and we noticed that the velocity profile is almost constant till it reaches the outlet.We observed from the plots that the maximum velocity is at the outlet i.e, 0.1867 m/sec.We have validated the velocity profile of the fully developed flow with its analytical profile by comparing both the profiles.
- From the plot of kinematic pressure we can observe that the numerical value of maximum kinematic pressure drop after the flow is fully developed is 0.0067 m^2/sec^2.We found that numerical and analytical pressure drop are almost same ,and we have concluded that numerical simulation are in close agreement with the analytical solution.
- we observed that numerical simulation value of shear stress are in close agreement with the analytical solution.Hence shear stress for fully developed flow also validated.
- Hence we can finally conclude that CFD simulations with correct operating conditions, exact boundary and initial conditions ,assumptions and correct approximations, we can get similar and close results that of analytical solutions.