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

Objective -

Perform simulation for the below mentioned angles using both symmetry and wedge boundary conditions and validate them with HP equations like you did in week 11 challenge.

Angles to test:

• 10o

• 25o
• 45o

Expected Results

1. Validate Hydro-dynamic length with the numerical result
2. Validate the fully developed flow velocity profile with its analytical profile
3. Validate maximum velocity and pressured drop for fully developed flow
4. Post-process Shear stress and validate wall shear stress for fully developed flow
5. Highligh the differences in both boundary conditions in terms of their usage and aplication

Calculation of endTime -

This an important parameter to follow in the simulation process. So that, Fluid flow compeletes through entire section of the pipe length.

matlab program for wedge 10 degrees :

clear all 
close all 
clc

%input & assumption values

%reynolds number
Re=2100;

%diameter of the pipe
D=0.004;

%entry length of the pipe
L=0.06*Re*D;

%total length of the pipe
L_p=L+0.7;

%radius of the pipe 
r=D/2;

%angle of wedge(angle must be less than 5)
theta=10;

%density
rho=997;

%dynamic viscosity
mu=0.89e-3;

%creating analytical velocity formula
%hagen poisuilleis equation
V_avg=(Re*mu)/(rho*D);
y=linspace(0,r,20);
u=2*V_avg*(1-(y.^2/r^2));

figure(1)
plot(y,u)
ylabel('velocity')
xlabel('radial distance')
legend('radial distance varies')

V_max =2*V_avg;

%calculate shear stress
dy = y(2)-y(1);
tau=(2*mu*V_max.*y)./(r^2);

figure(2)
plot(y,tau)
ylabel('shear stress')
xlabel('radial distance')
legend('shear stress development')

deltap=(32*mu*L_p*V_avg)/(D^2);

header1='/*--------------------------------*- C++ -*----------------------------------*\';
header2='  =========                 |';
header3='  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
header4='   \\    /   O peration     | Website:  https://openfoam.org';
header5='    \\  /    A nd           | Version:  10';
header6='     \\/     M anipulation  |';
header7='\*---------------------------------------------------------------------------*/';

f1=fopen("blockMeshDict","w");
fprintf(f1,'%s\n',header1);
fprintf(f1,'%s\n',header2);
fprintf(f1,'%s\n',header3);
fprintf(f1,'%s\n',header4);
fprintf(f1,'%s\n',header5);
fprintf(f1,'%s\n',header6);
fprintf(f1,'%s\n',header7);

fprintf(f1,'%s\n','FoamFile');
fprintf(f1,'%s\n','{');
fprintf(f1,'\tformat     ascii;\n');
fprintf(f1,'\tclass      dictionary;\n');
fprintf(f1,'\tobject     blockMeshDict;\n');
fprintf(f1,'%s\n','}');
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1,'%s\n\n',header8);

fprintf(f1,'%s\n\n','convertToMeters 1;');
fprintf(f1,'%s\n','vertices;');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t (%d %d %d)\n',0,0,0);
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1,'\t (%d %d %d)\n',L_p,0,0);
fprintf(f1,'\t (%d %d %d)\n',L_p,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1,'\t (%d %d %d)\n',L_p,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1,'%s\n\n',');');

fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (%d %d %d)\n',300,30,1);
fprintf(f1,' SimpleGrading (1 0.1 1)\n');
fprintf(f1,'%s\n\n',');');

fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t 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_p,r,0);
fprintf(f1,'%s\n\n',');');

fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t inlet \n');
fprintf(f1,'\t {\n');
fprintf(f1,'\t\t type patch;\n');
fprintf(f1,'\t\t faces\n');
fprintf(f1,'\t\t (\n');
fprintf(f1,'\t\t\t (0 1 2 0)\n');
fprintf(f1,'\t\t );\n');
fprintf(f1,'\t }\n');

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

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

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

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

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

fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',');');
header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1,'%s\n',header9);

wedge 10 degrees :

matlab program for 25 degrees :

clear all 
close all 
clc

%input & assumption values

%reynolds number
Re=2100;

%diameter of the pipe
D=0.004;

%entry length of the pipe
L=0.06*Re*D;

%total length of the pipe
L_p=L+0.7;

%radius of the pipe 
r=D/2;

%angle of wedge
theta=25;

%density
rho=997;

%dynamic viscosity
mu=0.89e-3;

%creating analytical velocity formula
%hagen poisuilleis equation
V_avg=(Re*mu)/(rho*D);
y=linspace(0,r,20);
u=2*V_avg*(1-(y.^2/r^2));

figure(1)
plot(y,u)
ylabel('velocity')
xlabel('radial distance')
legend('radial distance varies')

V_max =2*V_avg;

%calculate shear stress
dy = y(2)-y(1);
tau=(2*mu*V_max.*y)./(r^2);

figure(2)
plot(y,tau)
ylabel('shear stress')
xlabel('radial distance')
legend('shear stress development')

deltap=(32*mu*L_p*V_avg)/(D^2);

header1='/*--------------------------------*- C++ -*----------------------------------*\';
header2='  =========                 |';
header3='  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
header4='   \\    /   O peration     | Website:  https://openfoam.org';
header5='    \\  /    A nd           | Version:  10';
header6='     \\/     M anipulation  |';
header7='\*---------------------------------------------------------------------------*/';

f1=fopen("blockMeshDict","w");
fprintf(f1,'%s\n',header1);
fprintf(f1,'%s\n',header2);
fprintf(f1,'%s\n',header3);
fprintf(f1,'%s\n',header4);
fprintf(f1,'%s\n',header5);
fprintf(f1,'%s\n',header6);
fprintf(f1,'%s\n',header7);

