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

Simulation of Hagen Poiseuille Flow through a pipe in OpenFoam using Symmetry boundary condition (Part 2/2)

AAndria Alex Thomas12 Oct 2021
CFDMATLAB

OBJECTIVE

To simulate an axi-symmetric laminar flow by applying a symmetry boundary condition. This is done because simulating the domain as a whole will increase the computational time and might not even complete the simulation if the system resources are scarce. Angles to be considered are 10, 25 and 45 degrees. The results obtained are validated with the Hagen- Poiseuille's equations for the fully developed flow.

PARAMETERS GIVEN AND TAKEN INTO ACCOUNT FOR

Reynolds number: `2100` (laminar flow)

Working fluid: water

Density: `1000` kg/m^3

Dynamic viscosity `mu` : `0.89text(x)10^-3` Pa.s

Kinematic viscosity `nu` : `10^-6` m^2/s

Assume diameter of pipe: `0.03` m

MATLAB CODE

% Simulation of the flow through a pipe: symmetry boundary condition
% ---------------------------------------------------------------

  clear all
  close all
  clc

% Parameters
  Re= 2100;            % Reynolds number
  D= 0.03;             % Diameter of pipe
  r= D/2;              % Radius
  rho= 1000;           % Density
  mu= 0.00089;         % Dynamic viscosity
  nu= mu/rho;          % Kinematic viscosity
  theta= 10;           % symmetry angle in degrees
  L_h= 0.06*Re*D;      % Hydrodynamic entrance length for laminar flow

% Boundary conditions
  L= L_h+0.75;         % Total length of pipe
  V_avg= Re*mu/(rho*D);% Average velocity
  V_max= 2*V_avg;      % Maximum velocity
  dp= 8*mu*L*V_avg/(r^2);% Pressure

% Vertices
  y= r*cosd(theta/2);
  z= r*sind(theta/2);
  v0= [0 0 0];
  v1= [0 y -z];
  v2= [0 y z];
  v3= [L 0 0];
  v4= [L y -z];
  v5= [L y z];

  grading= 0.25;
  b1= 75;
  b2= 1;
  b3= 300;

% Creating bolckMeshDict file
  f1= fopen('blockMeshDict.txt','w');

  fprintf(f1,'%sn','/*--------------------------------*- C++ -*----------------------------------*');
  fprintf(f1,'%sn','  =========                 |');
  fprintf(f1,'%sn','  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
  fprintf(f1,'%sn','   \    /   O peration     | Website:  https://openfoam.org');
  fprintf(f1,'%sn','    \  /    A nd           | Version:  8');
  fprintf(f1,'%sn','     \/     M anipulation  |');
  fprintf(f1,'%snn','*---------------------------------------------------------------------------*/');

  fprintf(f1,'FoamFile n {n');
  fprintf(f1,'%16s %10s n','version','2.0;');
  fprintf(f1,'%16s %12s n','format','ascii;');
  fprintf(f1,'%16s %17s n','class','dictionary;');
  fprintf(f1,'%16s %20s n','object','blockMeshDict;');
  fprintf(f1,' }n');
  fprintf(f1,'%snn','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');

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

  fprintf(f1,'vertices n (n');
  fprintf(f1,'t(%d %d %d)n',v0);
  fprintf(f1,'t(%d %d %d)n',v1);
  fprintf(f1,'t(%d %d %d)n',v2);
  fprintf(f1,'t(%d %d %d)n',v3);
  fprintf(f1,'t(%d %d %d)n',v4);
  fprintf(f1,'t(%d %d %d)n',v5);
  fprintf(f1,' ); nn');

  fprintf(f1,'blocks n (n');
  fprintf(f1,'t %s (%d %d %d) %s (%d 1 1)n','hex (0 1 2 0 3 4 5 3)',b1,b2,b3,'simpleGrading',grading);
  fprintf(f1,' ); nn');

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

  fprintf(f1,'boundary n (n');
  fprintf(f1,'t%5sn','inlet','{');
  fprintf(f1,'t%12s %sn','type','patch;');
  fprintf(f1,'t%13sn %16sn','faces','(');
  fprintf(f1,'tt%14sn %17snn%13snn','(0 2 1 0)',');','}');

  fprintf(f1,'t%5sn','outlet','{');
  fprintf(f1,'t%12s %sn','type','patch;');
  fprintf(f1,'t%13sn %16sn','faces','(');
  fprintf(f1,'tt%14sn %17snn%13snn','(3 4 5 3)',');','}');

  fprintf(f1,'t%5sn','axis','{');
  fprintf(f1,'t%12s %sn','type','empty;');
  fprintf(f1,'t%13sn %16sn','faces','(');
  fprintf(f1,'tt%14sn %17snn%13snn','(0 3 3 0)',');','}');

  fprintf(f1,'t%5sn','walls','{');
  fprintf(f1,'t%12s %sn','type','wall;');
  fprintf(f1,'t%13sn %16sn','faces','(');
  fprintf(f1,'tt%14sn %17snn%13snn','(1 2 5 4)',');','}');

  fprintf(f1,'t%5sn','wedgefront','{');
  fprintf(f1,'t%12s %sn','type','symmetry;');
  fprintf(f1,'t%13sn %16sn','faces','(');
  fprintf(f1,'tt%14sn %17snn%13snn','(0 3 5 2)',');','}');

  fprintf(f1,'t%5sn','wedgeback','{');
  fprintf(f1,'t%12s %sn','type','symmetry;');
  fprintf(f1,'t%13sn %16sn','faces','(');
  fprintf(f1,'tt%14sn %17snn%13snn','(0 1 4 3)',');','}');
  fprintf(f1,' ); nn');

  fprintf(f1,'mergePatchPairs n (n');
  fprintf(f1,' ); nn');

  fprintf(f1,'// ************************************************************************* //');
  fclose(f1);

% Plotting the analytical velocity profile
  R=linspace(0,r,16);
  U=2*V_avg*(1-((R.^2)/(r^2)));
  figure(1)
  plot(R,U,'b')
  hold on
  plot(-R,U,'b')
  hold on
  xlabel('Radius in m')
  ylabel('Velocity')
  title('Analytical Velocity profile at fully developed region','FontSize',9)

blockMeshDict 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; 
          object       blockMeshDict; 
 }
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1; 
 
