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

Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation

SShruti Jadhav13 Oct 2022
CFDMATLAB
  • Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation.
  • 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
  • Highligh the differences in both boundary conditions in terms of their usage and aplication
  • Angles to test: 10,25,45
  • Objective-
  • To simulate an axis-symmetric laminar flow through the constant cross-sectional of the pipe by applying the symmetry boundary condition.
  • To simulating for the mentioned angles i.e for 10, 25, 45 degrees using both symmetry and wedge boundary conditions and validate them with HP equations.
  • The flow in which flow variables i.e, velocity and pressure , do not vary with the angular coordinate θ is known a axi-symmetric flow.
  • The wedge boundary condition is used in computational fluid dynamics to define an axisymmetric situation
  • The smooth and streamline flows are called laminar flow.

Reynolds number -

where, Re =Reynolds number
= Density of the fluid
= Flow speed
L = length
= Dynamic viscosity of the fluid.
  • The Hagen-Poiseuille equation describes the parabolic velocity profile of frictional, laminar pipe flows of incompressible Newtonian fluids.
  • Flows are basically caused by pressure differences that push fluid from places of higher pressure to places of lower pressure .
  • The length of the hydrodynamic inlet describes the formation of a velocity profile caused by viscous forces propagating from the pipe wall.
Formula- shear stress
where, F = Applied force
A= Cross sectional area
  • For flow through a circular pipe
  • Re is less than equal to 2300 , then the flow is laminar.
  • Re is less than equal to 4000 , then the flow is transition.
  • Re is greater than equal to 4000 , then the flow is turbulent.

  • Boundary Layer Theory:

    Consider a fluid entering a circular pipe with uniform velocity and no slip.

    The fluid velocity in the middle section must increase to keep the mass flow through the pipe constant, resulting in a velocity gradient.

    The hydrodynamic inlet region is the area from the inlet to the point at which the velocity profile is fully developed, and its length is the hydrodynamic inlet length.

  • Hagen- Poiseuille's law - The Hagen- Poiseuille's law gives the pressure drop in an incompressible and newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
  • Assumption- The fluid is incompressible and newtonian.

  • Matlab code
clear all
close all
clc

% Defining Input Values
% Laminar flow
  Re=2100;

% Wedge angle
  theta=10;  

% Diameter of pipe
  D = 0.04; 

% Radius of pipe
  r = D/2

% Hydrodynamic entry length L
 L=0.06*Re*D;

% Assumed characteristic length
 L_c=L+0.3;

% Dynamic viscosity of air
 mu = 1.81e-5

% Density of air
 rho = 1.204

% Kinematic viscosity
  nu=mu/rho 

% calculating average velocity
  v_avg=(mu*Re)/(rho*D)

% Calculating maximum velocity at the centerline of pipe
  v_max=2*v_avg

% Pressure drop for the laminar flow
  delta_p=(32*mu*v_avg*L)/D^2

R=linspace(0,r,400); 
for i=1:length(R)

    tau(i)=2*mu*v_max*(abs(R(i)))/r^2;

end

for j=1: length(R)

    vel_p(j)=2*v_avg*(1-(abs(R(j))^2/r^2));
end

figure(1)
plot(R,tau)
ylabel('Shear stress')
xlabel('Radial length')
 
figure(2)
plot(R,vel_p)
ylabel('Velocity')
xlabel('Radial length')



% Writing the blockMeshDict file

fileid = fopen('blockMeshDict','w');

fprintf(fileid,'This is a test file for blockMesh\n')
fprintf(fileid,'This is my second line\t\tdata\n');
fprintf(fileid,'The Reynolds number is  %f for this simulation',Re);

