Symmetry vs Wedge vs HP equation

Objective:

• To write a Matlab program that takes an angle as input and generates a blockMesh file for the given angle of 10,25,45.
• To compare the results obtained from solving the Hagen-Poiseuille Equation with simulation B.C. of wedge and symmetry.

Theory:

Using the Hagen-Poiseuille Equation,

Maximum velocity, V_max = 0.3742 m/s.

Pressure difference, ∆P = 2.4644e-4 N/m2.

Kinematic pressure difference, ∆P_k=2.4644e-7 m.

Code:

matlab code for mesh file:

clear all
close all
clc

% input
Re = 2100; % Reynolds number
mu = 8.9e-4; %Dynamic viscosity (n/m2.s)
rho = 999; % density (kg/m3)
v_kin = mu/rho;%kinematic viscosity (m2/s)
d = input('enter the diameter: '); %Diameter (m)
theta = input('enter value of theta:  '); % in degree
theta = theta/2;
theta_rad = theta*pi/180; % in radians
r = d/2;
V_avg = Re*v_kin/d;
V_max = V_avg*2;
Le = 0.06*Re*d %entry length of pipe
L = input('enter the length of the pipe :');
L0 = L-Le;
a1 = sprintf('kinematic viscosity (m2/s) is %d',v_kin);
a2 = sprintf('Average velocity is %d',V_avg);
disp(a1)
disp(a2)

% Using Hagen poiseuille Equation
l = linspace(0,L0,100);
del_P = ((8*mu*l.*V_avg)./r^2)/rho; %kinematic pressure difference
R = linspace(0,r,100);

u = V_max*(1-(R.^2./r^2)); % velocity along radius
sh_s = mu.*(diff(u)./diff(R));

del_P = flip(del_P);
sh_s1 = flip(sh_s);
u1 = flip(u);

sh_s = [sh_s sh_s1];
u = [u1 u];
R = linspace(-r,r,200);


%---------------------------------------------%

x = input('grading along x:  ');
y = input('grading along y:  ');
z = input('grading along z:  ');

f1 = fopen('blockMeshDict','w');
h1 ='/*--------------------------------*- C++ -*----------------------------------*\';
h2 ='  =========                 |';
h3 ='  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
h4 ='   \\    /   O peration     | Website:  https://openfoam.org';
h5 ='    \\  /    A nd           | Version:  7';
h6 ='     \\/     M anipulation  |';
h7 ='\*---------------------------------------------------------------------------*/';
h8 ='FoamFile';
h9 ='{';
h10 ='    version     2.0;';
h11 ='    format      ascii;';
h12 ='    class       dictionary;';
h13 ='    object      blockMeshDict;';
h14 ='}';
h15 ='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

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\n', h15);

fprintf(f1,'convertToMeters 1;\n\n');

fprintf(f1,'vertices\n');
fprintf(f1,'(\n');
fprintf(f1,'	(0 0 0) 			  // point0\n');
fprintf(f1,'	(%d %d 0)             // point1\n',-r*theta_rad,r*cosd(theta));
fprintf(f1,'	(%d %d 0)             // point2\n',r*theta_rad,r*cosd(theta));
fprintf(f1,'	(0 0 %d) 			  // point3\n',L);
fprintf(f1,'	(%d %d %d)            // point4\n',-r*theta_rad,r*cosd(theta),L);
fprintf(f1,'	(%d %d %d)            // point5\n\n',r*theta_rad,r*cosd(theta),L);
fprintf(f1,');\n\n');

fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,'hex (0 0 2 1 3 3 5 4) (1 35 500) simpleGrading (%d  %d  %d)\n',x,y,z);
fprintf(f1,');\n\n');

