Week 12 - Symmetry vs Wedge vs HP equation

Objective: To make the simulation work with Symmetry Boundary Condition. Write a Matlab program that takes an angle as input and generates a blockMesh file for the given angle. Angles to test: 10,25,45.

Theory:

SymmetryPlane BC: The symmetry boundary condition defines a mirror face/surface. The symmetry condition should only be used if the physical object or geometry and the expected flow field pattern of the developed solution are mirrored along that surface. This boundary condition helps to reduce the computational domain in size and enables the modeling of a sub-domain of the complete setup.
In general, the symmetry condition applies the following constraints on the flow variables:
The fluxes across a symmetry are zero. The normal components of velocity are set to zero. The normal components of all other variables are set to zero.
The 'symmetry' boundary condition can be applied to both planer or non-planer faces/surfaces on the domain boundaries.

• For this boundary condition no input variables are required and all relevant variable is assigned the symmetry condition.
•  This condition should not be used for axis-symmetric cases

Wedge BC: The wedge boundary condition defines an axis-symmetric boundary condition. The model and flow must be axis-symmetric along a central line such that all physical variables of the flow have the same value and distribution at a radius and angle.
The axis-symmetric 'wedge' boundary is specified by two planes that must be selected on separate sides of the domain on surfaces running along the central line.

•  In case of a 2-dimensional cylinder, the geometry is specified as a wedge of small-angle (< 5 deg) and 1 cell thick.

MATLAB CODE:

clear all
close all
clc

%input parameters

Re= 2100; % Reynold's Number 
D= 0.02; % Diameter of pipe in meters 
R=D/2; % Radius of pipe 
theta = input('Enter value for Wedge angle ='); % Wedge angle
mu=1.002e-3; %Dynamic viscocity of water @ 20 degree Celsius in SI units 
rho=998.21; %Density of pure water highest @ 20 degree Celsius in SI units
nu= mu/rho; %Kinematic viscocity of water @ 20 degree Celsius in SI units 

% Analytical calculation of Velocity,pressure and entrance length 
L_e= 0.06*Re*D; % Entrance length from inlet of pipe 
L_t= L_e + 0.25*L_e; % Assumed Total length of pipe 
v_avg= (Re*nu)/D; % Average velocity of water throughtout characteristic length 
v_max= 2*v_avg; % Average velocity of water at center of pipe characteristic length 
Delta_P= (32*mu*v_avg*L_t/D^2); % Pressure Drop
kin_P= Delta_P/rho; % Kinematic Pressure

% Finding 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_t 0 0]; 
v4= [L_t R*cosd(theta/2) -R*sind(theta/2)]; 
v5= [L_t R*cosd(theta/2) R*sind(theta/2)];

%% Script to form blockMeshDict file. 


f1 = fopen('blockMeshDict.txt','w'); %to create and write blockmesh file

fprintf(f1,'%s','/*---------------------*-C++ -*----------------------*');
fprintf(f1,'\n========             |\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: 7\n');
fprintf(f1,'   \\/    M anipulation | \n');
fprintf(f1,'%s','*-----------------------------------------------------*/');
fprintf(f1,'\nFoamFile\n');

b1=blanks(5); 
b2=blanks(6); 
b3=blanks(7); 
b4=blanks(4); 

fprintf(f1,'{\n'); 
fprintf(f1,'%s','     version',b1,'2.0;'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s','     format',b2,'ascii;');
fprintf(f1,'\n'); 
fprintf(f1,'%s','     class',b3,'dictionary;'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s','     object',b2,'blockMeshDict;'); 
fprintf(f1,'\n'); 
fprintf(f1,'}\n'); 
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * //');

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

% To form vertices 

fprintf(f1,'vertices\n'); 

fprintf(f1,'(\n'); 
fprintf(f1,b4,'%s','('); 
v_0= fprintf(f1,'(%f %f %f)\n',v0(1),v0(2),v0(3)); 
fprintf(f1,b4,'%s','('); 
v_1= fprintf(f1,'(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f1,b4,'%s','('); 
v_2= fprintf(f1,'(%f %f %f)\n',v2(1),v2(2),v2(3)); 
fprintf(f1,b4,'%s','('); 
v_3= fprintf(f1,'(%f %f %f)\n',v3(1),v3(2),v3(3)); 
fprintf(f1,b4,'%s','('); 
v_4= fprintf(f1,'(%f %f %f)\n',v4(1),v4(2),v4(3)); 
fprintf(f1,b4,'%s','('); 
v_5= fprintf(f1,'(%f %f %f)\n',v5(1),v5(2),v5(3)); 
fprintf(f1,'\n);\n'); 

% To create blocks

fprintf(f1,'\nblocks\n');

fprintf(f1,'(\n'); 
fprintf(f1,b4,'%s','('); 
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)'); 

% to create mesh 
m1= input('Enter no. of mesh element along X-axis:');
m2= input('Enter no. of mesh element along Y-axis:'); 
m3= input('Enter no. of mesh element along Z-axis:');

fprintf(f1,blanks(1),'%s','('); 
fprintf(f1,'(%3.0f %1.0f %3.0f)',m1,m2,m3); 
fprintf(f1,'%s',' simpleGrading (1 1 1)'); 
fprintf(f1,'\n);'); 

% To create edges 

fprintf(f1,'\nedges\n'); 
fprintf(f1,'(\n'); 
fprintf(f1,'%s arc 1 2 (0 %f 0)',blanks(1),R); 
fprintf(f1,'\n%s arc 4 5 (%f %f 0)',blanks(1),L_t,R); 
fprintf(f1,'\n);\n\n'); 

% To apply boundary conditions 

% For axis 
fprintf(f1,'boundary\n(\n'); 
fprintf(f1,'%s','  axis'); 
fprintf(f1,'\n  {\n'); 
fprintf(f1,'%s',blanks(10),'type empty;'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'faces'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'('); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(13),'(0 3 3 0)'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),');'); 
fprintf(f1,'\n }\n'); 

% For fixedwalls 
fprintf(f1,'%s','  fixedWalls'); 
fprintf(f1,'\n  {\n'); 
fprintf(f1,'%s',blanks(10),'type wall;'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'faces'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'('); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(13),'(2 5 4 1)'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),');'); 
fprintf(f1,'\n }\n'); 