vertices 
 (
	(0 0 0)
	(0 1.494292e-02 -1.307336e-03)
	(0 1.494292e-02 1.307336e-03)
	(4.530000e+00 0 0)
	(4.530000e+00 1.494292e-02 -1.307336e-03)
	(4.530000e+00 1.494292e-02 1.307336e-03)
 ); 

blocks 
 (
	 hex (0 1 2 0 3 4 5 3) (75 1 300) simpleGrading (2.500000e-01 1 1)
 ); 

edges 
 (
	 arc 1 2 (0 1.500000e-02 0)
	 arc 4 5 (4.530000e+00 1.500000e-02 0)
 ); 

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

            }

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

            }

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

            }

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

            }

	wedgefront
	    {
	        type symmetry;
	        faces
                (
		     (0 3 5 2)
                );

            }


	wedgeback
	    {
	        type symmetry;
	        faces
                (
		     (0 1 4 3)
                );

            }

 ); 

mergePatchPairs 
 (
 ); 

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

RESULTS FOR WEDGE ANGLE 10 DEGREES

  • Mesh generated

  • Velocity distribution from the inlet of the pipe

  • Velocity profiles

  • Velocity profile for fully developed flow

  • Kinematic Pressure

  • Shear stress distribution

RESULTS FOR WEDGE ANGLE 25 DEGREES

  • Mesh generated

  • Velocity distribution from the inlet of the pipe

  • Velocity profiles

  • Velocity profile for fully developed flow

RESULTS FOR WEDGE ANGLE 45 DEGREES

  • Mesh generated

  • Velocity distribution from the inlet of the pipe

  • Velocity profiles

  • Velocity profile for fully developed flow

COMPARISON OF WEDGE& SYMMETRY BOUNDARY CONDITIONS WITH HAGEN POISEUILLE EQUATION

Config/Deg Max.velocity by simulation Hagen Poiseuelle value Error
wed/3 0.123963 0.1246m/s 0.000637
wed/4 0.123966 0.000634
symm/10 0.123981 0.000619
symm/25 0.124102 0.000498
symm/45 0.124261 0.000339

Error is the least for 45 degree configuration. Thus, the more the simulating domain is similar to the actual geometry, the less is the error.

CONCLUSIONS

  • Symmetry boundary conditions is more accurate than wedge boundary conditions
  • Symmetry boundary conditions take more time to execute
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.