fprintf(f1,'edges\n');
fprintf(f1,'(\n');
fprintf(f1,'    arc 1 2 (0 %d 0)\n',r);
fprintf(f1,'    arc 4 5 (0 %d %d)\n',r,L);
fprintf(f1,');\n\n');
fprintf(f1,'boundary\n');
fprintf(f1,'(\n');
fprintf(f1,'    inlet\n');
fprintf(f1,'    {\n');
fprintf(f1,'        type patch;\n');
fprintf(f1,'        faces\n');
fprintf(f1,'        (\n');
fprintf(f1,'            (0 2 1 0)\n');
fprintf(f1,'        );\n');
fprintf(f1,'    }\n\n');
fprintf(f1,'    outlet\n');
fprintf(f1,'    {\n');
fprintf(f1,'        type patch;\n');
fprintf(f1,'        faces\n');
fprintf(f1,'        (\n');
fprintf(f1,'            (3 5 4 3)\n');
fprintf(f1,'        );\n');
fprintf(f1,'    }\n\n');
fprintf(f1,'    front\n');
fprintf(f1,'    {\n');
fprintf(f1,'        type symmetry;\n');
fprintf(f1,'        faces\n');
fprintf(f1,'        (\n');
fprintf(f1,'            (0 3 5 2)\n');
fprintf(f1,'        );\n');
fprintf(f1,'    }\n\n');
fprintf(f1,'    back\n');
fprintf(f1,'    {\n');
fprintf(f1,'        type symmetry;\n');
fprintf(f1,'        faces\n');
fprintf(f1,'        (\n');
fprintf(f1,'            (0 3 4 1)\n');
fprintf(f1,'        );\n');
fprintf(f1,'    }\n\n');
fprintf(f1,'    top\n');
fprintf(f1,'    {\n');
fprintf(f1,'        type wall;\n');
fprintf(f1,'        faces\n');
fprintf(f1,'        (\n');
fprintf(f1,'            (2 5 4 1)\n');
fprintf(f1,'        );\n');
fprintf(f1,'    }\n\n');
fprintf(f1,'    axis\n');
fprintf(f1,'    {\n');
fprintf(f1,'        type empty;\n');
fprintf(f1,'        faces\n');
fprintf(f1,'        (\n');
fprintf(f1,'            (0 3 3 0)\n');
fprintf(f1,'        );\n');
fprintf(f1,'    }\n\n');
fprintf(f1,');\n\n');
fprintf(f1,'mergePatchParis\n');
fprintf(f1,'(\n');
fprintf(f1,');\n\n');
h1 = '// ************************************************************************* //';
fprintf(f1,h1);
fclose(f1);

ControlDict file:

/*--------------------------------*- 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) 			  // point0
	(-1.963495e-03 4.619398e-03 0)             // point1
	(1.963495e-03 4.619398e-03 0)             // point2
	(0 0 2) 			  // point3
	(-1.963495e-03 4.619398e-03 2)            // point4
	(1.963495e-03 4.619398e-03 2)            // point5

);

blocks
(
hex (0 0 2 1 3 3 5 4) (1 35 500) simpleGrading (1  2.000000e-01  1)
);

edges
(
    arc 1 2 (0 5.000000e-03 0)
    arc 4 5 (0 5.000000e-03 2)
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (0 2 1 0)
        );
    }

    outlet
    {
        type patch;
        faces
        (
            (3 5 4 3)
        );
    }

    front
    {
        type symmetry;
        faces
        (
            (0 3 5 2)
        );
    }

    back
    {
        type symmetry;
        faces
        (
            (0 3 4 1)
        );
    }

    top
    {
        type wall;
        faces
        (
            (2 5 4 1)
        );
    }

    axis
    {
        type empty;
        faces
        (
            (0 3 3 0)
        );
    }

);

mergePatchParis
(
);

// ************************************************************************* //

Transport properties:

/*--------------------------------*- 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;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         11.11;

deltaT          1e-3;

writeControl    timeStep;

writeInterval   200;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


// ************************************************************************* //

/*--------------------------------*- 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;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 8.90891e-7;


// ************************************************************************* //

Pressure B.C.

/*--------------------------------*- 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;
    }
    
    axis
    {
	type		empty;
    }
    
    front
    {
        type            symmetry;
    }
    
    back 
    {
    	type 		 symmetry;
    }
    
    top
    {
    	type		 zeroGradient;
    }
}

// ************************************************************************* //

Velocity B.C.:

/*--------------------------------*- 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 1.871e-1);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 1.871e-1);
    }

    outlet
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }
    
    front
    {
        type            symmetry;
    }
    
    back
    {
        type            symmetry;
    }
    
    top
    {
        type            fixedValue;
        value		 uniform (0 0 0);
    }
    
   
}

// ************************************************************************* //

Simulation Output:

Symmetry Boundary condition

angle of 25
Time = 11.11

Courant Number mean: 0.0472556 max: 0.156026
smoothSolver: Solving for Ux, Initial residual = 0.00188682, Final residual = 9.79741e-06, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 4.0495e-09, Final residual = 4.0495e-09, No Iterations 0
smoothSolver: Solving for Uz, Initial residual = 8.95609e-11, Final residual = 8.95609e-11, No Iterations 0
DICPCG: Solving for p, Initial residual = 1.09433e-06, Final residual = 6.69831e-07, No Iterations 1
time step continuity errors : sum local = 1.78547e-11, global = -1.27127e-11, cumulative = -2.62903e-08
DICPCG: Solving for p, Initial residual = 6.70129e-07, Final residual = 6.70129e-07, No Iterations 0
time step continuity errors : sum local = 1.78627e-11, global = -1.27125e-11, cumulative = -2.6303e-08
ExecutionTime = 2208.54 s ClockTime = 2210 s

End

angle of 10

Time = 11.11

Courant Number mean: 0.0472556 max: 0.154523
smoothSolver: Solving for Ux, Initial residual = 0.00463987, Final residual = 5.91007e-07, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 2.51388e-09, Final residual = 2.51388e-09, No Iterations 0
smoothSolver: Solving for Uz, Initial residual = 1.13443e-10, Final residual = 1.13443e-10, No Iterations 0
DICPCG: Solving for p, Initial residual = 7.36023e-07, Final residual = 7.36023e-07, No Iterations 0
time step continuity errors : sum local = 1.88917e-11, global = 9.54027e-12, cumulative = -4.36834e-08
DICPCG: Solving for p, Initial residual = 7.37339e-07, Final residual = 7.37339e-07, No Iterations 0
time step continuity errors : sum local = 1.89254e-11, global = 9.54059e-12, cumulative = -4.36738e-08
ExecutionTime = 2261.23 s ClockTime = 2263 s

End

angle of 45
Time = 11.11

Courant Number mean: 0.0472557 max: 0.160073
smoothSolver: Solving for Ux, Initial residual = 0.000993496, Final residual = 4.85593e-06, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 2.54275e-09, Final residual = 2.54275e-09, No Iterations 0
smoothSolver: Solving for Uz, Initial residual = 9.48883e-11, Final residual = 9.48883e-11, No Iterations 0
DICPCG: Solving for p, Initial residual = 1.25196e-06, Final residual = 3.13538e-07, No Iterations 1
time step continuity errors : sum local = 9.13367e-12, global = -3.50806e-13, cumulative = -1.88217e-08
DICPCG: Solving for p, Initial residual = 3.0977e-07, Final residual = 3.0977e-07, No Iterations 0
time step continuity errors : sum local = 9.02389e-12, global = -3.52739e-13, cumulative = -1.88221e-08
ExecutionTime = 2515.43 s ClockTime = 2517 s

End

------------------------------------------------------------------------------------------------------------------

WEDGE

angle 10

Time = 11.11

Courant Number mean: 0.0519811 max: 0.170183
smoothSolver: Solving for Ux, Initial residual = 0.000632772, Final residual = 4.02303e-06, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 3.64413e-09, Final residual = 3.64413e-09, No Iterations 0
smoothSolver: Solving for Uz, Initial residual = 1.75738e-10, Final residual = 1.75738e-10, No Iterations 0
DICPCG: Solving for p, Initial residual = 1.09866e-06, Final residual = 2.46881e-07, No Iterations 1
time step continuity errors : sum local = 7.61292e-12, global = 6.68762e-12, cumulative = 4.91383e-08
DICPCG: Solving for p, Initial residual = 2.46064e-07, Final residual = 2.46064e-07, No Iterations 0
time step continuity errors : sum local = 7.58771e-12, global = 6.6852e-12, cumulative = 4.9145e-08
ExecutionTime = 1944.7 s ClockTime = 1945 s

End

ANGLE 25

Time = 11.11

Courant Number mean: 0.0519812 max: 0.171706
smoothSolver: Solving for Ux, Initial residual = 0.00131441, Final residual = 6.6321e-06, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 4.16072e-09, Final residual = 4.16072e-09, No Iterations 0
smoothSolver: Solving for Uz, Initial residual = 1.17138e-10, Final residual = 1.17138e-10, No Iterations 0
DICPCG: Solving for p, Initial residual = 1.2358e-06, Final residual = 2.01499e-07, No Iterations 1
time step continuity errors : sum local = 6.41224e-12, global = -5.18205e-13, cumulative = -4.02699e-08
DICPCG: Solving for p, Initial residual = 2.00924e-07, Final residual = 2.00924e-07, No Iterations 0
time step continuity errors : sum local = 6.39395e-12, global = -5.19397e-13, cumulative = -4.02704e-08
ExecutionTime = 1849.06 s ClockTime = 1849 s

End

Angle of 45

Time = 11.11

Courant Number mean: 0.0519813 max: 0.176142
smoothSolver: Solving for Ux, Initial residual = 0.00109939, Final residual = 6.7226e-06, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 3.25315e-09, Final residual = 3.25315e-09, No Iterations 0
smoothSolver: Solving for Uz, Initial residual = 9.4105e-11, Final residual = 9.4105e-11, No Iterations 0
DICPCG: Solving for p, Initial residual = 6.397e-07, Final residual = 6.397e-07, No Iterations 0
time step continuity errors : sum local = 2.21966e-11, global = 1.15224e-11, cumulative = 4.1569e-08
DICPCG: Solving for p, Initial residual = 6.40939e-07, Final residual = 6.40939e-07, No Iterations 0
time step continuity errors : sum local = 2.22396e-11, global = 1.15221e-11, cumulative = 4.15805e-08
ExecutionTime = 2145.86 s ClockTime = 2147 s

End

-----------------------------------------------------------------------------------------------------------------

WEDGE BOUNDARY CONDITION

ANGLE	CLOCK TIME
10	1945 s
25	1849 s
45	2147 s

SYMMETRY BOUNDARY CONDITION

ANGLE	CLOCK TIME
10	2315 s
25	2263 s
45	2517 s

Wedge:

Velocity along length:

We can see Entry length of 1.26m after that the change in velocity is minimum Thus the flow becomes a fully developed flow.

Velocity profile along Radius:

Angle 10 degree:

Angle 25 degree:

Angle 45 degree:

Symmetry:

Angle 10 degree:

Angle 25 degree:

Angle 45 degree:

Vmax:

Error = computation velocity/Analytical velocity

Angle	Wedge b.c.	Error %	symmetry b.c.	Error %
10	0.3729	0.34	0.37291	0.34
25	0.3732	0.26	0.37299	0.32
45	0.3731	0.29	0.39314	0.29

Conclusion:

• Thus Error difference in Vmax for Wedge boundary condition and Symmetry Boundary condition is negligible.
• Also, the time required for Symmetry is more compare to the Wedge Boundary condition.