fprintf(f1,'%s\n','FoamFile');
fprintf(f1,'%s\n','{');
fprintf(f1,'\tformat     ascii;\n');
fprintf(f1,'\tclass      dictionary;\n');
fprintf(f1,'\tobject     blockMeshDict\n');
fprintf(f1,'%s\n','}');
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1,'%s\n\n',header8);

fprintf(f1,'%s\n\n','convertToMeters 1;');
fprintf(f1,'%s\n','vertices;');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t (%d %d %d)\n',0,0,0);
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1,'\t (%d %d %d)\n',L_p,0,0);
fprintf(f1,'\t (%d %d %d)\n',L_p,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1,'\t (%d %d %d)\n',L_p,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1,'%s\n\n',');');

fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (%d %d %d)\n',300,30,1);
fprintf(f1,' SimpleGrading (1 0.1 1)\n');
fprintf(f1,'%s\n\n',');');

fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t 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_p,r,0);
fprintf(f1,'%s\n\n',');');

fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t inlet \n');
fprintf(f1,'\t {\n');
fprintf(f1,'\t\t type patch;\n');
fprintf(f1,'\t\t faces\n');
fprintf(f1,'\t\t (\n');
fprintf(f1,'\t\t\t (0 1 2 0)\n');
fprintf(f1,'\t\t );\n');
fprintf(f1,'\t }\n');

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

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

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

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

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

fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',');');
header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1,'%s\n',header9);

wedge 25 degrees :

matlab code for 45 degrees wedge:

clear all 
close all 
clc

%input & assumption values

%reynolds number
Re=2100;

%diameter of the pipe
D=0.004;

%entry length of the pipe
L=0.06*Re*D;

%total length of the pipe
L_p=L+0.7;

%radius of the pipe 
r=D/2;

%angle of wedge
theta=45;

%density
rho=997;

%dynamic viscosity
mu=0.89e-3;

%creating analytical velocity formula
%hagen poisuilleis equation
V_avg=(Re*mu)/(rho*D);
y=linspace(0,r,20);
u=2*V_avg*(1-(y.^2/r^2));

figure(1)
plot(y,u)
ylabel('velocity')
xlabel('radial distance')
legend('radial distance varies')

V_max =2*V_avg;

%calculate shear stress
dy = y(2)-y(1);
tau=(2*mu*V_max.*y)./(r^2);

figure(2)
plot(y,tau)
ylabel('shear stress')
xlabel('radial distance')
legend('shear stress development')

deltap=(32*mu*(L_p-L)*V_avg)/(D^2);
figure(3)
plot(L_p,deltap,'x')

header1='/*--------------------------------*- C++ -*----------------------------------*\';
header2='  =========                 |';
header3='  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
header4='   \\    /   O peration     | Website:  https://openfoam.org';
header5='    \\  /    A nd           | Version:  10';
header6='     \\/     M anipulation  |';
header7='\*---------------------------------------------------------------------------*/';

f1=fopen("blockMeshDict","w");
fprintf(f1,'%s\n',header1);
fprintf(f1,'%s\n',header2);
fprintf(f1,'%s\n',header3);
fprintf(f1,'%s\n',header4);
fprintf(f1,'%s\n',header5);
fprintf(f1,'%s\n',header6);
fprintf(f1,'%s\n',header7);

fprintf(f1,'%s\n','FoamFile');
fprintf(f1,'%s\n','{');
fprintf(f1,'\tformat     ascii;\n');
fprintf(f1,'\tclass      dictionary;\n');
fprintf(f1,'\tobject     blockMeshDict\n');
fprintf(f1,'%s\n','}');
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1,'%s\n\n',header8);

fprintf(f1,'%s\n\n','convertToMeters 1;');
fprintf(f1,'%s\n','vertices;');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t (%d %d %d)\n',0,0,0);
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1,'\t (%d %d %d)\n',L_p,0,0);
fprintf(f1,'\t (%d %d %d)\n',L_p,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1,'\t (%d %d %d)\n',L_p,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1,'%s\n\n',');');

fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (%d %d %d)\n',300,30,1);
fprintf(f1,' SimpleGrading (1 0.1 1)\n');
fprintf(f1,'%s\n\n',');');

fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t 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_p,r,0);
fprintf(f1,'%s\n\n',');');

fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t inlet \n');
fprintf(f1,'\t {\n');
fprintf(f1,'\t\t type patch;\n');
fprintf(f1,'\t\t faces\n');
fprintf(f1,'\t\t (\n');
fprintf(f1,'\t\t\t (0 1 2 0)\n');
fprintf(f1,'\t\t );\n');
fprintf(f1,'\t }\n');

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

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

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

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

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

fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',');');
header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1,'%s\n',header9);

wedge 45 degrees :

symmetry for 25 deg :

theory & conclusion

1. Tthe velocity and shear force values in all the angles in wedge and symmetry is same we are just try to vary the angles to chowcase how the boundary looks like.
2. whenever the sheaarforce increses at the wall region we will be able to see the velocity profile increases in the middle of the boundary region.
3. we can observe that radial distance increases when the velocity profile becomes decrease at the wall.
4. once the velocity profile reaches its high velocity it means the shear force at the wall is zero.velocity can reach its maximum at hydrodynamically fully developed flow region and also the radial distance also reaches its maximum in the same region.
5. we can see in both the wedge and symmetry conditions are same.the simulation time for symmetry take more time when its compared to wedge condition because the mesh of the symmetry boundary is finer when compared to the wedge boundary.