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

Simulation of Flow through a pipe in Matlab OpenFoam 2

JJagaveerapandiyan Selvaraj17 Jan 2020
CFD

Matlab code

i.e-wedge 45

clear all;
close all;
clc;

% % Diameter of pipe in inch to meter
% % P=3.14*d;
% % A=(3.14/4)*(d^2);
% % D=4*A/P=d;
d=8.563/1000; %% which is simply ID of 1/8\" pipe
r=d/2;
theta=45; % wedge angle
% % Arc length
% % property of water at 30c
al=2*3.14*r*(theta/360);
rho=995.65;
mu=0.798*10^-3;
nu=mu/rho;
Re=2100; %% laminar flow

L_ent=0.05*Re*d;
L=L_ent+0.1009; %% make L_ent

% % average velocity
v_avg=Re*nu/d;
% % maximum velocity
v_max=2*v_avg;
% % Hydrolic pressur drop
dP=32*mu*L*v_avg/d^2;
% % kinamatic pressure
dP_k=dP/rho; 
dx=L/(15*2000);
% % calculating vertex
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 0 0];
v4=[L r*cosd(theta/2) -r*sind(theta/2)];
v5=[L r*cosd(theta/2) r*sind(theta/2)];

% % creating blockmesh
header1=\"/*--------------------------------*- C++ -*----------------------------------*\\\";
header2=\"  =========                 |                                                  \";
header3=\"  \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox            \";
header4=\"   \\\\    /   O peration     | Website:  https://openfoam.org                   \";
header5=\"    \\\\  /    A nd           | Version:  7                                      \";
header6=\"     \\\\/     M anipulation  |                                                  \";
header7=\"\\*---------------------------------------------------------------------------*/\";
D1=\"    version     2.0;\";
D2=\"    format      ascii;\";
D3=\"    class       dictionary;\";
D4=\"    object      blockMeshDict;\";
D5=\"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\";

c0=\"convertToMeters 1;\";


blocks1=\"    hex (0 1 2 0 3 4 5 3) (15 1 2000) simpleGrading (1 1 1)\";
footer=\"// ************************************************************************* //\";

