Week 12 - Symmetry vs Wedge vs HP equation
INTRODUCTION: Continuing from the previous part of this challenge, we calculated numerical solution of the pipe flow simulation using the type ‘wedge’, type wedge is only accurate till small angles, say less than 5.If we increase the angle more than 5degree, eg 7degree keeping the type as wedge and perform the simulation, some of the terms start becoming undefined or not clearly defined. That said, we use type ‘symmetry’ for greater angles .
PROCEDURE:
The steps are just the same as week 11 challenge,
-Finding the vertex co-ordinates and typing code in octave.
-copying the blockmesh file parsed from octave code in blockMesh directory .
-Changing the type ‘wedge’ to type ‘symmetry’ .
-writing the delta t in controldict file, here I take it as dt=0.01sec and time of simulation as 20sec.
-running the icofoam solver.
-Running the parafoam command and plotting the results.
Here I show the results of the fluid flow development in pipe.
There is just the change in the front and back ‘type’ in the blockMesh directory as well as in the time step 0 files of pressure and velocity. Rest all is the same while performing the simulation.
Octave code for vertex and blockMeshDirectory:
%input theta in degrees
theta=input("enter theta");
%write simplgrading factor
%Simple Grading factor =SG
SG_x=input("enter x grading factor");
SG_y=input("enter y grading factor ");
SG_z=input("enter z grading factor");
radius=0.01;
Length=3.1;
%calculating vertex co-ordinates
vertex_0=[0 0 0];
vertex_1=[radius*cos((theta/2)*(pi/180)) -1*radius*sin((theta/2)*(pi/180)) 0];
vertex_2=[radius*cos((theta/2)*(pi/180)),radius*sin((theta/2)*(pi/180)),0];
vertex_3=[0,0,Length];
vertex_4=[radius*cos((theta/2)*(pi/180)),-1*radius*sin((theta/2)*(pi/180)) ,Length];
vertex_5=[radius*cos((theta/2)*(pi/180)),radius*sin((theta/2)*(pi/180)) ,Length];
f2=fopen('blockfoam_45s.txt','w');
fprintf(f2,'vertices \n');
fprintf(f2,'( \n');
fprintf(f2,' (%d %d %d) \n',vertex_0(1),vertex_0(2),vertex_0(3));
fprintf(f2,' (%1.7f %1.7f %d) \n',vertex_1(1),vertex_1(2),vertex_1(3));
fprintf(f2,' (%1.7f %1.7f %d) \n',vertex_2(1),vertex_2(2),vertex_2(3));
fprintf(f2,' (%d %d %d) \n',vertex_3(1),vertex_3(2),vertex_3(3));
fprintf(f2,' (%1.7f %1.7f %d) \n',vertex_4(1),vertex_4(2),vertex_4(3));
fprintf(f2,' (%1.7f %1.7f %d) \n',vertex_5(1),vertex_5(2),vertex_5(3));
fprintf(f2,'); \n');
%starting to write block
fprintf(f2,'\n');
fprintf(f2,'blocks \n');
fprintf(f2,'( \n');
fprintf(f2,' hex (0 1 2 0 3 4 5 3) (80 1 200) ')
fprintf(f2,'simpleGrading (%2.1f %2.1f %2.1f) \n',SG_x,SG_y,SG_z);
fprintf(f2,'); \n');
%starting to write edge
fprintf(f2,'\n');
fprintf(f2,'edges \n');
fprintf(f2,'( \n');
fprintf(f2,' arc 1 2 (0.01 0 0) \n');
fprintf(f2,' arc 4 5 (0.01 0 3.1) \n');
fprintf(f2,'); \n');
%starting to write boundary
fprintf(f2,'\n');
fprintf(f2,'boundary \n');
fprintf(f2,'( \n');
%for front
fprintf(f2,' front \n');
fprintf(f2,' { \n');
fprintf(f2,' type symmetry; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (5 2 0 3) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for back
fprintf(f2,' back \n');
fprintf(f2,' { \n');
fprintf(f2,' type symmetry; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (0 1 4 3) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for tank wall
fprintf(f2,' tankWall \n');
fprintf(f2,' { \n');
fprintf(f2,' type wall; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (5 4 1 2) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for inlet
fprintf(f2,' inlet \n');
fprintf(f2,' { \n');
fprintf(f2,' type patch; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (0 2 1 0) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for outlet
fprintf(f2,' outlet \n');
fprintf(f2,' { \n');
fprintf(f2,' type patch; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (3 5 4 3) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for axis
fprintf(f2,' axis \n');
fprintf(f2,' { \n');
fprintf(f2,' type empty; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (0 3 3 0) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
fprintf(f2,');');
%code
fclose(f2);blockMesh files for theta 10, 25 , 45:
theta 10:
/*--------------------------------*- 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 dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0.0099619 -0.0008716 0)
(0.0099619 0.0008716 0)
(0 0 3.1)
(0.0099619 -0.0008716 3.1)
(0.0099619 0.0008716 3.1)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (80 1 200) simpleGrading (0.1 1.0 1.0)
);
edges
(
arc 1 2 (0.01 0 0)
arc 4 5 (0.01 0 3.1)
);
boundary
(
front
{
type symmetry;
faces
(
(5 2 0 3)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
tankWall
{
type wall;
faces
(
(5 4 1 2)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
theta25:
/*--------------------------------*- 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 dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0.0097630 -0.0021644 0)
(0.0097630 0.0021644 0)
(0 0 3.1)
(0.0097630 -0.0021644 3.1)
(0.0097630 0.0021644 3.1)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (80 1 200) simpleGrading (0.1 1.0 1.0)
);
edges
(
arc 1 2 (0.01 0 0)
arc 4 5 (0.01 0 3.1)
);
boundary
(
front
{
type symmetry;
faces
(
(5 2 0 3)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
tankWall
{
type wall;
faces
(
(5 4 1 2)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
theta45:
/*--------------------------------*- 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 dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0.0092388 -0.0038268 0)
(0.0092388 0.0038268 0)
(0 0 3.1)
(0.0092388 -0.0038268 3.1)
(0.0092388 0.0038268 3.1)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (80 1 200) simpleGrading (0.1 1.0 1.0)
);
edges
(
arc 1 2 (0.01 0 0)
arc 4 5 (0.01 0 3.1)
);
boundary
(
front
{
type symmetry;
faces
(
(5 2 0 3)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
tankWall
{
type wall;
faces
(
(5 4 1 2)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Boundary conditions are the same
Pressure file at time step 0:
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
tankWall
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Velocity file at time step 0:
/*--------------------------------*- 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
{
inlet
{
type fixedValue;
value uniform (0 0 0.09345);
}
outlet
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
tankWall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Transport Properties are the same as week 11 challenge
Obtaining the results with symmetry BC:
The mesh looks as:
OBSERVATION:
1) Plotting the velocity over the whole length of cylinder.
Image: theta 10:
theta25:
theta45:
CONCLUSION: We can say that velocity becomes almost constant after a certain Length i.e after fully developed flow.
2) Plotting the velocity in the fully developed region, at length 2.7m ( as entrance length lasts till only 2.1m)
Showing the growth of velocity profile from developing to developed flow.
Comparing this velocity in fully developed region that we numerically solved with the analytical solution of hagen poisulle equation. Hagen poisulle equation only works in the fully developed flow region.
theta10:
velocity at 0.005m
velocity at 0.01m:
velocity at 0.1m:
velocity at 0.5m:
velocity at 1.5m:
velocity at 2.7m(fully developed flow):
Comparison of fully developed flow with hagen:
octave code for plotting the above comparison:
clear all
close all
clc
%velocity comparison at length 2.8m when theta is 2
A=csvread('VELOCITY_COMPARISON_WITH_HAGEN10s.csv')
%U_n is numerical velocity and U_a is analytical velocity
y=A(:,1);
U_n=A(:,2);
%velocity by hagen for fully developed flow
R=0.01;
%u_max is 2 times the average velocity
u_max=2*0.09345;
U_a=u_max*(1-((y.^2)/(R^2)));
plot(y,U_n,"color","red");
hold on
plot(y,U_a,"color","blue");
legend("blue- analytical curve, red-numerical curve")
set(gca,"xlabel",text("string","length","fontsize",25));
set(gca,"ylabel",text("string","velocity","fontsize",25));Excel file is saved as: 'VELOCITY_COMPARISON_WITH_HAGEN10S'
theta25:
Velocity for 0.005m:
velocity at 0.1m:
velocity at 1.5m:
velocity at 2.7m(fully developed region):
Comparison with analytical velocity:
theta 45:
velocity at 0.01m:
velocity at 0.1m:
velocity at 1.5m:
velocity at 2.7m(fully developed region):
CONCLUSION: Velocity profile is same by both analytical solution(hagen poiseuille) and numerical approach
Table below shows that error decreases as angle increases,
Type-wedge:
|
Angle |
2 |
4 |
5 |
|
Analytical max velocity |
0.1869 |
0.1869 |
0.1869 |
|
Numerical max velocity |
0.185304 |
0.185309 |
0.185312 |
|
error |
0.001596 |
0.001591 |
0.001588 |
Type-symmetry:
|
Angle |
10 |
425 |
45 |
|
Analytical max velocity |
0.1869 |
0.1869 |
0.1869 |
|
Numerical max velocity |
0.184915 |
0.185495 |
0.185965 |
|
error |
0.001985 |
0.001405 |
0.000935 |
3) Plotting the pressure along the length from 2.4m to 3m in the cylinder , which is again in the fully developed flow domain. We compare this pressure drop with drop as per calculated by the hagen equation .
theta10:
theta25:
theta45:
As the angle increases , the pressure drop becomes more than calculated as per the hagen equation. This can be shown by the table both for type wedge and for type symmetry.
Here: Pn-Pa= error ie Difference between Pressure drop by numerical solution and pressure drop by analytical solution.
Type-wedge:
|
Angle |
2 |
4 |
5 |
|
Pressure at 2.4 |
0.00469851 |
0.00470228 |
0.0047051 |
|
Presure at 3 |
0.000669956 |
0.000670559 |
0.000671009 |
|
Pressure drop |
0.004028554 |
0.004031721 |
0.004034091 |
|
Pressure drop from hagen equation |
0.003992184 |
0.003992184 |
0.003992184 |
|
Pn-Pa |
3.637E-05 |
3.9537E-05 |
4.1907E-05 |
Type-Symmetry:
|
Angle |
10 |
25 |
45 |
|
Pressure at 2.4 |
0.00473332 |
0.00492248 |
0.0054872 |
|
Presure at 3 |
0.000675408 |
0.000702232 |
0.00078316 |
|
Pressure drop |
0.004057912 |
0.004220248 |
0.00470404 |
|
Pressure drop from hagen equation |
0.003992184 |
0.003992184 |
0.003992184 |
|
Pn-Pa |
6.5728E-05 |
2.28E-04 |
7.12E-04 |
CONCLUSION:As the angle increases , the error increases,i.e Pressure drop is more accurate or more near to the value calculated by hagen -poiseuille equation for smaller angles.
4) Plotting the shear stress along with the radial distance .
theta10:
theta25:
theta45:
CONCLUSION: velocity profile is parabolic, so velocity gradient will be a straight line. Hence, this is correct.
5)CONCLUSION: Type-symmetry took some more time to perform the execution then type wedge.
Type-wedge:
|
Angle |
2 |
4 |
5 |
|
Simulation time |
442.41s |
518.16s |
413.19s |
|
Clocktime |
539s |
604s |
464s |
Type-symmetry:
|
Angle |
10 |
25 |
45 |
|
Simulation time |
491.2s |
617.5s |
487.14s |
|
Clocktime |
589s |
689s |
565s |
Hence, i have demonstrated the flow through a pipe for an axi-ymmetric cylinder usimg different angles using the type 'symmetry' instead of wedge. Profiles remains same for all the three angles.