Week 11 - Simulation of Flow through a pipe in OpenFoam

re=2100;
R=0.002;
D=2*R;
mu=0.89e-3;
L_e=0.06*re*D;
rho=997;
theta=4;
L=L_e+0.5;
vavg=(re*mu)/(rho*D);
vmax=2*vavg;
r=linspace(0,R,1000);
dp=(8*mu*L*vavg)/(R^2);
pk=(dp)/rho;
for i=1:length(r)
    tau(i)=(2*mu*vmax*r(i))/(R^2);
end
plot(r,tau)
xlabel('radial distance')
ylabel('shear stress')

header1='/*----------------------------------*- C++ -*------------------------------------*\'
header2='  =========                 |                                                                          ' 
header3='   \\      /  F ield                | OpenFOAM: The Open Source CFD Toolbox      '
header4='    \\    /   O peration        | Website:  https://openfoam.org                     '
header5='     \\  /    A nd                 | Version:  6                                                        '
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, 'FoamFile\n{\n')
fprintf(f1, '%12s \t 2.0;\n','version')
fprintf(f1, '%12s \t ascii;\n','format')
fprintf(f1, '%12s \t dictionary;\n','class')
fprintf(f1, '%12s \t blockMeshDict;\n','object')
fprintf(f1, '}\n')
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%s\n\n',header8)


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

fprintf(f1, 'vertices\n')
fprintf(f1, '(\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,0,0)
fprintf(f1, '\t (%d  %d %d)\n',L,R*cosd(theta/2),R*sind(theta/2))
fprintf(f1, '\t (%d  %d %d)\n',L,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1, ');\n\n')
fprintf(f1, 'blocks\n')
fprintf(f1, '(\n')
fprintf(f1, '\t hex (0 3 5 2 0 3 4 1) (%d %d %d)',400,40,1)
fprintf(f1, ' SimpleGrading (1 0.1 1)\n\n')
fprintf(f1, ');\n\n')
fprintf(f1, 'edges\n')
fprintf(f1, '(\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,R,0)
fprintf(f1, ');\n\n')
fprintf(f1, 'boundary\n')
fprintf(f1, '(\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, '\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, '\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, '\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, '\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, '\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, '); \n\n')
fprintf(f1, 'mergePatchPairs\n')
fprintf(f1, '(\n')
fprintf(f1, ');\n')
header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%s\n',header9)
fclose(f1)

velocity contour for the pipe at wedge angle 4 degree:

Grid generated at wedge angle 4 degree:

velocity profile at x=0.1m(yaxis =velocity magnitude, xaxis radial distance)

velocity profile at x=0.2m

velocity profile at x=0.35

velocity profile at x=0.504 m

velocity profile at x=0.6m

so we can see from the above plots that the velocity profile is developed completely at x=0.6m.

shear stress plot(matlab)

shearstress at x=0.6 of fully developed flow from paraview(x axis=radial distance,y axis=shear stress):

Report:

The code for the blockMeshDict was written in matlab and is attached as above.After this the changes where made in the controlDict file so as to limit the courant number and the field values of presssure ,velocity and transport properties were altered accordingly to get a stable solution.The wedge angle is taken as 4 degrees and the simulation is done accordingly.

When considering a flow in pipe there is a velocity profile devleoped due to the shear forces  which are present in pipe and this velocity profile develops through the length of the pipe

this figure shows the devlopement of the velocity profile from the entry to the exit of the pipe. In the plots given above the velocity profile for x=0.504 and x=0.6m is almost the same and the velocity profile remains same afterwards . Hence for reference the plot for velocity at x=0.6 is given. The plot in matlab is made by using the equation:

tau=nu(dv/dy) where nu is viscosity du/dy is the velocity gradient and tau is the shear stress.The enterance length is measured from the starting of the formation of velocity profile to the length when it fully develops and here we have considered laminar flow beacuse re=2100 is taken so L_e=0.06*re*D.

The graphs for shear stress was made in parafoam by the instructions given.