b4=blanks(4);
b8=blanks(8);
b12=blanks(12);
% % printing file as a file format
fp=fopen(\'blockMeshDict.txt\',\'wt\');
fprintf(fp,\'%s\\n\',header1);
fprintf(fp,\'%s\\n\',header2);
fprintf(fp,\'%s\\n\',header3);
fprintf(fp,\'%s\\n\',header4);
fprintf(fp,\'%s\\n\',header5);
fprintf(fp,\'%s\\n\',header6);
fprintf(fp,\'%s\\n\',header7);

fprintf(fp,\'FoamFile\\n{\\n\');
fprintf(fp,\'%s\\n\',D1);
fprintf(fp,\'%s\\n\',D2);
fprintf(fp,\'%s\\n\',D3);
fprintf(fp,\'%s\\n}\\n\',D4);
fprintf(fp,\'%s\\n\',D5);

fprintf(fp,\'\\n%s\\n\',c0);

% % printing vertices
fprintf(fp,\'\\nvertices\\n(\\n\');

fprintf(fp,\'\\t(%d %d %d)\\n\',v0);
fprintf(fp,\'\\t(%g %0.4f %0.4f)\\n\',v1);
fprintf(fp,\'\\t(%g %0.4f %0.4f)\\n\',v2);
fprintf(fp,\'\\t(%g %g %g)\\n\',v3);
fprintf(fp,\'\\t(%g %0.4f %0.4f)\\n\',v4);
fprintf(fp,\'\\t(%g %0.4f %0.4f)\\n\',v5);

% % printing blocks
fprintf(fp,\');\\n\\nblocks\\n(\\n\')
fprintf(fp,\'%s\\n);\\n\\n\',blocks1)

% % printing edges
fprintf(fp,\'edges\\n(\\n\');
fprintf(fp,\'%sarc 1 2 (0 %0.4f %0.4f)\\n\',b4,v2(3),v2(3));
fprintf(fp,\'%sarc 4 5 (%g %0.4f %0.4f)\\n\',b4,L,v2(3),v2(3));

% % printing boundary
fprintf(fp,\');\\n\\nboundary\\n(\\n\');
fprintf(fp,\'%saxis\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype empty;\\n%sfaces\\n%s(\\n%s(0 3 3 0)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);

fprintf(fp,\'%swall\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype wall;\\n%sfaces\\n%s(\\n%s(2 5 4 1)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);

fprintf(fp,\'%sinlet\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype patch;\\n%sfaces\\n%s(\\n%s(0 0 2 1)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);

fprintf(fp,\'%soutlet\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype patch;\\n%sfaces\\n%s(\\n%s(3 4 5 3)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);

fprintf(fp,\'%sFront\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype wedge;\\n%sfaces\\n%s(\\n%s(0 3 5 2)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);

fprintf(fp,\'%sBack\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype wedge;\\n%sfaces\\n%s(\\n%s(0 1 4 3)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);


fprintf(fp,\');\\n\\nmergePatchPairs\\n(\\n);\\n\\n\');
fprintf(fp,\'%s\',footer);

% % finding dt 
CFL=0.5;
dt=(CFL*dx)/v_avg;

theta varies 4,10,25,45 condition of front & back type varies as a wedge, symmetry to create code for solve for Matlab

output-file

wedge

/*--------------------------------*- 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 0 0)
	(0 0.0040 -0.0016)
	(0 0.0040 0.0016)
	(1.00002 0 0)
	(1.00002 0.0040 -0.0016)
	(1.00002 0.0040 0.0016)
);

blocks
(
    hex (0 1 2 0 3 4 5 3) (15 1 2000) simpleGrading (1 1 1)
);

edges
(
    arc 1 2 (0 0.0016 0.0016)
    arc 4 5 (1.00002 0.0016 0.0016)
);

boundary
(
    axis
    {
        type empty;
        faces
        (
            (0 3 3 0)
        );
    }
    wall
    {
        type wall;
        faces
        (
            (2 5 4 1)
        );
    }
    inlet
    {
        type patch;
        faces
        (
            (0 0 2 1)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (3 4 5 3)
        );
    }
    Front
    {
        type wedge;
        faces
        (
            (0 3 5 2)
        );
    }
    Back
    {
        type wedge;
        faces
        (
            (0 1 4 3)
        );
    }
);

mergePatchPairs
(
);

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

symmetry

/*--------------------------------*- 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 0 0)
	(0 0.0042 -0.0009)
	(0 0.0042 0.0009)
	(1.00002 0 0)
	(1.00002 0.0042 -0.0009)
	(1.00002 0.0042 0.0009)
);

blocks
(
    hex (0 1 2 0 3 4 5 3) (15 1 2000) simpleGrading (1 1 1)
);

edges
(
    arc 1 2 (0 0.004 0)
    arc 4 5 (1.00002 0.004 0)
);

boundary
(
    axis
    {
        type empty;
        faces
        (
            (0 3 3 0)
        );
    }
    wall
    {
        type wall;
        faces
        (
            (2 5 4 1)
        );
    }
    inlet
    {
        type patch;
        faces
        (
            (0 0 2 1)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (3 4 5 3)
        );
    }
    Front
    {
        type symmetry;
        faces
        (
            (0 3 5 2)
        );
    }
    Back
    {
        type symmetry;
        faces
        (
            (0 1 4 3)
        );
    }
);

mergePatchPairs
(
);

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

shear stress comparison

clear all;
close all; 
clc;
mu=0.798*10^-3;
resw=xlsread(\'s0.95.csv\');
vw=resw(:,1);
raw=resw(:,2);
ress=xlsread(\'w0.95.csv\');
vs=(ress(:,1));
ras=(ress(:,2));
thuw=zeros(length(vw),1);
thus=zeros(length(vs),1);
for i=2:length(vw)
thuw(i)=mu.*((vw(i-1)-vw(i))/(raw(i-1)-(raw(i))));
thus(i)=mu.*((vs(i-1)-vs(i))/(ras(i-1)-(ras(i))));
end
plot(abs(thuw),abs(raw),\'color\',\'b\');
hold on
plot(abs(thus),abs(ras),\'color\',\'r\');
legend(\'wedge\',\'symmetry\');
title(\'shear stress Vs spatial Radial at Fully Develope region\');
xlabel(\'shear stress\');
ylabel(\'spatial Radial\');

velocity profile comparison

clear all;
close all;
clc;

% % reading parafoam result in excel
% % wedge
res25w=xlsread(\'w0.25.csv\');
res50w=xlsread(\'w0.5.csv\');
res95w=xlsread(\'w0.95.csv\');
% % symmetry
res25s=xlsread(\'s0.25.csv\');
res50s=xlsread(\'s0.5.csv\');
res95s=xlsread(\'s0.95.csv\');

loglog(res25w(:,2),res25w(:,1),\'color\',\'b\');
hold on;
loglog(res50w(:,2),res50w(:,1),\'color\',\'m\');
hold on;
loglog(res95w(:,2),res95w(:,1),\'color\',\'r\');
hold on
loglog(res25s(:,2),res25s(:,1),\'-.\',\'color\',\'b\',\'LineWidth\',1);
hold on;
loglog(res50s(:,2),res50s(:,1),\'-.\',\'color\',\'m\',\'LineWidth\',1);
hold on;
loglog(res95s(:,2),res95s(:,1),\'-.\',\'color\',\'r\',\'LineWidth\',1);
legend(\'wedge 0.25m\',\'wedge 0.50m\',\'wedge 0.95m\',\'symmetry 0.25m\',\'symmetry 0.5m\',\'symmetry 0.95m\',\'Location\',\'best\');
title(\'spatial radius vs velocity magnitude\');
ylabel(\'velocity magnitude\');
xlabel(\'spatial radius\');
set(gca,\'ytick\',0:0.02:0.35);
xticks = get(gca,\'XTick\');
set(gca,\'XTickLabel\',xticks);

wedge & symmetry 4

side view OpenFoam result

velocity

Pressure

comparison result

Radial shear stress variation

velocity profile comparison

wedge & symmetry 10

side view OpenFoam result

velocity

Pressure

comparison result

Radial shear stress variation

velocity profile comparison

wedge & symmetry 25

side view OpenFoam result

velocity

Pressure

comparison result

Radial shear stress variation

velocity profile comparison

wedge & symmetry 45

side view OpenFoam result

velocity

Pressure

comparison result

Radial shear stress variation

velocity profile comparison

conclusion

by comparing the result with Hagen Poisuelle\'s equation with a different angle in both symmetry & wedge condition

45- almost the same result in profile and converge faster than all less computational expense

10-difference in result it takes computational expense than all simulation

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.