Hagen Poiseuille s Flow

Google drive link to all the blockMeshDict file for various angles:


Google drive link to all the plots and source code:




1. Reynolds Number = 1000 (For Laminar Flow)
2. Temperature= 20 degree celsius

at this temperature (From Saturated Water Tables Ref. Cengel and Cimbala)

Kinematic Viscosity = 1.004008e-06 m2/s
Density = 998 kg/m3

3. Length of pipe = 1 m
4. Diameter of pipe = 20 mm


Analytical Results:


1. Reynolds Number = (rho*V*D/mu)

2.Del_p=(32*mu*V*L/D2)   %Hagen Poiseuille's Equation 

We get: 

Pressure Drop = 4.024064e+00 pa
Kinematic Inlet pressure = 4.032128e-03 m2/s2
Kinematic Outlet pressure = 0 m2/s2
Max velocity = 1.004008e-01 m/s
Avg Velocity = 5.020040e-02 m/s 


1. The pressures(p) used to calculate the kinematic pressures(p/rho) were Gauge Pressures at inlet and outlet.

2. A pressure driven approach was used i.e the inlet and outlet kinematic pressures were fixed in initial conditions and the velocity was left to float.

3. The Transport Model used was Newtonian .


CFD Results:

1. Velocity Profile:

The Velocity profile obtained at the inlet and outlet are shown below (the plots shown are for wedge angle = 2 degrees )

The reflect filter in paraview was used to model the complete velocity profile across the diameter of pipe.


1. We see that the profile is fully developed and is a parabola which matches with the Hagen Poiseuille's Result.

2. Also the maximum velocity obtained from the above plots was 0.09912 m/s which has only 1.27 % error i.e almost equal to the result obtained by Hagen Poiseuille equation.


Further comparison for Maximum Velocity and Average Velocity with changing Wedge Angle :

We see that the average velocity remains mostly unchanged with changing the wedge angle while the maximum velocity approaches towards the analytical value as a result of decreasing the wedge angle.


2. Pressure Drop (Loss of Head in Pipe)

The Pressure Drop is the pipe along the length is linear.


3. Wall Shear Stress

As the flow is laminar, so the turbulence model used to compute the Wall Shear stress was Laminar.


4. Accuracy vs Wedge Angle

The percentage error in maximum velocity for decreasing wedge angle was plotted and it is seen that the percentage error decreases as we decrease the wedge angle.


5. Simulation Time vs Wedge Angle

It is seen that the simulation time rises at first when we decrease the wedge angle but it falls subsequently after that as we further decrease the wedge angle.



Program in MATLAB to generate blockMeshDict file for any wedge angle and grading schemes:

Following is the program to generate the blockMeshDict file automatically to generate the computational mesh

clear all
close all

theta=2;  %In degrees %Wedge angle
grading=0.25; %Grading factor along the walls of cylinder 

b1=400; %cells along the length of cylinder
b2=12; %cells along the radius of cylinder

%vertices of the wedge
v0=[0 0 0];
v1=[length 0 0];
v2=[length radius*cosd(theta/2) -radius*sind(theta/2)];
v3=[0 radius*cosd(theta/2) -radius*sind(theta/2)];
v4=[0 radius*cosd(theta/2) radius*sind(theta/2)];
v5=[length radius*cosd(theta/2) radius*sind(theta/2)];

s1='/*--------------------------------*- C++ -*----------------------------------*\';
s2='| =========                 |                                                 |';
s3='| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |';
s4='|  \\    /   O peration     | Version:  4.1                                   |';
s5='|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |';
s6='|    \\/     M anipulation  |                                                 |';

s8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

%wedge block
% hex (0 1 2 3 0 1 5 4) (80 20 1) simpleGrading(1 1 1)



fprintf(f1,'%16s \t 2.0;\n','version');
fprintf(f1,'%16s \t ascii;\n','format');
fprintf(f1,'%16s \t dictionary;\n','class');
fprintf(f1,'%16s \t blockMeshDict;\n}\n','object');

