BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Simulation of flow through a regular pipe in OpenFoam: Comparison of results obtained from Symmetry and Wedge boundary conditions with the Analytical solution (Hagen-Poiseuille Equation)

SShivanshu Kumar21 Sept 2020
CFDMATLAB

NOTE: This project is the continuation of the last project, "Simulation of Laminar Viscous Flow through a regular Pipe using Wedge Boundary conditions in OpenFOAM and post-processed in ParaView". It is highly recommended to read that project first, to have a better understanding of this one!

Problem Specification:

  • Run the simulation for the pipe flow with symmetry boundary conditions and wedge boundary conditions
  • Write a MATLAB program that takes an angle as input and generates a blockMeshDict file for the angles: Angles to test 10, 25, and 45 degrees
  • Compare the results obtained from the symmetry boundary condition and wedge boundary condition with an analytical solution (Hagen-Poiseuille's equation) and discuss your findings

Introduction:

In this project, we're going to run the pipe flow simulation-same as-that we did in the last project, but this time with different angles and boundary conditions. Here we're using symmetry boundary condition and wedge boundary condition for a wedge angle of 10, 25, and 45 degrees. After running the simulation for different angles and for different boundary conditions, we're going to compare and post-process the results in ParaView for different parameters like velocity, pressure, and shear-stress.

For laminar flow, `(Le)/D=0.06 Re`

where `Le=` Hydrodynamic entrance length

`D=` Diameter of pipe

`Re=` Reynolds number

The fully developed pipe flow is generally known as Hagen-Poiseuille flow. The velocity profile of a fully developed laminar horizontal circular pipe flow is given by:

`V=1/(4μ) (∂p)/(∂x)(y^2-R^2 )`

And, the dynamic pressure difference in the pipe is given by:

`∆p=(32μLV)/D^2`

where `∆p=` dynamic pressure (Pa)

`μ=` dynamic viscosity (Pa.s)

`y=` vertical distance from the pipe axis (m)

`R=` radius of the pipe (m)

`D=` diameter of the pipe (m)

`V=` average velocity

Pipe Geometry:

We considered a smooth cylindrical pipe of diameter 8mm and length 1500mm for this case study.

`L=1500 mm or 1.5m`

`D=8mm or 0.008m`

Preliminary Data:

Fluid: Water

Temperature: `298K`

Dynamic Viscosity @ 298K and 1 atm: `μ=0.001003 Pa.s`

Density of water: `ρ=998.2 (kg)/m^3`

Kinematic Viscosity: `ν=1.0048e-06 m^2/s`

Reynold’s number: `Re=2100`

Analytical Solution: (Hagen-Poiseuille)

Average Velocity, `V=(μ*Re)/(ρ*D)=(0.001003*2100)/(998.2*0.008)=0.2637 m/s`

Maximum Velocity, `V_max=2*V_(avg)=2*0.2637=0.5274 m/s`

Hydrodynamic entrance length, `Le=0.06*Re*D=0.06*2100*0.008=1.008 m`

Total length of the pipe, `L=1.5 m`

Pressure Difference for 1.1m to 1.5m, `∆p=(32*μ*L*V)/D^2 =(32*0.001003*0.4*0.2637)/(0.008*0.008)=52.89822Pa`

NOTE: For constant density problems we can simplify the system of equations. The incompressible solvers transform the pressure from static pressure, `p_s` [Pa] to kinematic pressure `p_k`, i.e. `p_k=p_s/ρ` `[m^2/s^2]`. Note that the kinematic pressure has the name `p` in OpenFOAM solvers, the same as the name for the static pressure `p`. To interpret the results from incompressible solvers, the kinematic pressure should be multiplied by the reference density (constant) value.

Therefore, Kinematic Pressure= Static Pressure/Density=52.89822/998.2=0.0529 `m^2/s^2` (approx.)

MATLAB code to generate shear stress plot (analytical solution):

Shear stress, `tau=mu(delu)/(dely)`

clear all
close all
clc

%Inputs
r=0.004; %radial length
y=linspace(0,r,100);%radius array
Y=linspace(0,r,99); %raidus array
u_max=0.5274; %maxm. velocity
mu=0.001003; %dynamic viscosity

for i=1:length(y)
    u(i)=u_max*(1-((y(i)/r)^2));
end

for i=1:length(y)-1
    du=u(i+1)-u(i);
    dy=y(i+1)-y(i);
    
    tau (i)=mu*(du/dy);
end

figure (1)
plot (Y,tau)
xlabel ('Radial Distance (m)')
ylabel ('Shear-stress (analytical)')
set (gca,'YDIR','reverse')
grid on

Pre-Processing:

Geometry and Mesh generation: Flow through a cylinder can be very challenging because of high resource requirements. To minimize the requirement of high computing resources and considering the axisymmetric property of pipe, only a wedge section of the domain was simulated.

A MATLAB code was written which will auto-generate the computational mesh and geometry for any length, radius, and wedge angle.

MATLAB code for generating blockMeshDict file:

This MATLAB code takes an angle as an input and generates a blockMeshDict file for the angles: 10deg, 25deg, and 45deg. Remember that the radius is 0.004m, and the total length is 1.5m. The blockMeshDict file is used to generate the geometry and computational mesh for the simulation using blockMesh command in OpenFoam.

NOTE: The below-mentioned MATLAB program is for symmetry boundary conditions. To change the number of cells along the y-axis and z-axis, we simply edit the ny and nz values accordingly. 

clear all
close all
clc

%Inputs

%Total length (m)
L=1.5; 
%Radius of the pipe (m)
r=0.004; 
%Wedge angle (deg)
theta=input('Enter the wedge angle=');
%No. of cells along x-axis
nx=500;
%no. of cells along y-axis
ny=40;
%no. of cells along z-axis
nz=1;

line_1='/*--------------------------------*- C++ -*----------------------------------*';
line_2='| =========                 |                                                 |';                
line_3='| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |';
line_4='|  \    /   O peration     | Website:  https://openfoam.org                  |';
line_5='|   \  /    A nd           | Version:  8                                     |';
line_6='|    \/     M anipulation  |                                                 |';
line_7='*---------------------------------------------------------------------------*/';


%blockMeshDict
f1=fopen('blockMeshDict','w')
fprintf(f1,'%sn',line_1);
fprintf(f1,'%sn',line_2);
fprintf(f1,'%sn',line_3);
fprintf(f1,'%sn',line_4);
fprintf(f1,'%sn',line_5);
fprintf(f1,'%sn',line_6);
fprintf(f1,'%sn',line_7);

%FoamFile Generation
fprintf(f1,'FoamFilen{n');
fprintf(f1,'%12st 2.0;n','version');
fprintf(f1,'%11st ascii;n','format');
fprintf(f1,'%10st dictionary;n','class');
fprintf(f1,'%11st blockMeshDict;n','object');
fprintf(f1,'}n');
line_8='// ************************************************************************* //';
fprintf(f1,'%snn',line_8);

fprintf(f1,'convertToMeters 1;nn');

%Vertices
fprintf(f1,'verticesn');
fprintf(f1,'(n');
fprintf(f1,'t (%d %d %d)n',0,0,0); %Point 0
fprintf(f1,'t (%d %d %d)n',0,r*cosd(theta/2),r*sind(theta/2)); %Point 1
fprintf(f1,'t (%d %d %d)n',0,r*cosd(theta/2),-r*sind(theta/2)); %Point 2
fprintf(f1,'t (%d %d %d)n',L,0,0); %Point 3
fprintf(f1,'t (%d %d %d)n',L,r*cosd(theta/2),r*sind(theta/2)); %Point 4
fprintf(f1,'t (%d %d %d)n',L,r*cosd(theta/2),-r*sind(theta/2)); %Point 5
fprintf(f1,');nn');

%Blocks
fprintf(f1,'blocksn');
fprintf(f1,'(n');
fprintf(f1,'     hex (0 3 5 2 0 3 4 1) (%d %d %d) ',nx,ny,nz);
fprintf(f1,'simpleGrading (1 0.1 1)nn');
fprintf(f1,');nn');

%Edges
fprintf(f1,'edgesn(nt arc %d %d (%d %d %d)n',1,2,0,r,0);
fprintf(f1,'t arc %d %d (%d %d %d)n',4,5,L,r,0);
fprintf(f1,');nn');

%Boundary
fprintf(f1,'boundaryn(n');

fprintf(f1,'t inletn'); %Inlet
fprintf(f1,'t{n');
fprintf(f1,'tt type patch;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt(n');
fprintf(f1,'ttt (0 1 2 0)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t outletn'); %Outlet
fprintf(f1,'t{n');
fprintf(f1,'tt type patch;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (3 5 4 3)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t Noslipwalln'); %Top-Noslipwall
fprintf(f1,'t{n');
fprintf(f1,'tt type wall;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (1 4 5 2)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t backn'); %Back
fprintf(f1,'t{n');
fprintf(f1,'tt type symmetry;n'); %replace wedge with symmtery for symmetry BC
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (0 2 5 3)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t frontn'); %Front
fprintf(f1,'t{n');
fprintf(f1,'tt type symmetry;n'); %replace wedge with symmtery for symmetry BC
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (0 3 4 1)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t axisn'); %Axis
fprintf(f1,'t{n');
fprintf(f1,'tt type empty;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (0 3 3 0)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');
fprintf(f1,');nn');

%mergePatchPairs
fprintf(f1,'mergePatchPairsn(n');
fprintf(f1,');n');

fprintf(f1,'%sn',line_8);

fclose(f1);

blockMeshDict: Symmetry Boundary conditions; nx=500, ny=40, nz=1; Wedge angle= 10deg

/*--------------------------------*- C++ -*----------------------------------*
| =========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Website:  https://openfoam.org                  |
|   \  /    A nd           | Version:  8                                     |
|    \/     M anipulation  |                                                 |
*---------------------------------------------------------------------------*/
FoamFile
{
     version	 2.0;
     format	 ascii;
     class	 dictionary;
     object	 blockMeshDict;
}
// ************************************************************************* //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 3.984779e-03 3.486230e-04)
	 (0 3.984779e-03 -3.486230e-04)
	 (1.500000e+00 0 0)
	 (1.500000e+00 3.984779e-03 3.486230e-04)
	 (1.500000e+00 3.984779e-03 -3.486230e-04)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 40 1) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 4.000000e-03 0)
	 arc 4 5 (1.500000e+00 4.000000e-03 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type symmetry;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type symmetry;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

blockMeshDict: Symmetry Boundary conditions; nx=500, ny=40, nz=1; Wedge angle= 25deg

/*--------------------------------*- C++ -*----------------------------------*
| =========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Website:  https://openfoam.org                  |
|   \  /    A nd           | Version:  8                                     |
|    \/     M anipulation  |                                                 |
*---------------------------------------------------------------------------*/
FoamFile
{
     version	 2.0;
     format	 ascii;
     class	 dictionary;
     object	 blockMeshDict;
}
// ************************************************************************* //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 0.00390518 0.000865758)
	 (0 0.00390518 -0.000865758)
	 (1.5 0 0)
	 (1.5 0.00390518 0.000865758)
	 (1.5 0.00390518 -0.000865758)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 40 1) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 0.004 0)
	 arc 4 5 (1.5 0.004 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type symmetry;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type symmetry;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

blockMeshDict: Symmetry Boundary conditions; nx=500, ny=40, nz=1; Wedge angle= 45deg

/*--------------------------------*- C++ -*----------------------------------*
| =========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Website:  https://openfoam.org                  |
|   \  /    A nd           | Version:  8                                     |
|    \/     M anipulation  |                                                 |
*---------------------------------------------------------------------------*/
FoamFile
{
     version	 2.0;
     format	 ascii;
     class	 dictionary;
     object	 blockMeshDict;
}
// ************************************************************************* //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 0.00369552 0.00153073)
	 (0 0.00369552 -0.00153073)
	 (1.5 0 0)
	 (1.5 0.00369552 0.00153073)
	 (1.5 0.00369552 -0.00153073)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 40 1) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 0.004 0)
	 arc 4 5 (1.5 0.004 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type symmetry;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type symmetry;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

blockMeshDict: Symmetry Boundary conditions; nx=500, ny=10, nz=3; Wedge angle= 10deg

/*--------------------------------*- C++ -*----------------------------------*
| =========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Website:  https://openfoam.org                  |
|   \  /    A nd           | Version:  8                                     |
|    \/     M anipulation  |                                                 |
*---------------------------------------------------------------------------*/
FoamFile
{
     version	 2.0;
     format	 ascii;
     class	 dictionary;
     object	 blockMeshDict;
}
// ************************************************************************* //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 3.984779e-03 3.486230e-04)
	 (0 3.984779e-03 -3.486230e-04)
	 (1.500000e+00 0 0)
	 (1.500000e+00 3.984779e-03 3.486230e-04)
	 (1.500000e+00 3.984779e-03 -3.486230e-04)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 11 4) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 4.000000e-03 0)
	 arc 4 5 (1.500000e+00 4.000000e-03 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type symmetry;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type symmetry;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

blockMeshDict: Symmetry Boundary conditions; nx=500, ny=15, nz=5; Wedge angle= 25deg

/*--------------------------------*- C++ -*----------------------------------*
| =========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Website:  https://openfoam.org                  |
|   \  /    A nd           | Version:  8                                     |
|    \/     M anipulation  |                                                 |
*---------------------------------------------------------------------------*/
FoamFile
{
     version	 2.0;
     format	 ascii;
     class	 dictionary;
     object	 blockMeshDict;
}
// ************************************************************************* //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 0.00390518 0.000865758)
	 (0 0.00390518 -0.000865758)
	 (1.5 0 0)
	 (1.5 0.00390518 0.000865758)
	 (1.5 0.00390518 -0.000865758)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 15 5) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 0.004 0)
	 arc 4 5 (1.5 0.004 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type symmetry;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type symmetry;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

blockMeshDict: Symmetry Boundary conditions; nx=500, ny=20, nz=7; Wedge angle= 45deg

/*--------------------------------*- C++ -*----------------------------------*
| =========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Website:  https://openfoam.org                  |
|   \  /    A nd           | Version:  8                                     |
|    \/     M anipulation  |                                                 |
*---------------------------------------------------------------------------*/
FoamFile
{
     version	 2.0;
     format	 ascii;
     class	 dictionary;
     object	 blockMeshDict;
}
// ************************************************************************* //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 0.00369552 0.00153073)
	 (0 0.00369552 -0.00153073)
	 (1.5 0 0)
	 (1.5 0.00369552 0.00153073)
	 (1.5 0.00369552 -0.00153073)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 20 7) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 0.004 0)
	 arc 4 5 (1.5 0.004 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type symmetry;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type symmetry;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

NOTE: The below-mentioned MATLAB program is for generating blockMeshDict file-wedge boundary conditions .

clear all
close all
clc

%Inputs

%Total length (m)
L=1.5; 
%Radius of the pipe (m)
r=0.004; 
%Wedge angle (deg)
theta=input('Enter the wedge angle=');
%No. of cells along x-axis
nx=500;
%no. of cells along y-axis
ny=40;
%no. of cells along z-axis
nz=1;

line_1='/*--------------------------------*- C++ -*----------------------------------*';
line_2='| =========                 |                                                 |';                
line_3='| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |';
line_4='|  \    /   O peration     | Website:  https://openfoam.org                  |';
line_5='|   \  /    A nd           | Version:  8                                     |';
line_6='|    \/     M anipulation  |                                                 |';
line_7='*---------------------------------------------------------------------------*/';


%blockMeshDict
f1=fopen('blockMeshDict','w')
fprintf(f1,'%sn',line_1);
fprintf(f1,'%sn',line_2);
fprintf(f1,'%sn',line_3);
fprintf(f1,'%sn',line_4);
fprintf(f1,'%sn',line_5);
fprintf(f1,'%sn',line_6);
fprintf(f1,'%sn',line_7);

%FoamFile Generation
fprintf(f1,'FoamFilen{n');
fprintf(f1,'%12st 2.0;n','version');
fprintf(f1,'%11st ascii;n','format');
fprintf(f1,'%10st dictionary;n','class');
fprintf(f1,'%11st blockMeshDict;n','object');
fprintf(f1,'}n');
line_8='// ************************************************************************* //';
fprintf(f1,'%snn',line_8);

fprintf(f1,'convertToMeters 1;nn');

%Vertices
fprintf(f1,'verticesn');
fprintf(f1,'(n');
fprintf(f1,'t (%d %d %d)n',0,0,0); %Point 0
fprintf(f1,'t (%d %d %d)n',0,r*cosd(theta/2),r*sind(theta/2)); %Point 1
fprintf(f1,'t (%d %d %d)n',0,r*cosd(theta/2),-r*sind(theta/2)); %Point 2
fprintf(f1,'t (%d %d %d)n',L,0,0); %Point 3
fprintf(f1,'t (%d %d %d)n',L,r*cosd(theta/2),r*sind(theta/2)); %Point 4
fprintf(f1,'t (%d %d %d)n',L,r*cosd(theta/2),-r*sind(theta/2)); %Point 5
fprintf(f1,');nn');

%Blocks
fprintf(f1,'blocksn');
fprintf(f1,'(n');
fprintf(f1,'     hex (0 3 5 2 0 3 4 1) (%d %d %d) ',nx,ny,nz);
fprintf(f1,'simpleGrading (1 0.1 1)nn');
fprintf(f1,');nn');

%Edges
fprintf(f1,'edgesn(nt arc %d %d (%d %d %d)n',1,2,0,r,0);
fprintf(f1,'t arc %d %d (%d %d %d)n',4,5,L,r,0);
fprintf(f1,');nn');

%Boundary
fprintf(f1,'boundaryn(n');

fprintf(f1,'t inletn'); %Inlet
fprintf(f1,'t{n');
fprintf(f1,'tt type patch;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt(n');
fprintf(f1,'ttt (0 1 2 0)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t outletn'); %Outlet
fprintf(f1,'t{n');
fprintf(f1,'tt type patch;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (3 5 4 3)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t Noslipwalln'); %Top-Noslipwall
fprintf(f1,'t{n');
fprintf(f1,'tt type wall;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (1 4 5 2)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t backn'); %Back
fprintf(f1,'t{n');
fprintf(f1,'tt type wedge;n'); %replace symmetry with wedge for wedge BC
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (0 2 5 3)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t frontn'); %Front
fprintf(f1,'t{n');
fprintf(f1,'tt type wedge;n'); %replace symmetry with wedge for wedge BC
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (0 3 4 1)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t axisn'); %Axis
fprintf(f1,'t{n');
fprintf(f1,'tt type empty;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (0 3 3 0)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');
fprintf(f1,');nn');

%mergePatchPairs
fprintf(f1,'mergePatchPairsn(n');
fprintf(f1,');n');

fprintf(f1,'%sn',line_8);

fclose(f1);

blockMeshDict: Wedge Boundary conditions; nx=500, ny=40, nz=1; Wedge angle= 10deg

/*--------------------------------*- C++ -*----------------------------------*
| =========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Website:  https://openfoam.org                  |
|   \  /    A nd           | Version:  8                                     |
|    \/     M anipulation  |                                                 |
*---------------------------------------------------------------------------*/
FoamFile
{
     version	 2.0;
     format	 ascii;
     class	 dictionary;
     object	 blockMeshDict;
}
// ************************************************************************* //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 3.984779e-03 3.486230e-04)
	 (0 3.984779e-03 -3.486230e-04)
	 (1.500000e+00 0 0)
	 (1.500000e+00 3.984779e-03 3.486230e-04)
	 (1.500000e+00 3.984779e-03 -3.486230e-04)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 40 1) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 4.000000e-03 0)
	 arc 4 5 (1.500000e+00 4.000000e-03 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type wedge;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type wedge;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

blockMeshDict: Wedge Boundary conditions; nx=500, ny=40, nz=1; Wedge angle= 25deg

/*--------------------------------*- C++ -*----------------------------------*
| =========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Website:  https://openfoam.org                  |
|   \  /    A nd           | Version:  8                                     |
|    \/     M anipulation  |                                                 |
*---------------------------------------------------------------------------*/
FoamFile
{
     version	 2.0;
     format	 ascii;
     class	 dictionary;
     object	 blockMeshDict;
}
// ************************************************************************* //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 3.905184e-03 8.657585e-04)
	 (0 3.905184e-03 -8.657585e-04)
	 (1.500000e+00 0 0)
	 (1.500000e+00 3.905184e-03 8.657585e-04)
	 (1.500000e+00 3.905184e-03 -8.657585e-04)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 40 1) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 4.000000e-03 0)
	 arc 4 5 (1.500000e+00 4.000000e-03 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type wedge;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type wedge;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

blockMeshDict: Wedge Boundary conditions; nx=500, ny=40, nz=1; Wedge angle= 45deg

/*--------------------------------*- C++ -*----------------------------------*
| =========                 |                                                 |
| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \    /   O peration     | Website:  https://openfoam.org                  |
|   \  /    A nd           | Version:  8                                     |
|    \/     M anipulation  |                                                 |
*---------------------------------------------------------------------------*/
FoamFile
{
     version	 2.0;
     format	 ascii;
     class	 dictionary;
     object	 blockMeshDict;
}
// ************************************************************************* //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 3.695518e-03 1.530734e-03)
	 (0 3.695518e-03 -1.530734e-03)
	 (1.500000e+00 0 0)
	 (1.500000e+00 3.695518e-03 1.530734e-03)
	 (1.500000e+00 3.695518e-03 -1.530734e-03)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 40 1) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 4.000000e-03 0)
	 arc 4 5 (1.500000e+00 4.000000e-03 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type wedge;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type wedge;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

Initial and Boundary Conditions: Symmetry boundary condition

Velocity, U

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  8
     \/     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.2637 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.2637 0 0);
    }
    
    outlet
    {
        type            zeroGradient;
    }

    Noslipwall
    {
        type            noSlip;
    }

    back
    {
        type            symmetry;
    }
    
    front
    {
        type            symmetry;
    }
    
    axis
    {
        type            empty;
    }
}

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

