Simulate Flow through pipe & verify results by coding on Matlab

AIM:

To Simulate flow through Pipe for different Reynolds number, to perform Grid dependency test and to verify the results by coding using MATLAB

DESCRIPTION AND GOVERNING EQUATIONS:

Types of flow:

Laminar flow:

Laminar flow is characterized by smooth streamlines and highly ordered motion. If the pipe is sufficiently long (relative to the entry length) then the entrance effects are negligible and therefore the flow is fully developed. Laminar flow occurs when the fluid flows in parallel layers without mixing. The velocity of the fluid is constant at any given moment. Since the flow is steady, there is no acceleration. The flow is laminar for cylindrical pipes when Reynolds number is less than 2300.

Laminar Flow Turbulent Flow

Transitional flow:

The transition from laminar to turbulent flow is not sudden. It occurs over some region in which turbulent flow in the center of the pipe and laminar flow is near the edges of the pipe. The flow fluctuates between laminar and turbulent flows before it becomes fully turbulent. The flow is considered transitional when Reynolds number is in between 2300 and 4000

Turbulent flow:

Turbulent flow is characterized by velocity fluctuations and highly disordered motion. Most flows encountered in practice are turbulent. Turbulent flow occurs when streamlines of the liquid are irregular and change over time. The paths of the fluid flow are also irregular and form tiny whirlpool regions. The flow is turbulent when Reynolds number is greater than 4000. In practice, most flows in engineering are turbulent. However, the theory of turbulent flow remains underdeveloped since this flow is a very complex mechanism dominated by fluctuations. Therefore, turbulent flow is analyzed by applying experimental measures.

Reynolds number:

Reynolds number can reveal whether flow is laminar or turbulent. The transition from laminar to turbulent flow depends on the surface roughness, flow velocity, geometry, surface temperature, and type of fluid, among others. Flow regime mainly depends on the ratio of inertial forces to viscous forces in the fluid. This ratio is called Reynolds number and is expressed as

Re = Inertial Force / Viscous Force

`Re = (ρV_(avg) D)/μ = (V_(avg) D)/ϑ`

Here

`ρ` is the density of the fluid (SI units: kg/m^3)

`V_(avg)` is the average flow speed (m/s)

D is the diameter of the tube (m)

μ is the dynamic viscosity of the fluid (Pa·s or N·s/m^2 or kg/(m·s))

`ϑ` is the kinematic viscosity of the fluid (m^2/s).

Critical Reynolds number Recr is the number at which the flow becomes turbulent.

The value of this number is different for different geometries and flow conditions. The generally accepted value of the critical Reynolds number for internal flow in a circular pipe is Recr = 2300. Under most practical conditions,

the flow in a circular pipe is:

Laminar when Re ≤ 2300

Transitional when 2300 ≤ Re ≤ 4000

Turbulent when Re ≥ 4000

Entrance region:

There is a no-slip region and hence the region near the pipe gets deaccelerated due to the viscous friction or viscous shear.

The region gets affected by the viscosity as you move downstream.

Velocity boundary layer: The region of the flow in which the effects of the viscous shearing forces caused by fluid viscosity are felt.

Boundary layer region: The viscous effects and the velocity changes are significant.

Irrotational(core) / Inviscid flow region: The frictional effects are negligible and the velocity remains essentially constant in the radial direction

Since in the velocity boundary layer region the flow gets deaccelerated that means in the inviscid region flow needs to be accelerated for the mass flow to keep going through and ultimately the boundary layer will grow and then it will get merged, and after that region, we would get a classic parabolic profile

In a fully developed region the velocity profile does not change and remains parabolic, the flow has no acceleration.in this region, the Analytical solution exists.

In developing region we do not have any analytical solution but we can predict the length of the developing region which is called the entrance length Le.

Entry Length (Le):

The length of the hydrodynamics entrance region is called entry length.

