Week 11 - Simulation of Flow through a pipe in OpenFoam
Hagen Poiseuille Flow
It is of utmost importance to develop formula for pressure difference for laminar flow in the pipe.
Reynolds number is the number which determines the characteristics of the fluid flow in a pipe and hence pressure difference formula has to be derived such that it is a function of Reynolds number.
This theoretical establishment is important to verify theory with the results of the simulation which will be done in Openfoam.
Since it is assumed that the flow is incompressible flow, then the flow will be governed by an incompressible continuity equation and incompressible momentum equation.
Both the governing equation will be in cyllindrical coordinates, since the flow is in a pipe which is an axisymmetrix geometry.
Now the respective notations for directions for cyllindrical coordinates must be understood.
x direction means in longitudinal direction of pipe i.e length wise.
r direction means in transverse direction to the length along the radius of the pipe.
`Theta` direction means in at transverse direction to radius and length of the pipe along the curvature of the pipe.
Continuity Equation:-
`(delrho)/(delt)+1/r*(del(rrhoV_r))/(delr)+1/r*(del(rhoV_theta))/(deltheta)+(del(rhoV_x))/(delx)=0`
Since this equation must describe a steady state of fluid flow
`(delrho)/(delt)=0`
Now since flow is assumed to be steady state, the fluid will not have flow in theta direction and gradient of velocity along x direction will also be zero
`v_theta=0`
`(delV_x)/(delx)=0`...since there is no acceleration of fluid in steady state in x-direction for the gradient of `V_x` to exist.
Simplifying the equation further, we get
`1/r*(del(rV_r))/(delr)=0`
`rho`is taken out and ignored since it is a incompressible flow
Now looking at the equation above and integrated, it is found that
`r*V_r=`constant
But it is to be understood that Velocity at the walls of the pipe does not have velocity in the r direction i.e `V_r=0`at r=R
Therefore putting above conclusion in the above equation we get,
`R*0=`constant
This implies that constant itself is zero. This inturn implies that `V_r`is zero everywhere.
So it is found that `V_r`and `V_theta`is zero everywhere.
Momentum Equation
In x-direction,
`rho((delV_x)/(delt)+V_r*(delV_x)/(delr)+V_theta/r*(delV_x)/(delr)+V_x*(delV_x)/(delr)=-(delp)/(delx)+mu(1/r*del/(delr)(rdelV_x)/(delr)+(1/r^2)*(del^2V_x)/(deltheta^2)+(del^2V_x)/(delx^2))`
Now applying following constrants to the above equation and their respective reason:-
Now since `V_r`and `V_theta`are zero everywhere, the only component existing for fluid velocity is `V_x`
But
`(delV_x)/(dt)=0`...steady state is assumed
`(delV_x)/(delx)=0`...no acceleration of velocity in steady state and hence no steady state. Therefore `(del^2V_x)/(delx^2)=0`
`(delV_x)/(deltheta)=0`...since flow is axisymmetric.
And therefore`(del^2V_x)/(deltheta^2)=0`
By substituting above constraints to momentum equation and simplifying we get,
`(delp)/(delx)=mu/rdel/(delr)(rdelV_x)/(delr)`
Now when pressure difference is applied across the length of the pipe, the pressure gradient becomes constant at each point in the pipe at steady state.
So lets put
`(delp)/(delx)=k=`constant
Now substituting in the derived equation from momentum equation, we get
`del/(delr)(rdelV_x)/(delr)=(rk)/mu`
Now to perform integration, the above partial differential equation must be converted to normal differential equation. now since `V_x`is not a function of theta and x for the reasons stated in the derivation above, this implies that it is a function of r only.
Therefore above equation becomes
`d/(dr)(rdV_x)/(dr)=(rk)/mu`
Now integrating both sides with respect to r, we get
`(rdelV_x)/(delr)=(r^2k)/(2*mu)+C_1`where C_1 is a constant of integration.
Now simplifying above equation, we get
`(delV_x)/(delr)=(rk)/(2*mu)+C_1/r`
Further integrating above equation
`V_x=(r^2k)/(4mu)+C_1*logr+C_2`where C_2 is second constant of integration
Now for fluid, the formula for shear stress in the x-direction with plane of normal in r-direction:-
`tau_(xr)=mu((delV_r)/(delx)+(delV_x)/(delr))`
but since `V_r=0`everywhere
`tau_(xr)=mu(delV_x)/(delr)`
substituting (delV_x)/(delr) into above equation
`tau_(xr)=mu((rk)/(2*mu)+C_1/r)`
now it must be understood when fluid flows over a surface, the velocity of the fluid is zero at the walls and shear stress due to transfer of momentum of stationary wall is zero at points where there is bulk flow of fluid or in other words at distance where velocity is supposed to be maximum and steady. This distance is R (radius of the pipe) from the wall. By the variable r, `tau_(xr)=0`at r=0 where `V_x`is maximum and momentum transfer is negligible.
Putting r=0 in the equation for `tau_(xr)` where the `tau_(xr)=0`
`0=mu(0.k)/(2*mu)+C_1/0`
The presence of term with C_1 makes the entire equation indeterminate, therefore C_1=0 to satisfy the equation.
But at r=R, `V_x=0`therefore substituting in the derived equation
`0=(R^2k)/(4mu)+0*logr+C_2`
therefore
`C_2=-(R^2k)/(4mu)`
therefore the equation for `V_x`becomes
`V_x=((r^2-R^2)k)/(4mu)`
Now since `V_x`is maximum at r=0
`V_(x max)=(-R^2k)/(4mu)`
the negative sign would give the impression that the velocity is minimum rather than maximum. But it must be understood that k which pressure gradient is itself negative, since pressure at the entrance is higher than that at exit.
`k=(dp)/(dx)=-(dp)/L`
therefore the `V_x`becomes
`V_(x max)=(R^2*dp)/(4muL)`
Now to calculate average velocity across the cross-section of pipe, formula for volumetric flow rate has to be developed.
`q=int_0^RVx.dA`
now as it is known `A(r)=pir^2`
therefore
`dA=2pir*dr`
substituting in the equation
`q=int_0^RVx*2pir*dr`
Now simplifying this
`q=2pi*int_0^R((r^2-R^2)k)/(4mu)*r*dr`
Now after integration, the equation becomes
`q=(2pik)/(4mu)*[(r^4)/(4)-(R^2*r^2)/(2)]_0^R`
So simplifying the equation,
`q=-(R^4pik)/(8mu)`
But formula for q with `V_(avg)`
`q=V_(avg)*pi*R^2`
Equating both the equations we get
`V_(avg)*pi*R^2=-(R^4pik)/(8mu)`
So the formula for V_avg becomes
`V_(avg)=-(R^2k)/(8mu)`
and since `k=-(dp)/L`
`V_(avg)=(R^2*dp)/(8muL)`
Now comparing `V_(avg)`with `V_(max)`
`V_(avg)=V_(x max)/2`
Now since formula for `V_(avg)`has been developed, it can be manipulated to get formula for dp as a function `V_(avg)`
`dp=(8*V_(avg)*mu*L)/R^2`
but formula of Reynolds number is
`Re=(rho*V_(avg)*(2R))/mu`
therefore manipulating the formula to get `V_(avg)`as a function of Re
`V_(avg)=(mu*Re)/(2rhoR)`
substituting this formula in the formula for dp we get
`dp=(4mu^2LRe)/(rho*R^3)`
And since kinematic pressure difference is required, therefore
`dp_k=(dp)/(rho)`
`dp_k=(4mu^2LRe)/(rho^2*R^3)`
And the `V_(x max)`formula as a function of Re
`V_(x max)=(mu*Re)/(rhoR)`
now further simplifying
`V_(x max)=(nu*Re)/(R)`
where `nu=mu/rho`= kinematic viscosity
Now when a flow starts from a reservoir body into the pipe, it happens untill some distance that the velocity profile develops. That distance is called Hydrodynamic Length
This image has been taken from Frank M White book for Fluid Dynamics
The empirical formula to calculate Hydrodynamic length is
`(Le)/D=0.06*Re`for laminar flow
Re=Reynold Number, D=diameter of the pipe and Le=Hydrodynamic length
Now since theoretical understanding has been established, specifications for the simulations has to be defined.
i)Reynolds Number=Re=1000
ii) Radius of the cross-section(R)=0.05 metre
iii)Dynamic Viscosity(`mu`)=1.002*10^-3 Pascal Second
i.v) Density(`rho`)= 998 `(kg)/m^3`
v) Kinematic Viscosity(`nu`)= `mu/rho`=1.004e-06
vi) Hydrodynamic Length= 0.06*Re*D=0.06*Re*2*r=6m
vii)Kinematic Pressure difference =`(4mu^2LeRe)/(rho^2*R^3)`=1.9354e-04 `m^2/(s^2)`
Now since pressure gradience across the length of the pipe must be negative, this dp value will be put as input at inlet and zero value will be put at outlet
viii)The ange of the wedge selected is 4 degree.
Expected output:-
xi) Maximum velocity=`(nu*Re)/R`=0.020080 m/s
x) Maximum Shear stress=`tau_(xr max)=(R*dp)/(2*Le)`=0.00080481 `N/m^2`
dx=0.05m
dt=0.01s
xi) Maximum predicted CFL= Vmax*dt/dx=0.0040160
Now for other part of challenge where a matlab code has been written to write a file for blockmesh for wedge if theta is changed to some value. Basically automation is desired in of writting a code fpr BlockMeshDict file. this challenge is extended to write file for controlDict, pressure file and transportProperties file. The code will automate writing process as and when the value of required properties are provided.
Now a code for functions will be written for different files
1)BlockMesh
function []=blockMesh_wedge_function(r,L,theta,text_address,no_x,no_y,no_z,x_grading,y_grading,z_grading)
h=r*cosd(theta);
b=r*sind(theta);
f1=fopen(text_address,'w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n')
fprintf(f1,' ========= |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox')
fprintf(f1,'\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org')
fprintf(f1,'\n')
fprintf(f1,' \\ / A nd | Version: 7')
fprintf(f1,'\n')
fprintf(f1,' \\/ M anipulation |')
fprintf(f1,'\n')
fprintf(f1,'/*---------------------------------------------------------------------------*/')
fprintf(f1,'\n')
fprintf(f1,'FoamFile')
fprintf(f1,'\n')
fprintf(f1,'{')
fprintf(f1,'\n')
fprintf(f1,' version 2.0;')
fprintf(f1,'\n')
fprintf(f1,' format ascii;')
fprintf(f1,'\n')
fprintf(f1,' class dictionary;')
fprintf(f1,'\n')
fprintf(f1,' object blockMeshDict;')
fprintf(f1,'\n')
fprintf(f1,'}')
fprintf(f1,'\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'convertToMeters 1;')
fprintf(f1,'\n')
fprintf(f1,'\n')
%Defining of vertices
fprintf(f1,'vertices')
fprintf(f1,'\n')
fprintf(f1,'(')
fprintf(f1,'\n')
fprintf(f1,' (0 0 0)')
fprintf(f1,'\n')
fprintf(f1,' (0 %10.10f %10.10f)',h,-b)
fprintf(f1,'\n')
fprintf(f1,' (0 %10.10f %10.10f)',h,b)
fprintf(f1,'\n')
fprintf(f1,' (%10.3f 0 0)',L)
fprintf(f1,'\n')
fprintf(f1,' (%10.3f %10.10f %10.10f)',L,h,-b)
fprintf(f1,'\n')
fprintf(f1,' (%10.3f %10.10f %10.10f)',L,h,b)
fprintf(f1,'\n')
fprintf(f1,');')
%End of Defining vertices
%Defining blocks
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'blocks')
fprintf(f1,'\n')
fprintf(f1,'(')
fprintf(f1,'\n')
fprintf(f1,' hex (0 3 4 1 0 3 5 2) (%.0f %.0f %.0f) simpleGrading (%f %f %f)',no_x,no_y,no_z,x_grading,y_grading,z_grading)
fprintf(f1,'\n')
fprintf(f1,');')
%Definition of blocks ends
%defining of edges starts
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'edges')
fprintf(f1,'\n')
fprintf(f1,'(')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,');')
%defining of edges ends
%defining of Boundary starts
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'boundary')
fprintf(f1,'\n')
fprintf(f1,'(')
fprintf(f1,'\n')
fprintf(f1,' inletWall')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type patch;')
fprintf(f1,'\n')
fprintf(f1,' faces')
fprintf(f1,'\n')
fprintf(f1,' (')
fprintf(f1,'\n')
fprintf(f1,' (0 0 2 1)')
fprintf(f1,'\n')
fprintf(f1,' );')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' outletWall')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type patch;')
fprintf(f1,'\n')
fprintf(f1,' faces')
fprintf(f1,'\n')
fprintf(f1,' (')
fprintf(f1,'\n')
fprintf(f1,' (3 4 5 3)')
fprintf(f1,'\n')
fprintf(f1,' );')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' fixedWalls')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type wall;')
fprintf(f1,'\n')
fprintf(f1,' faces')
fprintf(f1,'\n')
fprintf(f1,' (')
fprintf(f1,'\n')
fprintf(f1,' (4 1 2 5)')
fprintf(f1,'\n')
fprintf(f1,' );')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' front')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type wedge;')
fprintf(f1,'\n')
fprintf(f1,' faces')
fprintf(f1,'\n')
fprintf(f1,' (')
fprintf(f1,'\n')
fprintf(f1,' (0 3 5 2)')
fprintf(f1,'\n')
fprintf(f1,' );')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' Back')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type wedge;')
fprintf(f1,'\n')
fprintf(f1,' faces')
fprintf(f1,'\n')
fprintf(f1,' (')
fprintf(f1,'\n')
fprintf(f1,' (0 1 4 3)')
fprintf(f1,'\n')
fprintf(f1,' );')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' axis')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type empty;')
fprintf(f1,'\n')
fprintf(f1,' faces')
fprintf(f1,'\n')
fprintf(f1,' (')
fprintf(f1,'\n')
fprintf(f1,' (0 3 3 0)')
fprintf(f1,'\n')
fprintf(f1,' );')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,');')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'mergePatchPairs')
fprintf(f1,'\n')
fprintf(f1,'(')
fprintf(f1,'\n')
fprintf(f1,');')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'// ************************************************************************* //')
fclose(f1)
endfunction2) ControlDict
function []=ControlDict_function(tend,dt,interval,text_address)
f1=fopen(text_address,'w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n')
fprintf(f1,' ========= |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox')
fprintf(f1,'\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org')
fprintf(f1,'\n')
fprintf(f1,' \\ / A nd | Version: 7')
fprintf(f1,'\n')
fprintf(f1,' \\/ M anipulation |')
fprintf(f1,'\n')
fprintf(f1,'/*---------------------------------------------------------------------------*/')
fprintf(f1,'\n')
fprintf(f1,'FoamFile')
fprintf(f1,'\n')
fprintf(f1,'{')
fprintf(f1,'\n')
fprintf(f1,' version 2.0;')
fprintf(f1,'\n')
fprintf(f1,' format ascii;')
fprintf(f1,'\n')
fprintf(f1,' class dictionary;')
fprintf(f1,'\n')
fprintf(f1,' location "system";')
fprintf(f1,'\n')
fprintf(f1,' object controlDict;')
fprintf(f1,'\n')
fprintf(f1,'}')
fprintf(f1,'\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'application icoFoam;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'startFrom startTime;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'startTime 0;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'stopAt endTime;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'endTime %f;',tend)
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'deltaT %5.5f;',dt)
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'writeControl timeStep;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'writeInterval %f;',interval)
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'purgeWrite 0;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'writeFormat ascii;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'writePrecision 6;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'writeCompression off;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'timeFormat general;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'timePrecision 6;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'runTimeModifiable true;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'// ************************************************************************* //')
fclose(f1)
endfunction3)Pressure file
function []=Pressure_function_wedge(dp,text_address)
pend=0;
p0=pend+dp;
f1=fopen(text_address,'w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n')
fprintf(f1,' ========= |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox')
fprintf(f1,'\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org')
fprintf(f1,'\n')
fprintf(f1,' \\ / A nd | Version: 7')
fprintf(f1,'\n')
fprintf(f1,' \\/ M anipulation |')
fprintf(f1,'\n')
fprintf(f1,'/*---------------------------------------------------------------------------*/')
fprintf(f1,'\n')
fprintf(f1,'FoamFile')
fprintf(f1,'\n')
fprintf(f1,'{')
fprintf(f1,'\n')
fprintf(f1,' version 2.0;')
fprintf(f1,'\n')
fprintf(f1,' format ascii;')
fprintf(f1,'\n')
fprintf(f1,' class volScalarField;')
fprintf(f1,'\n')
fprintf(f1,' object p;')
fprintf(f1,'\n')
fprintf(f1,'}')
fprintf(f1,'\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'dimensions [0 2 -2 0 0 0 0];')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'internalField uniform 0;')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'boundaryField')
fprintf(f1,'\n')
fprintf(f1,'{')
fprintf(f1,'\n')
fprintf(f1,' inletWall')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type fixedValue;')
fprintf(f1,'\n')
fprintf(f1,' value uniform %5.10f;',p0)
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' outletWall')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type fixedValue;')
fprintf(f1,'\n')
fprintf(f1,' value uniform %1.3f;',pend)
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' fixedWalls')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type zeroGradient;')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' front')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type wedge;')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' Back')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type wedge;')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,' axis')
fprintf(f1,'\n')
fprintf(f1,' {')
fprintf(f1,'\n')
fprintf(f1,' type empty;')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,' }')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'// ************************************************************************* //')
fclose(f1)
endfunction4)Transport Properties
function []=transport_properties_function(nu,text_address)
f1=fopen(text_address,'w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n')
fprintf(f1,' ========= |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox')
fprintf(f1,'\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org')
fprintf(f1,'\n')
fprintf(f1,' \\ / A nd | Version: 7')
fprintf(f1,'\n')
fprintf(f1,' \\/ M anipulation |')
fprintf(f1,'\n')
fprintf(f1,'/*---------------------------------------------------------------------------*/')
fprintf(f1,'\n')
fprintf(f1,'FoamFile')
fprintf(f1,'\n')
fprintf(f1,'{')
fprintf(f1,'\n')
fprintf(f1,' version 2.0;')
fprintf(f1,'\n')
fprintf(f1,' format ascii;')
fprintf(f1,'\n')
fprintf(f1,' class dictionary;')
fprintf(f1,'\n')
fprintf(f1,' location "constant";')
fprintf(f1,'\n')
fprintf(f1,' object transportProperties;')
fprintf(f1,'\n')
fprintf(f1,'}')
fprintf(f1,'\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //')
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'nu [0 2 -1 0 0 0 0] %3.10f;',nu)
fprintf(f1,'\n')
fprintf(f1,'\n')
fprintf(f1,'// ************************************************************************* //')
fprintf(f1,'\n')
fprintf(f1,'\n')
fclose(f1)
endfunctionAll the specifications are used and a code is drafted to use all the above mentioned functions
clear all
close all
clc
%water is chosen
viscosity=1.002*10^-3;
density=998;
nu=viscosity/density
Re=1000;
R=0.05;
L=0.06*Re*2*R
dp=(4*(viscosity^2)*L*Re)/(density^2*(R^3))
theta=4;
%Expected result
vmax=(nu*Re)/R
dx=0.05
dy=0.002
dz=0.002
approx_CFL=0.01;
dt_approx=approx_CFL*dx/vmax
dt=0.01 %selected
predicted_CFL=vmax*dt/dx
%Meshing
no_x=ceil(L/dx)
no_y=ceil(R/dy)
no_z=ceil(2*R*sind(theta/2)/dz)
x_grading=1;
y_grading=0.5;
z_grading=1;
%ControlDict
tend=2000;
interval=1000;
text_address='/home/siddhesh/OpenFOAM/siddhesh-7/run/pipeflow/pipe_wedge/system/blockMeshDict';
blockMesh_wedge_function(R,L,theta,text_address,no_x,no_y,no_z,x_grading,y_grading,z_grading)
text_address='/home/siddhesh/OpenFOAM/siddhesh-7/run/pipeflow/pipe_wedge/system/controlDict';
ControlDict_function(tend,dt,interval,text_address)
text_address='/home/siddhesh/OpenFOAM/siddhesh-7/run/pipeflow/pipe_wedge/0/p';
Pressure_function_wedge(dp,text_address)
text_address='/home/siddhesh/OpenFOAM/siddhesh-7/run/pipeflow/pipe_wedge/constant/transportProperties';
transport_properties_function(nu,text_address)Now the respective files for each of the code are:-
1)blockMeshDict
/*--------------------------------*- 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.0498782025 -0.0034878237)
(0 0.0498782025 0.0034878237)
(6.000 0 0)
(6.000 0.0498782025 -0.0034878237)
(6.000 0.0498782025 0.0034878237)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (120 25 2) simpleGrading (1.000000 0.500000 1.000000)
);
edges
(
);
boundary
(
inletWall
{
type patch;
faces
(
(0 0 2 1)
);
}
outletWall
{
type patch;
faces
(
(3 4 5 3)
);
}
fixedWalls
{
type wall;
faces
(
(4 1 2 5)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
Back
{
type wedge;
faces
(
(0 1 4 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //2) ControlDict 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;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 2000.000000;
deltaT 0.01000;
writeControl timeStep;
writeInterval 1000.000000;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //3) Pressure 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inletWall
{
type fixedValue;
value uniform 0.0001935422;
}
outletWall
{
type fixedValue;
value uniform 0.000;
}
fixedWalls
{
type zeroGradient;
}
front
{
type wedge;
}
Back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //4) Velocity file (This is not written via a matlab code but is user specified prior to commencement of simulation)
/*--------------------------------*- 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 0 0);
boundaryField
{
inletWall
{
type zeroGradient;
}
outletWall
{
type zeroGradient;
}
fixedWalls
{
type noSlip;
}
front
{
type wedge;
}
Back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //5) Transport Properties 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;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.0000010040;
// ************************************************************************* // Output
Meshed Domain
Velocity Contour:
Velocity profile at x=6
Vmax = 0.01974 m/s
The Maximum courant number that appears in the simulation
Maximum Courant Number=0.003948
Shear Stress Plot
Maximum stress=0.000752 N/m^2
Conclusion:-
The Comparison between theoretical and simulation values
As the above comparison displays, the simulation results conform very closely with theoretical values. The shear stress plot is a linear plot for major part of the cross-section but changes its plot right at the pipe walls. The reason is that the velocity value sharply decreases near the walls which the mesh was not able to capture with utmost accuracy. To do that mesh needs to be finer and high simple-grading factor need to be used to capture sharp fall in velocity at pipe wall. But to reduce processing effort, such high quality meshing was not done.
Since this transient simulation is used to reach steady state, it is important to run the simulation for large number of time steps to ensure that the system has reached steady state. That is the reason the simulation was run for 2000 seconds. And the percent error between simulation and theoretical values are calculated in the table.
So in conclusion the wedge has been shown to be efficient substitute for a full cross-section pipe for axisymmetric simulation. Simulation results for wedge conforms well with theoretical values for the pipe.