Pressure, p

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  8
     \/     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;
    }
    
    Noslipwall
    {
        type            zeroGradient;
    }

    back
    {
        type            symmetry;
    }
    
    front
    {
        type            symmetry;
    }
    
    axis
    {
        type            empty;
    }
}

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

Initial and Boundary Conditions: Wedge boundary condition

Velocity, U

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  8
     \/     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.2637 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.2637 0 0);
    }
    
    outlet
    {
        type            zeroGradient;
    }

    Noslipwall
    {
        type            noSlip;
    }

    back
    {
        type            wedge;
    }
    
    front
    {
        type            wedge;
    }
    
    axis
    {
        type            empty;
    }
}

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

Pressure, p

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  8
     \/     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;
    }
    
    Noslipwall
    {
        type            zeroGradient;
    }

    back
    {
        type            wedge;
    }
    
    front
    {
        type            wedge;
    }
    
    axis
    {
        type            empty;
    }
}

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

NOTE: The controlDict, fvSchemes, fvSolution, and transportProperties files are going to be same for each simulation.

controlDict

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  8
     /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         20;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

fvSchemes

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  8
     /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}


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

fvSolution

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  8
     /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}


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

transportProperties: (Kinematic Viscosity=1.0048e-06)

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  8
     /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 1.0048e-06;


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