fprintf(fileid,'\n \t\t(%f\t\t%f\t\t%f)', 0, 0, 0)
fprintf(fileid,'\n \t\t(%f\t\t%f\t\t%f)', r*sind(theta/2), r*cosd(theta/2),0)
fprintf(fileid,'\n \t\t(%f\t\t%f\t\t%f)', -r*sind(theta/2), r*cosd(theta/2),0)
fprintf(fileid,'\n \t\t(%f\t\t%f\t\t%f)', 0, 0, L)
fprintf(fileid,'\n \t\t(%f\t\t%f\t\t%f)', r*sind(theta/2), r*cosd(theta/2),L)
fprintf(fileid,'\n \t\t(%f\t\t%f\t\t%f)', -r*sind(theta/2), r*cosd(theta/2),L)


%blocks
fprintf(fileid,'%s \n','blocks')
fprintf(fileid,'%s \n','(')
fprintf(fileid,'\t hex(0 3 4 1 0 3 5 2) (%d %d %d) simpleGrading (%d %d %d)\n', 200,20,1,1,0.1,1)
fprintf(fileid,'%s \n',');')

%edges
fprintf(fileid,'\n %s \n','edges')
fprintf(fileid,'%s \n','(')
fprintf(fileid,'\t arc 1 2 (%d %d %d) \n',0,r,0)
fprintf(fileid,'\t arc 1 2 (%d %d %d) \n',L,r,0)
fprintf(fileid,'\n %s',');')

%boundary
fprintf(fileid,'\n %s \n','boundary')
fprintf(fileid,'%s \n','(')
fprintf(fileid,'\n %s \n','axis')
fprintf(fileid,'{ \n \t type empty; \n')
fprintf(fileid,'\t faces \n( \n')
fprintf(fileid,'\t (%d %d %d %d) \n',0,3,3,0)
fprintf(fileid,'\t ) ; \n }\n')

%wall
fprintf(fileid,'\n %s \n','wall')
fprintf(fileid,'%s \n','{')
fprintf(fileid,'\t type wall; \n')
fprintf(fileid,'\t faces \n ( \n')
fprintf(fileid,'\t (%d %d %d %d) \n',1,4,5,2)
fprintf(fileid,'\t ) ; \n }\n')

%inlet
fprintf(fileid,'\n %s \n','inlet')
fprintf(fileid,'%s \n','{')

fprintf(fileid,'\t type patch; \n')
fprintf(fileid,'\t faces \n ( \n')
fprintf(fileid,'\t (%d %d %d %d) \n',0,1,2,0)
fprintf(fileid,'\t ) ; \n }\n')

%outlet
fprintf(fileid,'\n %s \n','outlet')
fprintf(fileid,'%s \n','{')
fprintf(fileid,'\t type patch; \n')
fprintf(fileid,'\t faces \n ( \n')
fprintf(fileid,'\t (%d %d %d %d) \n',3,5,4,3)
fprintf(fileid,'\t ); \n }\n')

%front 
fprintf(fileid,'\n\t front \n { \n')
fprintf(fileid,'\t type wedge ; \n')
fprintf(fileid,'\t faces \n ( \n')   
fprintf(fileid,'\t (%d %d %d %d) \n',0,3,4,1)
fprintf(fileid,'\t ) ; \n \t }\n')

%back
fprintf(fileid,'\n\t back \n { \n')
fprintf(fileid,'\t type wedge ; \n')
fprintf(fileid,'\t faces \n ( \n')   
fprintf(fileid,'\t (%d %d %d %d) \n',0,2,5,3)
fprintf(fileid,'\t ) ; \n \t }\n')

fprintf(fileid,');')
fclose(fileid)
  • Stepwise Explaination of the code-

  • Assign initial Hagen-Pousville flow rates values.

  • % Laminar flow

Re=2100;

