Week 11 - Simulation of Flow through a pipe in OpenFoam
Week 11 - Simulation of Flow through a pipe in OpenFoam
GIVEN PROBLEM:
1. simulate an axi-symmetric flow by applying the wedge boundary condition.
2. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.
3. You can do a transient or steady state simulation and you will validate your results with the Hagen- Poiseuille's equations for the fully developed flow.
4. You will also write code in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.
5. You can create a wedge from a hex block by collapsing its bottom face (by merging both sets of adjacent vertices)
6. Wedge angle should be less than 5 degree
7. You can assume the pipe diameter
THEORY:
OpenFOAM (acronym for Open Source Field Operation and Manipulation) is a C++ toolbox for the development of customised numerical solvers, and pre/postprocessing utilities for the solution of continuum mechanics problems, most prominently the CFD problems.
icoFoam solves the incompressible laminar Navier-stokes equations using the PISO algorithm. The code is inherently transient, requiring an initial condition (such as zero velocity) and boundary conditions.
The problem statement for this challenge; we have a pipe geometry. But, since, simulating a complete pipe can be resource extensive and the fact that pipes are '"axis-symmetrical", we can deduce it into a wedge element and still obtain the same result.
The geometric parameters are already known and hence, we only need to define the angle of the wedge formed. Also as we are treating the flow to be laminar, the Reynolds number is also known. The assumptions are:
- ρ=997kg/m3
- μ=8.89×10-4kg/ms
- Re=2100
From Reynold's number formula, we have,
Re=ρvD/μ→vD=1.8726×10-3
We need to choose velocity and diameter which can satisfy the above condition. Also, we need to ensure that pipe length is greater than the entrance length.
Le=0.006ReD→Le=12.6D & L>Le
Thus, we arrive at,
D=0.025m
L=0.315m
Now, we only need to select the wedge angle and thus obtain a blockMeshDict file for the same. In the blockMeshDict file, we'll apply the wedge boundary conditions and study for the case of 4 degree. The file would be autogenerated using MATLAB and appropriate initial conditions would be applied to simulate the flow.
Expected Velocity Profile.
HOW TO FIND CHANGE IN PRESSURE THEORATICALLY:
In standard fluid-kinetics notation.
Δp is the pressure difference between the two ends,
L is the length of pipe,
μ is the dynamic viscosity,
Q is the volumetric flow rate,
R is the pipe radius,
A is the cross section of pipe.
MATLAB CODE:
% Script to generate blockMeshDict for any wedge angle for a pipe with diameter D and lenght L
clear all
close all
clc
% Input parameters
theta = 4; % required wedge angle
Re = 2100; % reynolds number
D = 0.025;% Diameter of pipe in meters
L = 0.06*Re*D; % length of pipe in meters
R = D/2; % radius
% co-ordinates for the vertices will be as follows using trignometric relations
c0 = [0 0 0]; % point of origin;
c1 = [0 R*(cosd(theta/2)) R*(sind(theta/2))];
c2 = [0 R*(cosd(theta/2)) -R*(sind(theta/2))];
c3 = [L 0 0];
c4 = [L R*(cosd(theta/2)) R*(sind(theta/2))];
c5 = [L R*(cosd(theta/2)) -R*(sind(theta/2))];
%%strings for blockMeshDict file
h1 = ('/*--------------------------------*- C++ -*----------------------------------*\');
h2 = ('| ========= |');
h3 = ('| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |');
h4 = ('| \\ / O peration | Version: 10 |');
h5 = ('| \\ / A nd | Web: www.OpenFOAM.org |');
h6 = ('| \\/ M anipulation | |');
h7 = ('\*---------------------------------------------------------------------------*/');
h8 = ('// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
%% Generating a blockMeshDict file
% Header
name = sprintf('blockMeshDict');
file = fopen(name,'w');
fprintf(file,'%s\n',h1);
fprintf(file,'%s\n',h2);
fprintf(file,'%s\n',h3);
fprintf(file,'%s\n',h4);
fprintf(file,'%s\n',h5);
fprintf(file,'%s\n',h6);
fprintf(file,'%s\n',h7);
fprintf(file,'FoamFile\n{\n');
fprintf(file,'%16s \t 2.0;\n','version');
fprintf(file,'%16s \t ascii;\n','format');
fprintf(file,'%16s \t dictionary;\n','class');
fprintf(file,'%16s \t blockMeshDict;\n}\n','object');
fprintf(file,'%s\n\n',h8);
% vertices
fprintf(file,'convertToMeters 1;\n\n');
fprintf(file,'vertices\n(\n');
fprintf(file,'\t(%d %d %d)\n',c0(1),c0(2),c0(3));
fprintf(file,'\t(%d %d %d)\n',c1(1),c1(2),c1(3));
fprintf(file,'\t(%d %d %d)\n',c2(1),c2(2),c2(3));
fprintf(file,'\t(%d %d %d)\n',c3(1),c3(2),c3(3));
fprintf(file,'\t(%d %d %d)\n',c4(1),c4(2),c4(3));
fprintf(file,'\t(%d %d %d)\n);\n\n',c5(1),c5(2),c5(3));
fprintf(file,'blocks\n(\n');
fprintf(file,'\t hex (0 3 5 2 0 3 4 1) (400 40 1) simplegrading (1 0.1 1)\n');
fprintf(file,');\n');
fprintf(file,'edges\n(\n');
fprintf(file,'\t arc 1 2 (%d %d %d)\n',0, R, 0);
fprintf(file,'\t arc 4 5 (%d %d %d)\n',L, R, 0);
fprintf(file,');\n');
fprintf(file,'boundary\n(\n');
fprintf(file,'\t Inlet\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type patch;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (0 1 2 0)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n\n');
fprintf(file,'\t Outlet\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type patch;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (3 4 5 3)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,'\t Top\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type patch;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (1 2 5 4)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,'\t wedgeback\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type wedge;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (0 1 4 3)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,'\t wedgefront\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type wedge;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (0 2 5 3)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,'\t axis\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type empty;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (0 3 3 0)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,');\n');
fprintf(file,'mergePatchPairs\n');
fprintf(file,'(\n');
fprintf(file,');\n');
fprintf(file,h8);
fclose(file);blockMeshDict file:
*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 10 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 10.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 1.249239e-02 4.362437e-04)
(0 1.249239e-02 -4.362437e-04)
(3.150000e+00 0 0)
(3.150000e+00 1.249239e-02 4.362437e-04)
(3.150000e+00 1.249239e-02 -4.362437e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (400 40 1) simplegrading (1 0.1 1)
);
edges
(
arc 1 2 (0 1.250000e-02 0)
arc 4 5 (3.150000e+00 1.250000e-02 0)
);
boundary
(
Inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
Outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
Top
{
type patch;
faces
(
(1 2 5 4)
);
}
wedgeback
{
type wedge;
faces
(
(0 1 4 3)
);
}
wedgefront
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
controlDict file:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 10 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 10.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 1;
deltaT 0.02;
writeControl timeStep;
writeInterval 5;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
velocity file:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 10 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 10.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.075 0 0);
boundaryField
{
Inlet
{
type fixedValue;
value uniform (0.075 0 0);
}
Outlet
{
type zeroGradient;
}
axis
{
type empty;
}
Top
{
type fixedValue;
value uniform (0 0 0);
}
wedgefront
{
type wedge;
}
wedgeback
{
type wedge;
}
}
// ************************************************************************* //
pressure file:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 10 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 10.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.0125;
}
axis
{
type empty;
}
Top
{
type zeroGradient;
}
wedgefront
{
type wedge;
}
wedgeback
{
type wedge;
}
}
// ************************************************************************* //
physicalproperties:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 10.0;
format ascii;
class dictionary;
location "constant";
object physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.892e-6;
// ************************************************************************* //
output:
Mesh file:
velocity profiles:
velocity @ 1m:
pressure Distribution:
Shear stress Distribution:
Results:
maximum velocity pressure loss max shear stress
numerical: .083 (m/s) .019 pa .053 pa
Theoratical: .083(m/s) .019 pa .051 pa
CONCLUSION:
- We simulate a axi symetrical wedge as the simulation of the whole pipe requires high processing.
- The velocity profiles changes until the flow is fully developed.
- The velocity profile obtained is parabolic in nature with max velocity attained at the centre of the pipe.
- Shear Stress increases the friction on the walls of the pipe because of the which the velocity is maximum at the center of the pipe.
- The Shear stress plot is obtained after post processing and it is found that the shear stress is maximum at the walls and zero at the centre of the pipe.