Post-Processing:

NOTE: Remember that the nx, ny, and nz denotes the number of cells along x-axis, y-axis, and z-axis respectively.

Symmetry Boundary Condition: nx=500, ny=40, nz=1; wedge angle= 10deg

Velocity Contour: Hydrodynamic entrace region (near inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.4m

Pressure Contour:

Pressure-drop Plot:

x=1.1m to x=1.5m (Fully developed flow region)

Shear Stress (tau) Plot: x=1.4m

Symmetry Boundary Condition: nx=500, ny=40, nz=1; wedge angle= 25deg

Velocity Contour: Hydrodynamic entrace region (near inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.4m

Pressure Contour:

Pressure-drop Plot:

x=1.1m to x=1.5m (Fully developed flow region)

Shear Stress (tau) Plot:

Symmetry Boundary Condition: nx=500, ny=40, nz=1; wedge angle= 45deg

Velocity Contour: Hydrodynamic entrace region (near inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.4m

Pressure Contour:

Pressure-drop Plot:

x=1.1m to x=1.5m (Fully developed flow region)

Shear Stress (tau) Plot: x=1.4m

Symmetry Boundary Condition: nx=500, ny=10, nz=3; wedge angle= 10deg

Velocity Contour: Hydrodynamic entrace region (near inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.4m

Pressure Contour:

Pressure-drop Plot:

x=1.1m to x=1.5m (Fully developed flow region)

Shear Stress (tau) Plot: x=1.4m

Symmetry Boundary Condition: nx=500, ny=15, nz=5; wedge angle= 25deg

Velocity Contour: Hydrodynamic entrace region (near inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.4m

Pressure Contour:

Pressure-drop Plot:

x=1.1m to x=1.5m (Fully developed flow region)

Shear Stress (tau) Plot: x=1.4m

Symmetry Boundary Condition: nx=500, ny=20, nz=7; wedge angle= 45deg

Velocity Contour: Hydrodynamic entrace region (near inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.4m

Pressure Contour:

Pressure-drop Plot:

x=1.1m to x=1.5m (Fully developed flow region)

Shear Stress (tau) Plot: x=1.4m

Wedge Boundary Condition: nx=500, ny=40, nz=1; wedge angle= 10deg

Velocity Contour: Hydrodynamic entrace region (near inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.4m

Pressure Contour:

Pressure-drop Plot:

x=1.1m to x=1.5m (Fully developed flow region)

Shear Stress (tau) Plot: x=1.4m

Wedge Boundary Condition: nx=500, ny=40, nz=1; wedge angle= 25deg

Velocity Contour: Hydrodynamic entrace region (near inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.4m

Pressure Contour:

Pressure-drop Plot:

x=1.1m to x=1.5m (Fully developed flow region)

Shear Stress (tau) Plot: x=1.4m

Wedge Boundary Condition: nx=500, ny=40, nz=1; wedge angle= 45deg

Velocity Contour: Hydrodynamic entrace region (near inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.4m

Pressure Contour:

Pressure-drop Plot:

x=1.1m to x=1.5m (Fully developed flow region)

Shear Stress (tau) Plot: x=1.4m

CONCLUSION:

  • For the same parameters, the results are the same for both wedge and symmetry boundary conditions.
  • Overall accuracy seems to decrease with increase in wedge angle.
  • The error in velocity appears to increase as the error in pressure decreases.
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.