It may also be taken as the distance from the fluid entrance to 2% of the fully developed wall shear stress value.

In laminar flow, the hydrodynamic entry length is given as 𝐿ℎ−𝑙𝑎𝑚𝑖𝑛𝑎𝑟 ≅ 0.05𝑅𝑒𝐷

In turbulent flow: 𝐿ℎ−𝑡𝑢𝑟𝑏𝑢𝑙𝑒𝑛𝑡 ≅ 10𝐷

OBJECTIVES:

1) To Run a pipe flow simulation with an inlet Reynolds number of 100,1000 and 10,000.

For each of these cases to do the following

1. To Place line probes at 95%, 90%, and 85% of the pipe length.
2. To Compare the normalized velocity profile at each of these locations
3. To Normalize the velocity profile by the inlet velocity.

2) To verify the results by using an analytical solution by coding in MATLAB.

USER INPUTS:

Geometry:

Diameter of the pipe = 0.04m

Length of the pipe = 1m

Thickness = 0.001m

Lids are created at the ends of the pipe to make it watertight

Line Sketch is created for line probes at the inlet, 0.95m, 0.9 and 0.85m from the inlet along the diameter of the pipe

General Settings:

Analysis Type: Internal

Steady-state analysis ( since the solution does not change with time here)

Fluids:

Liquids > Water

Since the initial temperature is T=293.2K and Pressure P=101325Pa

Density=1000kg/m^3

Dynamic viscosity=1e-3

Flow Type= Laminar and Turbulent.

Wall conditions: Adiabatic Wall

Roughness: 0micrometer

Initial conditions:

Pressure: P=101325Pa

initial temperature is T=293.2K

Mesh:

Global mesh > edit definition > manual

Basic Mesh

Nx=10

Ny=10

Nz=50

Later after the Mesh/ Grid dependency test is done and a Mesh is picked based on the results

MATHEMATICAL MODEL:

1) Governing Equations

2) Boundary conditions.

1) Governing Equations:

1) Conservation of mass (continuity equation):

`∇ overlineV = 0`

2) Conservation of momentum (`overlineF =moverlinea`)

`ρ(overlineV .∇) overlineV =-∇p+μ∇^2 overlineV`

there are three equations one in the axial, radial and along theta

we are having 4 unknown quantities which can be obtained by solving above 4 equations

`V_r (r,z,θ) , V_z (r,z,θ) , V_θ (r,z,θ) , p(r,z,θ)`

Assumptions:

1) Steady

2) incompressible

3) Laminar

In the case of the turbulent flow, the equations will almost be the same, we have to consider the fluctuating components.

The mathematical model will be solved along with the boundary conditions to get a system of algebraic equations in the cell centers, and inversion of these equations will give us the cell center values of and after obtaining the values we can do the post-processing.

Hit Run calculation

After solving we will get selected variables at selected points and from that we can do Post-processing.

2) Boundary conditions:

Inlet boundary:

For Re=100; Inlet Velocity= 0.0025m/s

For Re=1000; Inlet Velocity= 0.025m/s

For Re=10000; Inlet Velocity= 0.25m/s

Outlet boundary:

Static pressure P=101325 Pa

Pipe wall:

Real Wall condition

GRID DEPENDENCY TEST:

This test is performed by taking the Case 1 when Re=100 and velocity is 0.0025m/s

Project name > New Parametric study

Input Variables:

Mesh divisions along x,y, and z are selected.

Output Variables:

Various cutplots and xy plots are selected.

Scenario:

RESULTS AND PLOTS FORM GRID DEPENDENCY TEST:

VELOCITY CUTPLOT FOR DESIGN POINT 1:

PRESSURE CUTPLOT FOR DESIGN POINT 1:

VELOCITY CUTPLOT FOR DESIGN POINT 2:

PRESSURE CUTPLOT FOR DESIGN POINT 2:

VELOCITY CUTPLOT MESH FOR DESIGN POINT 3:

PRESSURE CUTPLOT FOR DESIGN POINT 3:

XY PLOTS:

Velocity Plot along Line Probes at different Locations:

1) 0.95m from Inlet:

2) 0.90m from Inlet:

1) 0.85m from Inlet:

SHEAR STRESS PLOT ALONG THE WALL:

CENTERLINE VELOCITY PLOT:

CENTERLINE PRESSURE PLOT:

By performing a Grid dependency test, we can see from cutplots that velocity profile for design point 1 is quite nonuniform and contours are quite choppy, but in case of design point 2 the velocity profiles become more uniform and in case of design point 3 the plots improve further.

And from the x-y plots of design point 1 we can see a difference in values of the velocities due to coarser mesh when compared to the case of design points 2 and 3. Due to these reasons, we have chosen the Mesh of design point 2 for further study as it is best in terms of the results and it will take less computational time compared to the mesh from design point 3.

We would use the Mesh of design point 2 for further study as it gives the best result in terms of the results and the less amount of computational time required.

PARAMETRIC STUDY:

VARIATION OF THE INLET VELOCITY TO GET THE FLOWS WITH DIFFERENT REYNOLDS NUMBER:

Project name > New Parametric study

Input Variables:

Inlet velocity is selected

Output Variables:

Various cutplots and xy plots are selected.

Scenario:

Inlet velocity is varied to get the Re number as 100,1000,10000

DESIGN POINT 1: Inlet velocity is 0.0025m/s

Velocity cutplots:

Pressure Cutplot:

DESIGN POINT 2: Inlet velocity is 0.025m/s

Velocity cutplots:

Pressure cutplot:

DESIGN POINT 3: Inlet velocity is 0.25m/s

Velocity cutplots:

Pressure cutplot:

Normalised Velocity curves:

Normalised velocity V/s Normalised Radius:

1) Line probe 0.95m from inlet:

2) Line probe 0.9m from inlet:

3) Line probe 0.85m from inlet:

ANIMATIONS:

Velocity cutplot:

Flow simulation:

MATLAB CODES:

CALCULATION:

FOR LAMINAR FLOW THROUGH PIPES:

1) The relationship between local velocity and the maximum centerline velocity in case of laminar flow is

`u=u_max (1-(r/R)^2)`

we can see that the velocity profile is of parabolic form.

2) Relationship between the mean and maximum velocity is given as

`u ̅=u_max/2`

FOR TURBULENT FLOW THROUGH PIPES:

For turbulent pipe flow among the numerous empirical velocity profiles, the simplest and the best known is the power-law velocity profile expressed as

`u=u_max (1-r/R)^(1/n)`

where R is the pipe radius and n is a parameter (exponent) which depends of the Reynolds number i.e.

n=n(Re).

An empirical approximation of this dependence can be expressed as

`n(Re)=c.ln(Re)`

where c=0.62 for Re ≤ 2*(Re)^5

and c=0.65 for Re > 2*(Re)^5

Based on this approach on obtained the following relationship between the average velocity and maximum velocity for turbulent flow

`(u ̅)/u_max =(2(n)^2)/((n+1).(2n+1))`

It is to mention that

-the value n=7 is applicable to a wide range of pipe flows and is the one commonly used resulting

`(u ̅)/u_max =0.817`

-the empirical power law cannot be used to calculate the wall shear stress

CODE 1: Re=100

% To Calculate the centerline velocity and the velocity at any point %
% Re 100

clc
clear all
close all

u_avg= 0.0025                                                              %average inlet velocity in m/sec
R=0.02;                                                                    %Radius of the pipe
r= linspace(-0.02,0.02,100);


for i=1:100
    u_max =2*(u_avg);                                                      % maximum velocity in lamimar flow
    u(i) = u_max *(1-(r(i)/R)^2);                                          % local velocity in laminar flow
    % r is the cetain distance form the pipe center
