Comparison between Boundary Conditions using OpenFOAM
Comparison between Boundary Conditions: Symmetry, Wedge and Hagen Poiseuille's equation
AIM:
To write a code for the simulation of a laminar pipe flow and compare the results for Wedge Boundary Condition, Symmetry Boundary Condition and theoretical Hagen Poiseuille solution.
OBJECTIVES OF THE PROJECT:
To create the BlockMeshDict file for symmetry angles 10,25 and 45 degrees using MATLAB.
To simulate the flow using OpenFoam.
To obtain the velocity profile for each developed flow condition.
GOVERNING EQUATIONS AND THEORY:
The condition at a symmetry boundary are: 1. No flow across boundary, 2. No scalar flux across boundary.
To implement this condition, normal velocities are set to zero at a symmetry boundary and the values of all other properties just outside the solution domain are equated to their values at the nearest node just inside the domain.
The wedge boundary condition is applied to the section of the pipe flow shown above, for wedge angle less than 5 degrees. For higher angles the "Symmetry Boundary Condition" is used. The `text(symmetryplane)` patch type is used for a planar patch which is a symmetry plane, and the `text(symmetry)` patch type is used for any (non-planar) patch type which use the symmetry plane condition.
BODY OF THE CONTENT:
MATLAB CODE TO GENERATE BLOCKMESHDICT FILE:
%Creating a blockmestdict file for Hagen-Poiseuille Law and Symmetry
%Condition Comparison
clear all, close all,clc
%Flow parameters input
Re=2000;
rho=1000;
D=0.02;
R=D/2;
L=0.06*Re*D;
mu=0.00089;
nu=mu/rho;
%Wedge angle for pipe
theta=10;
%Boundary Conditions
v_avg=Re*nu/D;
v_max=2*v_avg;
L=L+0.3;
dp=(8*v_avg*mu*L)/R^2;
%Vertices
v0=[0 0 0];
v1=[0 R*cosd(theta/2) -R*sind(theta/2)];
v2=[0 R*cosd(theta/2) R*sind(theta/2)];
v3=[L 0 0];
v4=[L R*cosd(theta/2) -R*sind(theta/2)];
v5=[L R*cosd(theta/2) R*sind(theta/2)];
%Grading factors
grading=0.2;
b1=20;
b2=3;
b3=500;
%BlockMesh Dict File
f1=fopen('blockMeshDict.txt','w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*n');
fprintf(f1,'| ========= | |n');
fprintf(f1,'| / F ield | OpenFOAM: The Open Source CFD Toolbox |n');
fprintf(f1,'| / O peration | Website: https://openfoam.org |n');
fprintf(f1,'| / A nd | Version: 8 |n');
fprintf(f1,'| / M anipulation | |n');
fprintf(f1,'*---------------------------------------------------------------------------*/n');
fprintf(f1,'FoamFilen{n');
fprintf(f1,'%16s %10sn','version','2.0;');
fprintf(f1,'%16s %12sn','format','ascii;');
fprintf(f1,'%16s %17sn','class','dictionary;');
fprintf(f1,'%16s %20sn','object','blockMeshDict;');
fprintf(f1,'}n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //n');
fprintf(f1,'convertToMeters 1;nn');
fprintf(f1,'verticesn(n');
fprintf(f1,'t(%f %f %f)n',v0);
fprintf(f1,'t(%f %f %f)n',v1);
fprintf(f1,'t(%f %f %f)n',v2);
fprintf(f1,'t(%f %f %f)n',v3);
fprintf(f1,'t(%f %f %f)n',v4);
fprintf(f1,'t(%f %f %f)n',v5);
fprintf(f1,');nn');
fprintf(f1,'blocksn(n');
fprintf(f1,'hex (0 1 2 0 3 4 5 3) (%d %d %d) simpleGrading (%d 1 1)n',b1,b2,b3,grading);
fprintf(f1,');nn');
fprintf(f1,'edgesn(n');
fprintf(f1,'t arc 1 2 (0 %d 0)n',R);
fprintf(f1,'t arc 4 5 (%d %d 0)n',L,R);
fprintf(f1,');nn');
fprintf(f1,'boundaryn(n');
fprintf(f1,'t %5snt{n','axis');
fprintf(f1,'t %13s %s;n','type','empty');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 3 3 0)',');','}');
fprintf(f1,'t %5snt{n','wall');
fprintf(f1,'t %13s %s;n','type','wall');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(2 5 4 1)',');','}');
fprintf(f1,'t %5snt{n','inlet');
fprintf(f1,'t %13s %s;n','type','patch');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 2 1 0)',');','}');
fprintf(f1,'t %5snt{n','outlet');
fprintf(f1,'t %13s %s;n','type','patch');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(3 4 5 3)',');','}');
fprintf(f1,'t %5snt{n','wedgefront');
fprintf(f1,'t %13s %s;n','type','symmetry');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 3 5 2)',');','}');
fprintf(f1,'t %5snt{n','wedgeback');
fprintf(f1,'t %13s %s;n','type','symmetry');
fprintf(f1,'t %14sn %18s;n','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 1 4 3)',');','}');
fprintf(f1,');nn');
fprintf(f1,'mergePatchPairs n (n );n');
fprintf(f1,'// ************************************************************************* //');
fclose(f1);BLOCKMESHDICT FILE CREATED:
/*--------------------------------*- 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.000000 0.000000 0.000000)
(0.000000 0.009962 -0.000872)
(0.000000 0.009962 0.000872)
(2.700000 0.000000 0.000000)
(2.700000 0.009962 -0.000872)
(2.700000 0.009962 0.000872)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (20 3 500) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (2.700000e+00 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgefront
{
type symmetry;
faces
(
(0 3 5 2)
);
}
wedgeback
{
type symmetry;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //CONTROLDICT FILE:
/*--------------------------------*- 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 10;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //INITIAL CONDITIONS (PRESSURE):
/*--------------------------------*- 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
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0171;
}
wall
{
type zeroGradient;
}
wedgefront
{
type symmetry;
}
wedgeback
{
type symmetry;
}
}
// ************************************************************************* //INITIAL CONDITIONS (VELOCITY):
/*--------------------------------*- 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.089 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.089 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
wedgefront
{
type symmetry;
}
wedgeback
{
type symmetry;
}
}
// ************************************************************************* //OUTPUT:
Simulation Result at the output of pipe:
For 10 degree symmetry angle
For 25 degree symmetry angle
For 45 degree symmetry angle
Velocity Profile Variation at different lengths:
MATLAB was used to extract data from csv files and plot output.
CODE:
%Program to extract and plot the velocity profiles at various input lengths
clear all, close all, clc
%The data is extracted from the respective csv files.
T_10_5=readtable('10degree,0.5.csv','Range','A2:A1001');
T_25_5=readtable('25degree,0.5.csv','Range','A2:A1001');
T_45_5=readtable('45degree,0.5.csv','Range','A2:A1001');
%Table format cannot be used for plotting and is converted to array
V_10_5=table2array(T_10_5);
V_25_5=table2array(T_25_5);
V_45_5=table2array(T_45_5);
T_10_25=readtable('10degree,2.5.csv','Range','A2:A1001');
T_25_25=readtable('25degree,2.5.csv','Range','A2:A1001');
T_45_25=readtable('45degree,2.5.csv','Range','A2:A1001');
%Table format cannot be used for plotting and is converted to array
V_10_25=table2array(T_10_25);
V_25_25=table2array(T_25_25);
V_45_25=table2array(T_45_25);
%The radius is defined as the y-axis parameter for the pipe flow
r=linspace(0,0.01,1000);
%r2 is defined to mirror the pipe flow, for the same velocity profile in
%the lower half of the pipe.
r2=linspace(0,-0.01,1000);
figure(1)
subplot(1,2,1)
plot(V_10_5,r,'r')
hold on;
plot(V_10_5,r2,'r','HandleVisibility','off')
hold on;
plot(V_25_5,r,'b')
hold on;
plot(V_25_5,r2,'b','HandleVisibility','off')
hold on;
plot(V_45_5,r,'g')
hold on;
plot(V_45_5,r2,'g','HandleVisibility','off')
xlabel('Axial Direction')
ylabel('Diameter (m)')
title('Velocity Profile at 0.5m from inlet')
subplot(1,2,2)
plot(V_10_25,r,'r')
hold on;
plot(V_10_25,r2,'r','HandleVisibility','off')
hold on;
plot(V_25_25,r,'b')
hold on;
plot(V_25_25,r2,'b','HandleVisibility','off')
hold on;
plot(V_45_25,r,'g')
hold on;
plot(V_45_25,r2,'g','HandleVisibility','off')
legend('10 degree','25 degree','45 degree')
xlabel('Axial Direction')
ylabel('Diameter (m)')
title('Velocity Profile at 2.5m from inlet')OUTPUT:
CONCLUSION:
Velocity profile obtained:
`V_(0.5m)(10^0)=0.1467 text(m/s)` | `V_(2.5m)(10^0)=0.1594 text(m/s)`
`V_(0.5m)(25^0)=0.1468 text(m/s)` | `V_(2.5m)(25^0)=0.1637 text(m/s)`
`V_(0.5m)(45^0)=0.1469 text(m/s)` | `V_(2.5m)(45^0)=0.1661 text(m/s)`
Theoretical Velocity: `V_(max)=2.V_(avg)=2(0.089)=0.178 text(m/s)`
The velocities obtained for the developed flow show some error, as the flow was simulated for 10s, and if continued to higher end time durations can result in a much satisfactory solution.
It can be observed from the velocity profiles that for all the symmetry angles, the profile at 0.5m is same, and shows variation further as the flow gets developed.
The Courant Number for all the simulations was maintained at 0.114, in the stable regime.
It can be observed from the developed flow velocity profile that the lower symmetry angle simulations show some error. The `45^0` symmetry simulation shows best results.
Comparing the results for the wedge boundary condition, at an angle of `2.5^0`, after 10 seconds of simulation the `V_(2.5m)(2.5^0)=0.1695 text(m/s)`.
The flow development is better for the wedge boundary condition, but the accuracy of the flow can be gauged only in the fully developed regime.