fprintf(f1,'convertToMeters 0.001;\n\n');
fprintf(f1,'\t(%d %d %d)\n',v0(1),v0(2),v0(3));
fprintf(f1,'\t(%d %d %d)\n',v1(1),v1(2),v1(3));
fprintf(f1,'\t(%d %d %d)\n',v2(1),v2(2),v2(3));
fprintf(f1,'\t(%d %d %d)\n',v3(1),v3(2),v3(3));
fprintf(f1,'\t(%d %d %d)\n',v4(1),v4(2),v4(3));
fprintf(f1,'\t(%d %d %d)\n);\n\n',v5(1),v5(2),v5(3));

fprintf(f1,'\t hex (0 1 2 3 0 1 5 4) (%d %d %d)\n',b1,b2,1);
fprintf(f1,'\t simpleGrading\n\t (1 %d 1)\n);\n\n',grading);


fprintf(f1,'\t%s\n\t{\n\t\t type %s;\n\t\t faces \n\t\t (\n\t\t\t (%d %d %d %d) \n\t\t );\n\t}\n\n','inlet','patch',0,4,3,0);
fprintf(f1,'\t%s\n\t{\n\t\t type %s;\n\t\t faces \n\t\t (\n\t\t\t (%d %d %d %d) \n\t\t );\n\t}\n\n','outlet','patch',1,2,5,1);
fprintf(f1,'\t%s\n\t{\n\t\t type %s;\n\t\t faces \n\t\t (\n\t\t\t (%d %d %d %d) \n\t\t );\n\t}\n\n','Wall','wall',2,3,4,5);
fprintf(f1,'\t%s\n\t{\n\t\t type %s;\n\t\t faces \n\t\t (\n\t\t\t (%d %d %d %d) \n\t\t );\n\t}\n\n','axis','empty',0,1,1,0);
fprintf(f1,'\t%s\n\t{\n\t\t type %s;\n\t\t faces \n\t\t (\n\t\t\t (%d %d %d %d) \n\t\t );\n\t}\n\n','wedgeFront','wedge',0,1,2,3);
fprintf(f1,'\t%s\n\t{\n\t\t type %s;\n\t\t faces \n\t\t (\n\t\t\t (%d %d %d %d) \n\t\t );\n\t}\n\n','wedgeBack','wedge',0,1,5,4);





Projects by Utkarsh Garg

The backward facing step flow project was performed as a part of skill lync challenge. The simulation was performed on three base grid sizes: 1. 3e-4m 2. 2.5e-4m 3. 2e-4m The above grid sizes were chosen due to lack of computation power.   The Results a Read more

The simulation was run by refining the grid size successively in three steps as follows: Grid Size 1. 2e-4 m 2. 1.5e-4 m 3. 1.2e-4 m   The post-processing results of the three cases are as follows: 1. 2e-4 m Velocity Magnitude Mass Flow Rate Spec Read more

CENTRIFUGAL PUMP DESIGN AND PERFORMANCE ANALYSIS   1. Geometry: PUMP   IMPELLER      2. Performance Analysis of the Centrifugal Pump Performance Analysis of a Centrifugal Pump Velocity normal to face (Outlet Velocity) [m/ Read more

1. Geometry: All dimensions in mm.   2. Boundary Conditions: Total pressure inlet = 4 bar Static pressure outlet = 1 bar Real wall = all inner surfaces of the cylinder, valve, and intake.   3. Grid dependence study: Goal (Value) Design Point Read more

Project Report   Assumptions (same as in previous challenge): 1. Reynolds Number = 1000 (For Laminar Flow) 2. Temperature= 20 degree Celsius 3. Length of pipe = 1 m 4. Diameter of pipe = 20 mm 5. At 20 Degrees Celsius (Ref: Saturated Water tables from Cengel Read more

Google Drive Link for all sorce code files and functions files https://drive.google.com/open?id=1UE-_9nn7VMuXmUqr3XWFgwbxaaGaNpdX   Outputs: 1. Jacobi Method The program also showed the following output in the command window \"The Jacobi Solver could not co Read more

Google drive link to all the contours and plot: https://drive.google.com/open?id=1fBI5Zcassvv2ut9o2Qsy7ktrIq3LQx10 The link contains all the plots for various angles of attack and the pressure and velocity distribution for all cases.   The airfoil used was NACA Read more

Google drive link containing source code for both functions and all the plots https://drive.google.com/open?id=1q-2do7AV4YP2U5ilL54n37U3CXuNoeNK   Observations and Conclusions: 1. Normalized Mass flow rate We see that in case of conservative form of governing Read more


The End