% Wedge angle
theta=2;

  • Calculation of properties such as dynamic viscosity, average speed, maximum speed, inlet length.

  • Change in pressure value is calculated with the hagen-pousville formula.

  • Linspace command is used to store the array function.

  • Calculate the shear stress and velocity profile using a for loop.

  • Plotting values ​​to obtain an analytical solution.

  • Shear stress vs radial length and velocity vs radial length graph are plotted.

  • Create a BlockMeshDict file using the write function:

  • Assign vertices at all points of the wedge. (v1, v2, v3, v4, v5)

  • Create a controlDict file using the write function

  • Drafting the blockMeshDict file using the write function in MATLAB.

  • The drafting should be done with a similar pattern as the blockMeshDict file used in openFoam.

  • Writing the vertices of the wedge .

  • Assigning edges, walls, borders, entry, exit, front and back of blocks.

  • Create a pressure and velocity file using the write function and make changes accordingly.
  • Using paraview mesh are created and graphs are plotted.

  • Flow simulation procedure:

    A. To select the solver and copy the folder

    1. Open terminal command in ubuntu.

    2. open tutorials, code: tut

    3. Open the incompressible folder and select icofoam, then cavity. (code: icoFoam)

    4. Copy the cavity folder to the startup folder. using code: 'cp -r pipr/$FOAM_RUN'

    5. Type the run command to open the run folder.

    6. Open the pipe. (code: cd pipe)

    ( Also copy and paste the blockMeshDict file created by MATLAB into the system folder)

    7. Open the system file.

    8. Run the blockMeshDict code file. (code: blockMesh)

  • Results

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

convertToMeters 1;

vertices
(
            (0.000000		0.000000		0.000000)
 		(0.000000		0.019988		0.000698)
 		(0.000000		0.019988		-0.000698)
 		(0.000000		0.000000		0.000000)
 		(0.000000		0.019988		0.000698)
 		(0.000000		0.019988	      -0.000698)
)
blocks 
( 
	 hex(0 3 4 1 0 3 5 2) (200 20 1) simpleGrading (1 1 1)
); 

 edges 
( 
	 arc 1 2 (0 2.000000e-02 0) 
	 arc 5 4 (6.000000e-02 2.000000e-02 0) 

 );
 boundary 
( 

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

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

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

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

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


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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         30;

deltaT          0.0005;

writeControl    runTime;

writeInterval   0.1;

purgeWrite      10;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


// ************************************************************************* //
  • Pressure file
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet

    axis
    {
        type            empty;
    }

    {
        type            fixedValue;
        value           uniform 0;
    }

    wall
    {
        type            zeroGradient;
    }
    
    front
    {
        type            symmetry;
    }
   
    Back
    {
        type             symmetry;
    }

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

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

internalField   uniform (0 0 0.020);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0.020);
    }

   outlet
    {

       type             zeroGradient;
    }

    axis
    {
        type            empty;
    }
    
    wall

    {
        type            noSlip;
    }

    front

    {
        type            symmetry;
    }

    Back
    {
        type            symmetry;
    }
}

// ************************************************************************* //
  • Plots for wedge angle 10 degree-

  • Velocity plot

  • Pressure plot

  • Results for wedge angle 25 degree-
  • Velocity plot

  • Shear stress plot

  • Velocity plot

  • Pressure plot

  • Shear stress plot

  • Results for 45 degree wedge angle-
  • Velocity plot

  • Shear stress plot

  • Shear Stress plot

  • Conclusion-
  • We can see that there is not ,much in the maximum velocities.
  • There is not much difference between the analytical and numerical results, and also with increasing wedge angle, the Courant number increases slightly.
  • we can say that the symmetric boundary condition plays an important role in terms of accuracy in the result.
  • The only observed difference is the execution time is slightly higher than the wedge-type. While the angle increases, the error in maximum speed decreases.
  • Comparison between the two Boundary Conditions.
  • We performed simulation for two types of boundary conditions,i.e, wedge boundary condition and symmetry boundary condition.
  • The Courant number increases slightly in the case of a wedge boundary condition and therefore the symmetry boundary condition will be slightly more accurate. The Symmetry boundary condition treats the mesh as 3D and the calculation is performed along all three global axes. During the simulation of the wedge boundary condition, the simulation is not performed over the wedge arc, but while in the case of the symmetry boundary condition, the calculation is performed over the wedge arc; therefore, the simulation time is higher in this case.
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.