end
u_max

figure(1)
plot(r,u,'linewidth',3)
title('Re=100 Velcoity plot')
ylabel('Local Veclotiy [m/s]')
xlabel('Diameter [m]')

ones_matix= ones(100);
u_avg= u_avg*(ones_matix);


figure(2)
hold on
plot(r,u,'linewidth',3)
plot(r,u_avg,'linewidth',3)
title('Re=100 Comparison of Velcoity profiles')
xlabel('Diameter [m]')
ylabel('Veclotiy [m/s]')
hold off


%for Normalised Velocity profile
figure(3)
plot(r/0.02,u/0.0025,'linewidth',3)
title('Re=100 Normalised Velcoity Profile')
ylabel('u/(u-avg)')
xlabel('r/R')

RESULTS AND PLOTS FROM CODE 1:

Velocity Plot in the developed region V/s Diameter of Pipe:

Comparison of Velocity Plot in the developed region V/s Diameter of Pipe:

Normalised Velocity Plot in the developed region V/s Diameter of Pipe:

CODE 2: Re=1000

% To Calculate the centerline velocity and the velocity at any point %
% Re 1000

clc
clear all
close all

u_avg= 0.025         %average inlet velocity in m/sec
R=0.02;    %Radius of the pipe
r= linspace(-0.02,0.02,100);


for i=1:100
    u_max =2*(u_avg);                    % maximum velocity in lamimar flow
    u(i) = u_max *(1-(r(i)/R)^2);        % local velocity in laminar flow
    % r is the cetain distance form the pipe center
end
u_max

figure(1)
plot(r,u,'linewidth',3)
title('Re=1000 Velcoity plot')
ylabel('Local Veclotiy [m/s]')
xlabel('Diameter [m]')

ones_matix= ones(100);
u_avg= u_avg*(ones_matix);

figure(2)
hold on
plot(r,u,'linewidth',3)
plot(r,u_avg,'linewidth',3)
title('Re=1000 Comparison of Velcoity profiles')
xlabel('Diameter [m]')
ylabel('Veclotiy [m/s]')
hold off

%for Normalised Velocity profile
figure(3)
plot(r/0.02,u/0.025,'linewidth',3)
title('Re=1000 Normalised Velcoity profile')
ylabel('u/(u-avg)')
xlabel('r/R')

RESULTS AND PLOTS FROM CODE 2:

Velocity Plot in the developed region V/s Diameter of Pipe:

Comparison of Velocity Plot in the developed region V/s Diameter of Pipe:

Normalised Velocity Plot in the developed region V/s Diameter of Pipe:

CODE 3: Re=10,000

% To Calculate the centerline velocity and the velocity at any point %
% Re 10000

%For turbulent pipe flow among the numerous empirical velocity profiles, the simplest and the best known is the power-law velocity profile expressed as
clc
clear all
close all

u_avg= 0.25         %average inlet velocity in m/sec
R=0.02;    %Radius of the pipe
r= linspace(-0.02,0.02,100);

Re=10000;
c=0.62;       %since Re≤2*Re^5 condition satisfies

%An empirical approximation of this dependence can be expressed as
n=c*log(Re)   


for i=1:100
    % Based on this approach on obtained the following relationship between the average velocity and maximum velocity for turbulent flow
    u_max =( u_avg*((n+1)*(2*n+1)/(2*n^2))); % maximum velocity in lamimar flow
    
    % local velocity in laminar flow
    u(i) = (u_max *(1-abs(r(i)/R))^(1/n) );  % local velocity in laminar flow
    % r is the cetain distance form the pipe center
    % R is the pipe radius and n is a parameter (exponent) which depends of the Reynolds number i.e. n=n(Re).
end
u_max

figure(1)
plot(r,u,'linewidth',3)
title('Re=10000 Velcoity plot')
ylabel('Local Veclotiy [m/s]')
xlabel('Diameter [m]')