% For front
fprintf(f1,'%s','  front'); 
fprintf(f1,'\n  {\n'); 
fprintf(f1,'%s',blanks(10),'type symmetry;'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'faces'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'('); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(13),'(0 3 5 2)'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),');'); 
fprintf(f1,'\n  }\n'); 

% For back 
fprintf(f1,'%s','  back'); 
fprintf(f1,'\n  {\n'); 
fprintf(f1,'%s',blanks(10),'type symmetry;'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'faces'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'('); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(13),'(0 1 4 3)'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),');'); 
fprintf(f1,'\n }\n'); 

%For inlet
fprintf(f1,'%s',' inlet');
fprintf(f1,'\n  {\n');
fprintf(f1,'%s',blanks(10),'type patch;'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'faces'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'('); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(13),'(0 0 2 1)'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),');'); 
fprintf(f1,'\n }\n\n'); 

% For outlet 
fprintf(f1,'%s','  outlet'); 
fprintf(f1,'\n  {\n'); 
fprintf(f1,'%s',blanks(10),'type patch;'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'faces'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),'('); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(13),'(3 4 5 3)'); 
fprintf(f1,'\n'); 
fprintf(f1,'%s',blanks(10),');'); 
fprintf(f1,'\n  }\n\n'); 

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

fprintf(f1,'\nmergePatchPairs'); 
fprintf(f1,'\n(\n);\n\n'); 

fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'); 


fclose(f1);

blockMeshDict 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.000000 0.000000 0.000000)
    (0.000000 0.009763 -0.002164)
    (0.000000 0.009763 0.002164)
    (3.150000 0.000000 0.000000)
    (3.150000 0.009763 -0.002164)
    (3.150000 0.009763 0.002164)

);

blocks
(
    hex (0 1 2 0 3 4 5 3) ( 50 1 100) simpleGrading (1 1 1)
);
edges
(
  arc 1 2 (0 0.010000 0)
  arc 4 5 (3.150000 0.010000 0)
);

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

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


);

mergePatchPairs
(
);

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

pressure script

/*--------------------------------*- 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
{
    axis
    {
        type            empty;
    }

    fixedWalls
    {
        type            zeroGradient;
    }

    front
    {
        type            symmetry;
    }
   
    back
    {
        type             symmetry;
    }
 
    inlet
    {
        type            zeroGradient;
    }
   
    outlet
    {
        type             fixedValue;
        value            uniform 0.026661;
    }
}

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

velocity script

/*--------------------------------*- 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.1054 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }

    fixedWalls
    {
        type            noSlip;
    }

    front
    {
        type            symmetry;
    }
    
    back
    {
        type            symmetry;
    }

    inlet
    {
        type            fixedValue;
        value           uniform (0.1054 0 0);
    }

    outlet
    {
         type           zeroGradient;
    }
}

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

transport property script

/*--------------------------------*- 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] 1e-6;


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

controlDict script

/*--------------------------------*- 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         50;

deltaT          0.1;

writeControl    timeStep;

writeInterval   25;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

wedge angle = 10

Plots:

From the plots we see that the flow is fully developed at X = 2.52m.

Similarly the simulation is carried out for other wedge angles and the results are as below.

Results:

symmetry BC

wedge BC

Conclusion:

From the above results we can conclude that symmetry boundary condition is very much accurate than compared to wedge bundary condition as the error percentage obtained in the symmetry BC is very less. Also, symmetry BC allows us to opt for higher wedge angles than lower wedge angles(less than 5 degrees) used in wedge BC since being an axis symmetric BC.