Week 11 - Simulation of Flow through a pipe in OpenFoam
Question:
In this challenge, you will simulate an axis-symmetric flow by applying the wedge boundary condition. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100. You can do a transient or steady state simulation and you will validate your results with Hagen- Poiseuille's equations for the fully developed flow.
Aim:
write code in Matlab to automate the generation of the blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc, and output a fully functional blockMeshDict file.
Expected Results
- Validate Hydro-dynamic length with the numerical result
- Validate the fully developed flow velocity profile with its analytical profile
- Validate maximum velocity and pressured drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
EQUATION:
In standard fluid-kinetics notation.
where:
Δp = pressure difference between the two ends,
L = length of pipe,
μ = dynamic viscosity,
Q = volumetric flow rate,
R = pipe radius,
A = cross section of pipe.
The equation does not hold close to the pipe entrance.
Given parameters.
- Wedge angle less than 5 degrees,
- Flow is axis symmetrical.
- Fluid is incompressible and the flow is laminar.
Expected Velocity Profile.
Assuming the following:
- The distance the fluid has to travel to become fully developed is given by the Entry Length formula.
- Given Reynolds number is 2100. This implies flow is laminar.
- Hence Entry length formula for Laminar is given by:Ld=06.ReLd=0.06.Re
- where L is effective length; d is dia of pipe and Re is reyonolds number.
- Assuming the dia of pipe as 10mm; substituting the values, we get L = 1.2m;
Matlab program:
%Automatically generate blockmeshDict file for a wedge angle provided by user.
clc
clear all
close all
% variable values
D = 0.01;
r = D/2;
L = 1.2;
theta = input(' input the value of wedge angle=');
% String for header,
h1 =' /*---------------------------------*- C++ -*-------------------------------------------------*\ ';
h2 =' | ========== | |';
h3 =' | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |';
h4 =' | \\ / O peration | Version: 7 |';
h5 =' | \\ / A nd | Web: www.OpenFOAM.org |';
h6 =' | \\ / M anipulation | |';
h7 =' \*------------------------------------------------------------------------------------------*/ ';
h8 =' FoamFile';
h9 ='{';
h10 =' version 2.0;';
h11 =' format ascii;';
h12 =' class dictionary;';
h13 =' object blockMeshDict;';
h14 ='}';
h15 =' // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * \\';
conv ='convertToMeters 1;';
sp1 = blanks(5); %sp = space
sp2 = blanks(10);
sp3 = blanks(15);
v_0 = [0 0 0];
v_1 = [r*sind(theta/2) r*cosd(theta/2) 0];
v_2 = [-r*sind(theta/2) r*cosd(theta/2) 0];
v_3 = [0 0 L];
v_4 = [r*sind(theta/2) r*cosd(theta/2) L];
v_5 = [-r*sind(theta/2) r*cosd(theta/2) L];
block ='hex (0 0 1 2 3 3 4 5) (1 50 200) simpleGrading (1 0.5 1)';
edge1 ='arc 1 2';
edge2 ='arc 4 5';
e1 = [0 r 0];
e2 = [0 r L];
inlet = '(0 0 2 1)';
outlet = '(3 4 5 3)';
wedge1 = '(0 1 4 3)';
wedge2 = '(3 5 2 0)';
wall =' (4 1 2 5)';
axis =' (0 3 3 0)';
final = '// ************************************************************************************************************************************************** //';
f1 = fopen('blockMeshDict','wt');
fprintf(f1,'%s\n',h1);
fprintf(f1,'%s\n',h2);
fprintf(f1,'%s\n',h3);
fprintf(f1,'%s\n',h4);
fprintf(f1,'%s\n',h5);
fprintf(f1,'%s\n',h6);
fprintf(f1,'%s\n',h7);
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n',h9);
fprintf(f1,'%s\n',h10);
fprintf(f1,'%s\n',h11);
fprintf(f1,'%s\n',h12);
fprintf(f1,'%s\n',h13);
fprintf(f1,'%s\n',h14);
fprintf(f1,'%s\n',h15);
fprintf(f1,'\n');
fprintf(f1,'%s\n',conv);
fprintf(f1,'\n');
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s (%f %f %f) \n',sp1,v_0(1),v_0(2),v_0(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_1(1),v_1(2),v_1(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_2(1),v_2(2),v_2(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_3(1),v_3(2),v_3(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_4(1),v_4(2),v_4(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_5(1),v_5(2),v_5(3));
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s \n',sp1,block);
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n','edges');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s %s (%f %f %f) \n',sp1,edge1,e1(1),e1(2),e1(3));
fprintf(f1,'%s %s (%f %f %f) \n',sp1,edge2,e2(1),e2(2),e2(3));
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n','boundary');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s %s \n',sp1,'inflow');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type patch;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,inlet);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'outflow');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type patch;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,outlet);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'front');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type wedge;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,wedge1);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'back');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type wedge;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,wedge2);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'wall');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type wall;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,wall);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'axis');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type empty;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,axis);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n','mergePatchPairs');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n',final);
fclose(f1);
- VELOCITY BOUNDARY CONDITIONS.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ 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
{
inflow
{
type fixedValue;
value uniform (0 0 0.1869);
}
outflow
{
type zeroGradient;
}
wall
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
- PRESSURE BOUNDARY CONDITIONS.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inflow
{
type zeroGradient;
}
outflow
{
type fixedValue;
value uniform 0.0638749;
}
wall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Meshed Geometry:
The isometric view of the mesh is shown. Mesh grading is given only for the walls.
Velocity profiles:
Shear Stress :-
shear stress = (dynamic Viscosity)*(du/dy)
Conclusion:
We simulated an axis symetrical wedge as the simulation of the whole pipe requires high processing. The velocity profiles changes as the flow are fully developed for the entry length of the flow is 1.2m. And after this distance from the inlet, the velocity we get parabola type shape and doesn't change any more. Shear Stress increases with the friction on pipe walls because velocity is maximum at the center of the pipe.