ones_matix= ones(100);
u_avg= u_avg*(ones_matix);


figure(2)
hold on
plot(r,u,'linewidth',3)
plot(r,u_avg,'linewidth',3)
title('Re=10000 Comparison of Velcoity profiles')
xlabel('Diameter [m]')
ylabel('Veclotiy [m/s]')
hold off

%for normalised velocity profile
figure(3)
plot(r/0.02,u/0.25,'linewidth',3)
title('Re=10000 Normalised Velcoity profile')
ylabel('u/(u-avg)')
xlabel('r/R')

RESULTS AND PLOTS FROM CODE 3:

Velocity Plot in the developed region V/s Diameter of Pipe:

Comparison of Velocity Plot in the developed region V/s Diameter of Pipe:

Normalised Velocity Plot in the developed region V/s Diameter of Pipe:

CONCLUSION:

1) By performing a Grid dependency test, we can see from cutplots that velocity profile for design point 1 is quite nonuniform and contours are quite choppy, but in case of design point 2 the velocity profiles become more uniform and in case of design point 3 the plots improve further. And from the x-y plots of design point 1 we can see a difference in values of the velocities due to coarser mesh when compared to the case of design points 2 and 3. Due to these reasons, we have chosen the Mesh of design point 2 for further study as it is best in terms of the results and it will take less computational time compared to the mesh from design point 3. We would use the Mesh of design point 2 for further study as it gives the best result in terms of the results and the less amount of computational time required

2) From the Shear stress plot obtained we can see that at the entrance there is a small amount of shear stress which decreases as the flow becomes developed and get converged to a certain value.

3) From the Centreline velocity plot, we can see that the velocity starts from an average value at the entrance and increases to a certain maximum value and becomes stable in the developed region.

4) From the velocity cutplots, we can see that as the velocity increase the effect of the boundary layer formation decreases, and hence from the velocity cutplot of 0.0025m/s we can see that a larger boundary layer is formed as compared to the velocity cutplot of 0.25m/s.

5) From all the pressure cutplots we can see that at the entrance the pressure increases which helps the fluid to move to the outlet.

6) Form the normalised velocity plots we can see that the for the laminar flows the ratio of the local velocity to the average velocity (u_max)/(u_avg)=2

And in case of turbulent flows the value of (u_max)/(u_avg)=1.25

And the values obtained for the maximum normalised velocity and the plots are the same as that obtained from the MATLAB

And by probing at different locations we can see that that the velocity profiles increases and becomes more and parabolic as we move towards the outlet in case of laminar flow

And in case of the turbulent flow, the profile becomes less flat as we probe the lines towards the outlet.

7) Comparison of the centreline / maximum velocity values from Solidworks flow simulation and MATLAB

SOLIDWORKS FLOWSIMULATION	MATLAB
0.004917584	0.005
0.045901732	0.05
0.315147018	0.3195

The Results and Plots obtained from SOLIDWORKS MOTION are consistent with that obtained from MATLAB CODING and hence we can say that the simulation that we performed is correct and this verifies the results.

REFERENCES:

1) Practical Analysis Of Turbulent Flow In A Pipe Using Computational Fluid Dynamics V.RajaKiranKumar1 , Dr.H.N.Suresh Kumar2 , Ch.Kiran Kumar3

http://www.ijeijournal.com/papers/Vol.3-Iss.12/L030127781.pdf

2) Computational Fluid Dynamics in Turbulent Flow Applications Alejandro Alonzo-García, Claudia del Carmen Gutiérrez-Torres and José Alfredo Jiménez-Bernal http://dx.doi.org/10.5772/63831

3) FLUID MECHANICS & HYDRAULIC MACHINES by K SUBRAMANYA

https://books.google.co.in/books/about/FLUID_MECHANICS_HYDRAULIC_MACHINES.html?id=gYINgN8iyv0C&redir_esc=y