Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective: Simulation of flow through pipe.
You need to calculate the length of the pipe
- Calculate length of the pipe using the entry length formula for laminar flow through a pipe
- Show that entry length is sufficient to produce a fully developed flow. You will have to compare velocity profiles at different positions along the length of the pipe.
Given: 1. Working fluid: Water.
2. Reynolds number based on pipe diameter and inlet velocity = 2100
3. wedge angle should be less than 5 degrees
Theory:
Hagen poiseuille\'s Law:
In non ideal fluid dynamics, Hagen-poiseuille equation is a physical law that gives pressure drop in an incompressible and newtonian fluid in a laminar flow through a long cylinder pipe of constant cross section.
1. Shear Stress distribution(τ)(τ):
τ=(−∂p∂x)⋅r2τ=(-∂p∂x)⋅r2
where, −∂p∂x-∂p∂x indicates decrease in pressure in the direction of flow. The pressure must decrease because pressure force is only means avaliable to compensate for resistance to flow, the Potential and kinetic energy remain constant.
In laminar flow, shear stress is entirly due to viscous action.
2. Velocity distributution(U):
U=(14μ)(−∂p∂x)(R2−r2)U=(14μ)(-∂p∂x)(R2-r2)
where, R is inner radius of pipe and r is radius of pipe less than R.
Maximum velocity occurs at axis and is given by:
Umax=(14μ)(−∂p∂x)R2Umax=(14μ)(-∂p∂x)R2
3. Pressure drop in pipe:
P1−P2=32μ⋅Uavg⋅LD2P1-P2=32μ⋅Uavg⋅LD2
where, U_avg= 0.5 * U_max, L is total length of pipe
Consider a flow entering a pipe. Let us think of the entering flow being uniform, so inviscid. As soon as the flow \'hits\' the pipe many changes take place. The most important of these is that viscosity imposes itself on the flow and the \"No Slip\" condition at the wall of the pipe comes into effect. Consequently the velocity components are each zero on the wall, ie., u = v = 0. The flow adjacent to the wall decelerates continuously. We have a layer close to the body where the velocity builds up slowly from zero at wall to a uniform velocity towards the center of the pipe. This layer is what is called the Boundary Layer. Viscous effects are dominant within the boundary layer. Outside of this layer is the inviscid core where viscous effects are negligible or absent.
The boundary layer is not a static phenomenon; it is dynamic. it grows meaning that its thickness increases as we move downstream. From Fig. it is seen that the boundary layer from the walls grows to such an extent that they all merge on the centreline of the pipe. Once this takes place, inviscid core terminates and the flow is all viscous. The flow is now called a Fully Developed Flow. The velocity profile becomes parabolic. Once the flow is fully developed the velocity profile does not vary in the flow direction. In fact in this region the pressure gradient and the shear stress in the flow are in balance. The length of the pipe between the start and the point where the fully developed flow begins is called the Entrance Length. Denoted by Le, the entrance length is a function of the Reynolds Number of the flow. In general,
LeD≈0.06⋅ReLeD≈0.06⋅Re, for a Laminar Flow
LeD≈4.4⋅Re16LeD≈4.4⋅Re16, for a Turbulent Flow
At critical condition, i.e., Re =2300, the Le/d for a laminar flow is 138. Under turbulent conditions it ranges from 18(at Re = 4000) to 95 (at Re=108)
Application decides whether a long enough entrance length is required or a shorter one is required. Take wind tunnel for instance. Here the aim is to have a uniform flow over the model in the test section. So one desires to have a longer entrance length.
MATLAB code to create blockMeshDict file based on requirement of wedge angle:
Here, we consider wedge angle is 4 degree.
clear all
close all
clc
%% Input Parameters %%
Re= 2100; % Reynold\'s Number
D= 0.02; % Diameter of pipe in meters
R=D/2; % Radius of pipe
theta= input(\'Enter value for Wedge angle =\') % Wedge angle
mu=1.002e-3; %Dynamic viscocity of water @ 20 degree Celsius in SI units
rho=998.21; %Density of pure water highest @ 20 degree Celsius in SI units
nu= mu/rho; %Kinematic viscocity of water @ 20 degree Celsius in SI units
% Analytical calculation of Velocity,pressure and entrance length
L_e= 0.06*Re*D; % Entrance length from inlet of pipe
L_t= L_e + 0.25*L_e; % Assumed Total length of pipe
v_avg= (Re*nu)/D; % Average velocity of water throughtout characteristic length
v_max= 2*v_avg; % Average velocity of water at center of pipe characteristic length
Delta_P= (32*mu*v_avg*L_t)/D^2; % Pressure Drop
kin_P= Delta_P/rho; % Kinematic Pressure
% Finding vertices
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_t 0 0]
v4= [L_t R*cosd(theta/2) -R*sind(theta/2)]
v5= [L_t R*cosd(theta/2) R*sind(theta/2)]
%% Script to form blockMeshDict file.
f1= fopen(\'blockMeshDict.txt\',\'w\') % To create and write blockmesh file
fprintf(f1,\'%s\',\'/*--------------------------------*- C++ -*----------------------------------*\\\');
fprintf(f1,\'\\n========= |\\n\');
fprintf(f1,\'\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox\\n\');
fprintf(f1,\' \\\\ / O peration | Website: https://openfoam.org\\n\');
fprintf(f1,\' \\\\ / A nd | Version: 6\\n\');
fprintf(f1,\' \\\\/ M anipulation |\\n\');
fprintf(f1,\'%s\',\'\\*---------------------------------------------------------------------------*/\');
fprintf(f1,\'\\nFoamFile\\n\');
b1=blanks(5);
b2=blanks(6);
b3=blanks(7);
b4=blanks(4);
fprintf(f1,\'{\\n\');
fprintf(f1,\'%s\',\' version\',b1,\'2.0;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',\' format\',b2,\'ascii;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',\' class\',b3,\'dictionary;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',\' object\',b2,\'blockMeshDict;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'}\\n\');
fprintf(f1,\'%s\',\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\');
fprintf(f1,\'\\n\\nconvertToMeters 1;\\n\\n\');
% To form vertices
fprintf(f1,\'vertices\\n\');
fprintf(f1,\'(\\n\');
fprintf(f1,b4,\'%s\',\'(\');
v_0= fprintf(f1,\'(%f %f %f)\\n\',v0(1),v0(2),v0(3));
fprintf(f1,b4,\'%s\',\'(\');
v_1= fprintf(f1,\'(%f %f %f)\\n\',v1(1),v1(2),v1(3));
fprintf(f1,b4,\'%s\',\'(\');
v_2= fprintf(f1,\'(%f %f %f)\\n\',v2(1),v2(2),v2(3));
fprintf(f1,b4,\'%s\',\'(\');
v_3= fprintf(f1,\'(%f %f %f)\\n\',v3(1),v3(2),v3(3));
fprintf(f1,b4,\'%s\',\'(\');
v_4= fprintf(f1,\'(%f %f %f)\\n\',v4(1),v4(2),v4(3));
fprintf(f1,b4,\'%s\',\'(\');
v_5= fprintf(f1,\'(%f %f %f)\\n\',v5(1),v5(2),v5(3));
fprintf(f1,\'\\n);\\n\');
% To create blocks
fprintf(f1,\'\\nblocks\\n\');
fprintf(f1,\'(\\n\');
fprintf(f1,b4,\'%s\',\'(\');
fprintf(f1,\'%s\',\'hex (0 1 2 0 3 4 5 3)\');
% to create mesh
m1= input(\'Enter no. of mesh element along X-axis:\')
m2= input(\'Enter no. of mesh element along Y-axis:\')
m3= input(\'Enter no. of mesh element along Z-axis:\')
fprintf(f1,blanks(1),\'%s\',\'(\');
fprintf(f1,\'(%3.0f %1.0f %3.0f)\',m1,m2,m3);
fprintf(f1,\'%s\',\' simpleGrading (1 1 1)\');
fprintf(f1,\'\\n);\');
% To create edges
fprintf(f1,\'\\nedges\\n\');
fprintf(f1,\'(\\n\');
fprintf(f1,\'%s arc 1 2 (0 %f 0)\',blanks(1),R);
fprintf(f1,\'\\n%s arc 4 5 (%f %f 0)\',blanks(1),L_t,R);
fprintf(f1,\'\\n);\\n\\n\');
% To apply boundary conditions
% For axis
fprintf(f1,\'boundary\\n(\\n\');
fprintf(f1,\'%s\',\' axis\');
fprintf(f1,\'\\n {\\n\');
fprintf(f1,\'%s\',blanks(10),\'type empty;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'faces\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'(\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(13),\'(0 3 3 0)\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\');\');
fprintf(f1,\'\\n }\\n\');
% For fixedwalls
fprintf(f1,\'%s\',\' fixedWalls\');
fprintf(f1,\'\\n {\\n\');
fprintf(f1,\'%s\',blanks(10),\'type wall;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'faces\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'(\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(13),\'(2 1 4 5)\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\');\');
fprintf(f1,\'\\n }\\n\');
% For front
fprintf(f1,\'%s\',\' front\');
fprintf(f1,\'\\n {\\n\');
fprintf(f1,\'%s\',blanks(10),\'type wedge;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'faces\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'(\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(13),\'(1 0 3 4)\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\');\');
fprintf(f1,\'\\n }\\n\');
% For back
fprintf(f1,\'%s\',\' back\');
fprintf(f1,\'\\n {\\n\');
fprintf(f1,\'%s\',blanks(10),\'type wedge;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'faces\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'(\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(13),\'(2 5 3 0)\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\');\');
fprintf(f1,\'\\n }\\n\');
% For inlet
fprintf(f1,\'%s\',\' inlet\');
fprintf(f1,\'\\n {\\n\');
fprintf(f1,\'%s\',blanks(10),\'type patch;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'faces\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'(\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(13),\'(0 1 2 0)\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\');\');
fprintf(f1,\'\\n }\\n\\n\');
% For outlet
fprintf(f1,\'%s\',\' outlet\');
fprintf(f1,\'\\n {\\n\');
fprintf(f1,\'%s\',blanks(10),\'type patch;\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'faces\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\'(\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(13),\'(3 5 4 3)\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\',blanks(10),\');\');
fprintf(f1,\'\\n }\\n\\n\');
fprintf(f1,\'\\n);\\n\');
fprintf(f1,\'\\nmergePatchPairs\');
fprintf(f1,\'\\n(\\n);\\n\\n\');
fprintf(f1,\'%s\',\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\');
fclose(f1)