Week 11 - Simulation of Flow through a pipe in OpenFoam
close all
clear all
clc
s = 1; %scaling factor
L = 5; %Length of pipe
r = 0.015; %Radius of pipe
theta = 2; %Angle of wedge wrt axis
rho = 997; %density of working fluid
Re =2100; %Renolds number
dx = 20; %Number of cells in X-direction
dy = 1; %Number of cells in y-direction
dz = 500; %Number of cells in z-direction
Xsg = 0.3; %Simple grading along X-direction
Ysg = 1; %Simple grading along Y-direction
Zsg = 1; %Simple grading along Z-direction
t = 210; %End time
dt = 0.001; %Time step
i =50; %writing interval
mu = 0.997*10^(-3); %dynamic viscosity
FOAMscripts(s,L,r,theta,rho,Re,dx,dy,dz,Xsg,Ysg,Zsg,t,dt,i,mu)
Function code to create all required files:
function FOAMscripts(s,L,r,theta,rho,Re,dx,dy,dz,Xsg,Ysg,Zsg,t,dt,i,mu)
nue = mu/rho; %kinematic viscosity
v = (Re*nue)/(2*r); %Velocity
x1 = r*cosd(theta);
y1 = r*sind(theta);
z1 = 0;
x2 = r*cosd(360-theta);
y2 = r*sind(360-theta);
z2 = L;
%common header lines
line1='/*--------------------------------*-C++ -*--------------------------------*';
line2=' ========= | ';
line3=' / F ield | OpenFOAM: The Open Source CFD Toolbox ';
line4=' / O peration | Website: https://openfoam.org ';
line5=' / A nd | Version: 8 ';
line6=' / M anipulation | ';
line7='*-------------------------------------------------------------------------*/';
line8='FoamFile';
line9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
mkdir 'Flow_through_a_pipe';
mkdir 'Flow_through_a_pipe' 'system';
mkdir 'Flow_through_a_pipe' 'constant';
mkdir 'Flow_through_a_pipe' '0';
%% blockmeshdict
f=fopen('Flow_through_a_pipesystemblockMeshDict','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat tascii;ntclass ttdictionary;ntobject tblockMeshDict;n}n');
%starting line
fprintf(f,'%s',line9);
%defining scaling factor
fprintf(f,'nnconvertToMeters %d;nn',s);
%defining vertices
fprintf(f,'verticesn(nnt(0 0 0)nt(%f %f %f)nt(%f %f %f)nt(0 0 %f)nt(%f %f %f)nt(%f %f %f)n);n',x1,y1,z1,x1,y1,z2,z2,x2,y2,z1,x2,y2,z2);
%defining blocks
fprintf(f,'nblocksn(nthex (0 4 1 0 3 5 2 3) (%d %d %d) simpleGrading (%f %f %f)n);n',dx,dy,dz,Xsg,Ysg,Zsg);
%defining edges
fprintf(f,'nedgesn(ntarc 1 4 (%f %f %f)ntarc 2 5 (%f %f %f)n);n',x1/2,y1/2,z1,x2/2,y2/2,z2);
%defining boundaries
fprintf(f,'nboundaryn(nntaxisnt{ntttype empty;nttfacesntt(nttt(0 3 3 0)ntt);nt}n');
fprintf(f,'tinletnt{ntttype patch;nttfacesntt(nttt(0 4 1 0)ntt);nt}n');
fprintf(f,'toutletnt{ntttype patch;nttfacesntt(nttt(3 5 2 3)ntt);nt}n');
fprintf(f,'tleftwallnt{ntttype wedge;nttfacesntt(nttt(4 0 3 5)ntt);nt}n');
fprintf(f,'trightwallnt{ntttype wedge;nttfacesntt(nttt(0 1 2 3)ntt);nt}n');
fprintf(f,'touterwallnt{ntttype wall;nttfacesntt(nttt(1 4 5 2)ntt);nt}nnn);n');
%defining patch pairs
fprintf(f,'nmergePatchPairsn(n);n');
%ending line
fprintf(f,'n%s',line9);
fclose(f);
%% control dict
f=fopen('Flow_through_a_pipesystemcontrolDict','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttdictionary;ntlocationt"system";ntobject ttcontrolDict;n}n');
%starting line
fprintf(f,'%sn',line9);
fprintf(f,'napplicationticoFoam;nnstartFromtstartTime;nnstartTimet0;nnstopAtttendTime;nnendTimett%f;nndeltaTtt%fs;',t,dt);
fprintf(f,'nnwriteControlttimeStep;nnwriteIntervalt%d;nnpurgeWritet0;nnwriteFormattascii;nnwritePrecisiont6;nnwriteCompressiontoff;nntimeFormattgeneral;nntimePrecisiont6;nnrunTimeModifiablettrue;',i);
%ending line
fprintf(f,'nn%s',line9);
fclose(f);
%% fvSchemes
f=fopen('Flow_through_a_pipesystemfvSchemes','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttdictionary;ntlocationt"system";ntobject ttfvSchemes;n}n');
%starting line
fprintf(f,'%sn',line9);
fprintf(f,'nddtSchemesn{ntdefaultttEuler;n}nngradSchemesn{ntdefaultttGauss linear;ntgrad(p)ttGauss linear;n}nndivSchemesn{ntdefaultttnone;ntdiv(phi,U)tGauss linear;n}nn');
fprintf(f,'laplacianSchemesn{ntdefaultttGauss linear orthogonal;n}nninterpolationSchemesn{ntdefaultttlinear;n}nnsnGradSchemesn{ntdefaultttorthogonal;n}');
%ending line
fprintf(f,'nn%s',line9);
fclose(f);
%% fvSolution
f=fopen('Flow_through_a_pipesystemfvSolution','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttdictionary;ntlocationt"system";ntobject ttfvSolution;n}n');
%starting line
fprintf(f,'%sn',line9);
fprintf(f,'nsolversn{ntpnt{nttsolverttPCG;nttpreconditionertDIC;ntttolerancet1e-06;nttrelToltt0.05;nt}n');
fprintf(f,'ntpFinalnt{ntt$p;nttrelToltt0;nt}nntUnt{nttsolverttsmoothSolver;nttsmoothertsymGaussSeidel;ntttolerancet1e-05;nttrelToltt0;nt}n}n');
fprintf(f,'nPISOn{ntnCorrectorst2;ntnNonOrthogonalCorrectorst0;ntpRefCellt0;ntpRefValuet0;n}');
%ending line
fprintf(f,'nn%s',line9);
fclose(f);
%% transport properties
f=fopen('Flow_through_a_pipeconstanttransportProperties','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttdictionary;ntlocationt"constant";ntobject tttransportProperties;n}n');
%starting line
fprintf(f,'%sn',line9);
%defining kinematic viscosity of water
fprintf(f,'nnutt[0 2 -1 0 0 0 0] %f;n',nue);
%ending line
fprintf(f,'n%s',line9);
fclose(f);
%% pressure
f=fopen('Flow_through_a_pipep','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttvolScalarField;ntobject ttP;n}n');
%starting line
fprintf(f,'%sn',line9);
%defining pressure conditions
fprintf(f,'ndimensionst[0 2 -2 0 0 0 0];nninternalFieldtuniform 0;nnboundaryFieldn{n');
fprintf(f,'tinletnt{ntttypetzeroGradient;nt}ntoutletnt{ntttypetfixedValue;nttvaluetuniform 0;nt}');
fprintf(f,'ntaxisnt{ntttypetempty;nt}ntleftwallnt{ntttypetwedge;nt}ntrightwallnt{ntttypetwedge;nt}ntouterwallnt{ntttypetzeroGradient;nt}nn}');
%ending line
fprintf(f,'n%s',line9);
fclose(f);
%% velocity
f=fopen('Flow_through_a_pipeU','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttvolVectorField;ntobject ttU;n}n');
%starting line
fprintf(f,'%sn',line9);
%defining velocity conditions
fprintf(f,'ndimensionst[0 1 -1 0 0 0 0];nninternalFieldtuniform (0 0 0);nnboundaryFieldn{n');
fprintf(f,'toutletnt{ntttypetzeroGradient;nt}ntinletnt{ntttypetfixedValue;nttvaluetuniform (0 0 %f);nt}',v);
fprintf(f,'ntaxisnt{ntttypetempty;nt}ntleftwallnt{ntttypetwedge;nt}ntrightwallnt{ntttypetwedge;nt}ntouterwallnt{ntttypetnoSlip;nt}nn}');
%ending line
fprintf(f,'n%s',line9);
fclose(f);