Week 4.1: Project - Steady state simulation of flow over a throttle body

Throttle 

A throttle is a mechanism by which the fluid flow is managed by construction or obstruction. An engines power can be

increased or decreased by the restriction of inlet air due to the presence of throttle. The term throttle can be used to increase

the efficiency of the car, the car accelerator peddal can be used as throttle to increase the efficiency of the engine. It is used

as throttle in the field of aviation, known as thrust lever in jet engine powered aircraft. For steam locomotive the valve which

controls the regulation of gas is called as reguator.

Application of throttle

The throttle is a mechanism which is used to increase or decrease the power of engine by restricting the flow of inlet gases

coming through it. There are numerous applications of the throttle mechanism which are listed below

Internal Combustion Engines

Fuel Injected Engines

Reciprocating Engine Aircraft

Steam locomotives

Automobiles

Case-setup

The elbow. STL file is loaded into the Converge-CFD set-up. 

The internal flow simulation does not involve the usage of the external boundary, the external boundary is of special interest

when there is a CHT simulation setup. The external boundary can be deleted by selecting the triangle option from the ribbon

and then hitting the "D" button of the keyboard by enabling the hotkeys. The other alternative is selecting the "Repair" option

from the "Geometry Dock" window and selecting the entity and hitting "Apply". 

The diagnosis dock is selected from the "View" option and the "Find" option is selected. It is found that there are 211 open

edges in the geometry ("Intersections(0)","Nonmanifold Problems(0)","Open Edges(211)", "Overlapping Tris(0)","Normal

Orientation(0)","Isolated tris (0)".

To remove open edges from the geometry, the "Repair" option in the "Geometry Dock" is selected, under "Patch" the "Open

Edges" are selected viz (inlet and outlet) of the geometry and then "Boundary Flagging" is done to make "Open Edges" in the

resulting geometry to "0".

The "Normal Toggle" is selected from the ribbon, upon selecting the normal toggle it was observed that the normal is pointing

outside the volume, technically the normal should point inside the volume where the fluid is flowing. Therefore in the

geometry dock, the "Transform" option is clicked. The "Normal" tab is selected from the geometry dock and one of the

triangles containing the normal pointing outward is selected. The "Apply" option is selected for removing the normals.

 

The "Case Setup Dock" is selected from the "View" and "Begin Case Setup" is executed.

The "Time" based "Application" type is selected. The "Material" is selected and the predefined mixture is selected as air. The

"Reactant mechanism" is checked off because it is not a combustion problem. The species is selected and the "Apply" is

clicked. Under "Gas simulation" the Equation of state is selected as "Redlich Kwong", the critical temperature is 133K and the

critical pressure is 3770000Pa. The Turbulent Prandtl number is 0.9 and the Turbulent Schmidt number is 0.78 under Global

Transport parameters. The $O_2$

and the $N_2$

are selected as species. Under "Run parameters" the "Steady-state solver" is selected. The "Temporal type" is locked for the

steady-state. The "simulation mode" is selected as "Full Hydrodynamic" because the geometry is simple so while creating the

mesh inside the geometry it solves the NS equation as well, but if the geometry is complex "No hydrodynamic solver" is

selected as such if there is an error it will point out immediately while creating the mesh, otherwise if the hydrodynamic

solver is selected then it will be very tedious to identify the error in the simulation case set up.

Under "Simulation time parameters" end time is chosen as 15000 cycles (steady-state). The initial and minimum time step is

1e-09 s. The maximum time step is 2s, and the maximum convection CFL limit is 1 and the rest are default values.

The "Pressure" "Equation" is selected under the "Solver parameter"  and the "Preconditioner" is selected to "None". The

"Maximum convection CFL limit (final stage)" is set to 0.5. 

Base-grid

 

Boundary conditions

The boundary conditions for the different named selection are tabulated below

In the "Region and Initialization" the mass fraction of $O_2$

and $N_2$

were created and named as "Volumetric Region User Defined" and the same is updated in the boundary condition region

name. These are the initial condition and will be washed out after the solution has converged.

The monitor variables are added in the "Steady-state monitor". Avg Velocity,  Mass flow rate, and Total Pressure has been

added and the target area is selected as output. The minimum number of cycles to be executed for the steady-state solver is

5000 which signifies that the solver will run for a minimum of 5000 cycles if the solution reaches convergence value before

5000, then the simulation will stop at 5000, but if the solution reaches a convergence value after 5000 then it will terminate

at that moment. The solver checks the difference in the value of two consecutive sample sizes (1000) and matches it with the

tolerance value to stop the solution. The converge-CFD software saves a lot of time and computational energy by tracking the

solution when it reaches a steady-state.  

Setting up of Turbulence model

The Realisable $k-\epsilon$

model is selected from the case-setup tree for the above-said simulation. Since the size of the eddies is restricted near the

wall and the maximum size of the eddies is formed away from the wall. It is well suitable for resolving flows in the

logarithmic flow region where y+ ranges from 30 to 300 and flows involving a high Reynolds number.

Determination of flow velocity

First of all, let us assume that the flow is incompressible, the Bernoulli Equation for the incompressible flow states that the

sum of the mechanical, potential, and kinetic energy remains constant, so any increase in one form may result in the

decrease in another form. Therefore for the above surface, the equation can be written as 

Taking the inlet of the pipe as a reference 

Equation 1

$P_{s}1+\frac{1}{2}\rho v_1^2+\rho gh=P_{total}$

 

where $P_{s1}$

is the static pressure at the inlet.

$\frac{\rho v_1^2}{2}$

is the dynamic pressure at the inlet.

$\rho gh$

is the hydrostatic pressure at the inlet, which can be neglected as such the reference line for calculation is taken where the

inlet is positioned.

The sum of static and dynamic pressure at the inlet is termed stagnation pressure which is defined as the pressure at a point

where the fluid is brought completely at rest.

The sum of static, dynamic, hydrostatic pressure is the total pressure at the inlet which is $1.5e+05Pa$

.

Therefore static pressure at the inlet is unknown which is given by 

Equation 2

$P_{s1}=P_{total}-\frac{1}{2}\rho v_1^2$

 

Similarly at the outlet of the pipe,

Equation 3 

$P_{s2}+\frac{1}{2}\rho v_2^2=P_{total}$

 

$P_{s2}$

is the static pressure at the outlet`

$\rho$

is the density of air

$v_2$

 is the velocity at the outlet boundary condition.

The sum of the flow across a streamline is constant, mechanical, and potential energy across a streamline is constant.

Equation 4

$P_{s1}+\frac{\rho v_1^2}{2}+\rho g{h}_1+h_{pump}=P_{s2}+\frac{\rho v_2^2}{2}+\rho g{h}_2+h_{turbine}+H_L$

 

where subscripts 1 & 2 represent the condition at the inlet and the outlet respectively.

The challenge does not mention the pumping head at the inlet, and the turbine head at the outlet respectively. Therefore it

can be omitted. The pumps however are useful in driving the fluid flow from the inlet to the respective destination by

providing the useful head, whereas the turbine is used for extracting the head from the fluid by the turbine.

$H_L$

represents the irreversible head loss between inlet and outlet due to all components of the piping system other than pump or

turbine.

The total head loss can be equated to the sum of major head loss and the minor head loss

Equation 5

$H_{Ltotal}=H_{Lmajor}+H_{Lminor}$

 

$H_{Lmajor}=f\frac{L}{D}\frac{v_1^2}{2g}$

 

$H_{Lminor}=K_L\frac{v_1^2}{2g}$

 

The minor loss coefficient depends on the geometry of the component.

The minor loss due to elbow is 

Equation 6

$H_{Lelbow}=K_{Lelbow}\frac{v_1^2}{2g}$

 

For ${90}^{\circ }elbow$

 $K_{Lelbow}$ is $0.3$

 

The minor loss due to throttle ( the throttle will function as butterfly valve)

Equation 7

$H_{Lthrottle}=K_{Lthrottle}\frac{v_1^2}{2g}$

 

For throttle fully open $K_{Lthrottle}=0.05$

 

Head loss due to obstruction in pipe

The throttle placed in the direction of fluid flow will act as an obstruction to the flow of fluid 

In sections 1-1 the obstruction will have the maximum cross-section area which is "a". "A" is the area of the pipe. The flow

area at section 1-1 is "A-a", the liquid starts to contract after section 1-1. In section 2 the fluid will attain its full velocity

which is "v". The head loss is given by 

Equation 8

$h=\frac{{({}_c-v)}^2}{2g}$

 

Equation 9

The coefficient of contraction= Area at Vena Contracta/Area of Orifice

Applying continuity equation in section C-2

Equation 10

$A_c{v}_c=Av$

 

Equation 11

$v_c=\frac{Av}{C_c(A-a)}$

 

Putting the value of $v_c$

in Equation 8

Equation 12

$h=\frac{v^2}{2g}{[\frac{A}{C_c(A-a)}-1]}^2$

 

"v" in the equation refers to the inlet velocity. The equation can be rewritten as 

$h=\frac{v_1^2}{2g}{[\frac{A}{C_c(A-a)}-1]}^2$

 

The area of obstruction when the throttle valve is fully open is $1.6063e-5m^2$

 

The coefficient of contraction is calculated as $0.93066$

 

The coefficient loss due to obstruction is $0.02388$

 

Summing up all the values of coefficient loss in equation 5 is given as 

Equation 13

$H_{Ltotal}=(f\frac{L}{D}+0.3+0.4+0.02388)\ast \frac{v_1^2}{2g}$

 

Length of pipe $=1.7732396565229E-01m$

 

Diameter of pipe $=0.0171742m$

 

Friction factor for the circular pipe is given as $f=\frac{64}{R_e}$

 

Reynolds number for the flow is given as $R_e=\frac{\rho \ast v_1\ast D}{\mu }$

 

Rewriting Equation 13 

$H_{Ltotal}=(\left(\frac{64\ast \mu }{\rho \ast v_1\ast D}\right)\frac{L}{D}+0.3+0.4+0.02388)\ast \frac{v_1^2}{2g}$

 

Rewriting Equation 4 

$P_{total}=P_{s2}+\rho g{h}_2+\frac{\rho v_2^2}{2}+(\left(\frac{64\ast \mu }{\rho \ast v_1\ast D}\right)\frac{L}{D}+0.3+0.4+0.02388)\ast \frac{v_1^2}{2g}$

 

Equation 14

$150000=100000+1.177\ast 9.81\ast 0.119998+\frac{1.177v_2^2}{2}+(\left(\frac{64\ast 1.846e-05}{1.177\ast v_1\ast 0.017142}\right)\frac{0.177323}{0.0171742}+0.3+0.4+0.02388)\ast \frac{v_1^2}{2\ast 9.81}$

 

where $\rho$

at 300K and 1atm pressure $=1.177\frac{kg}{m^3}$

 

Dynamic viscosity of the air at 300K $1.846e-05\frac{kg}{m-\sec }$

 

$P_{t}otal=150000Pa$

 

$P_{s1}=100000Pa$

 

As such there are two unknowns and one equation which forced us to assume that let $v_1=v_2$

 

Solving the above equation will give the two values of $v_1$

 

$0.6253v_1^2+0.0204v_1-49998.61$

 

Therefore $v_1=282.755,-282.787$

. Rejecting the negative value of $v_1$

because the fluid flows into the pipe at the inlet.

Calculating Mach number for the flow

The Mach number is the dimensionless number used in the fluid flow problem where compressibility is significant.

Equation 15

$M=\frac{v}{c}$

 

where $v$

is the velocity of the fluid in $m{\sec }^{-1}$

 

$c$

is the speed of the sound in $m{\sec }^{-1}$ which is $346.3m{\sec }^{-1}$at $300K$

 

$M=\frac{282.755}{346.3}$

 

$M=0.816$

 

While all flows are compressible, the flows are usually treated as incompressible when ($M<0.3$). Thus in our case, the

Mach number is greater than 0.3 which contradicts our assumption that the flow is incompressible.

Calculating the density of air at the inlet using Redlich Kwong Equation of state 

A cubic equation of state implies an equation which when expanded would contain the volume terms raised to the power first,

second and third respectively. Many of the common two parameters cubic equations can be expressed by the equations

Equation 16

$P=\frac{RT}{V-b}-\frac{a}{T^{0.5}V(V+b)}$

 where $a=0.42748\ast \frac{R^2{T}_c^{2.5}}{P_c}$

$b=0.08664\ast \frac{R{T}_c}{P_c}$

$R=287.05\frac{J}{KgK}$

Critical Temperature $T_c$=$132.63K$

Critical Pressure $P_c$$=3.7858MPa$ 

Absolute temperature = $300K$

Therefore $a=1.88626$

 $b=0.00087128471$

 Putting the values of a and b in Equation 16

$150000=\frac{86115}{V-0.000871}-\frac{1.88626}{17.32V(V+0.000871)}$

 

This is a cube root equation that involves three values of "V"

$2598000V^3-1491511.8V^2-1299.19V+0.00164293$

 

Therefore $V=0.57497,-0.00087,0$

 

Thus the molar volume of the air at the inlet calculated from the Redlich Kwong equation is $0.57497\frac{m^3}{kg}$

neglecting the other two because the molar volume can never be negative and zero.

Therefore density of the air at the inlet is given by ${\rho }_1=\frac{1}{molarvolume}$

$=\frac{1}{0.57497}$$=1.739\frac{kg}{m^3}$

 

Calculating velocity at the inlet using the following relations for compressible flows

In an isentropic process for the compressible flow, the relationship between static pressure and the total pressure is given by

Equation 17

$P_{static}=P_{total}{\{1+\frac{\gamma -1}{2}\ast M^2\}}^{-\left(\frac{\gamma -1}{\gamma }\right)}$

 

where Mach number, $M=\frac{v}{a}$

and $a=\sqrt{\gamma RT}$

 

Rewriting the Equation 9

$P_{1static}=P_{total}{\{1+\frac{\gamma -1}{2}\ast \frac{{}^{}}{\gamma RT}\}}^{-\left(\frac{\gamma }{\gamma -1}\right)}$

 

The total pressure at the inlet boundary is 150000Pa. Let us assume that the static pressure at the inlet boundary to be

140000Pa to calculate the appropriate value of inlet velocity.

${\left(\frac{140000}{150000}\right)}^{\frac{-1}{3.5}}=\{1+\frac{1.4-1}{2}\ast \frac{V^2}{1.4\cdot 287.05\cdot 300}\}$

 Velocity at the inlet =109.54 $\frac{m}{\sec }$

 Using the value of inlet velocity to calculate the total head loss in Equation 13

Therefore $H_L=446.07m$

 Please note that this is an approximation value not the converged value of the velocity at the inlet. The velocity at the inlet is

useful for determining the Reynolds number and corresponding first cell height of the wall.

Calculating Reynolds number for the elbow pipe

Equation 18

Reynolds number is given by $R_e=\frac{\rho \ast V\ast D}{\mu }$

 

where $\rho$

at 300K and 1atm pressure $=1.177\frac{kg}{m^3}$

 

$D$

is the inlet diameter of the pipe, which is $0.0171742m$

 

 Dynamic viscosity of the air at 300K $1.846e-05\frac{kg}{m-\sec }$

 

Therefore $R_e=119954.848474$

 

For an internal flow if $R_e\ge 4000$

then the flow is turbulent.

 Calculating Hydrodynamic Entrance Length $L_e$

for the turbulent flow.

The entrance length in the turbulent flow is much shorter and as expected its dependence on the Reynolds number is weaker.

The non-dimensional hydrodynamic entrance length is approximated as 

Equation 19

$L_e=10\ast D$

$L_e=0.171742m$

Using the value of $L_e$

to calculate the local Reynolds number which is $1199548.48474$

Calculating Turbulent boundary layer thickness

The boundary layer thickness is defined as the thickness at which the viscous velocity becomes 99% of the freestream

velocity.

The boundary layer thickness ($\delta$) calculated from one-seventh power law combined with empirical data for turbulent flow through smooth pipes is given by 

Equation 20 

$\frac{\delta }{x}=\frac{0.38}{R{e}_x^{0.2}}$

 

Therefore $\delta =0.00397060377m$

 

Calculating Skin friction Coefficient

Equation 21

$C_f=\frac{0.059}{{979651}^{0.2}}$

 

$C_f=0.003589$

 

Calculating Wall Shear Stress

Equation 22

${\tau }_w=0.5\ast \rho \ast U^2\ast C_f$

 

${\tau }_w=25.3434N{m}^{-2}$

 

Calculating Frictional Velocity 

Equation 23

$U_t=\sqrt{\frac{{\tau }_w}{\rho }}$

 

$U_t=4.6402m{\sec }^{-1}$

 

Calculating Reynolds number for the throttle body

Equation 24

The throttle will act as a flat plate, the characteristics length of the throttle plate will be along the direction of the

flow.Therefore Reynolds number for the throttle body is 111817. Since the calculated Reynolds number is less

than $R_{ex,cr}=5\ast {10}^5$

and much lesser than the $R_{ex,transition}=50\ast {10}^5$

, therefore the boundary layer around the throttle body is laminar.

Calculating Laminar Boundary Layer Thickness

The laminar boundary layer thickness for the flat plate is given by 

Equation 25

$\frac{\delta }{x}=\frac{4.91}{\sqrt{R{e}_x}}$

$\delta =0.000235m$

 

Calculating Skin friction Coefficient for a laminar flat plate 

Equation 26

$C_{flam}=\frac{0.664}{\sqrt{R_{ex}}}$

$C_{flam}=0.001985$

Calculating Wall Shear Stress for a laminar flat plate 

Equation 27

${\tau }_w=0.5\ast \rho \ast U^2\ast C_f$

 

${\tau }_w=14.01N{m}^{-2}$

 

Calculating Frictional Velocity for a laminar flat plate 

Equation 28

$U_t=\sqrt{\frac{{\tau }_w}{\rho }}$

 

$U_t=3.45m{\sec }^{-1}$

 

Summarizing all the values

 

Calculating first cell height using fixed embedment from Converge-CFD

Sample Calculation

Equation 29

$Embed=\frac{base-grid-size}{2^{Embed-scale}}$

 

Keeping Embed scale to 1, for element size $4e-3$

, the value of Embed is `

$Embed=\frac{0.004}{2}$

 

$Embed=0.002m$

 

First cell height $=0.002m$

 

While dealing with the control volume, $y+$

deals with the first cell centroid from the wall. The first cell height and the first cell centroid can be mathematically expressed

as 

$△y=2y_p$

 

Therefore,

$y_p=\frac{0.002}{2}$

 

$y_p=0.001m$

 

$y+=\frac{\rho \ast U_t\ast y_p}{\mu }$

 

For base grid size $4mm$

 

$y+=295.856$

 

Summarizing values for y+ for different mesh size along pipe elbow

Summarizing values for embed scale and embed layers for different mesh size along throttle body

Results

Mesh-0

Base-grid-size : 4e-3m

The mesh is coarse, and it does not involve the application of fixed embedding.

Mesh-1

Mesh-2

Mesh-3

Mesh-4

Mesh-0

Contour plots

Velocity 

Pressure

Mass flow rate 

Velocity 

Pressure

Mass flow rate 

Total cell count

Mesh-1

Contour plots

Velocity

Pressure

Mass flow rate 

 

Velocity 

Pressure

Mass flow rate 

Total cell count

 

Mesh-2

Contour plots

Velocity

Pressure

Mass flow rate 

 

Velocity

Pressure

Mass flow rate

Total cell count

 

Mesh-3

Contour plots

Velocity

Pressure

Mass flow rate 

 

 

Velocity

Pressure

Mass flow rate 

Total cell count

Mesh-4

Contour plots

Velocity

Pressure

Mass flow rate 

 

Velocity

Pressure 

Mass flow rate 

Total cell count

Steady-state animation of flow over a throttle body

The vector plot was made using Glyph from ParaView where velocity was selected for the arrow, and for the slice, the pressure was selected. The vector plot helps us in identifying the region of separation, recirculation, and flow reattachment to the boundary layer over the throttle body and wherever there is a change in the curvature region of the pipe.

 

Summarizing all the values

Conclusion

The steady-state simulation of flow over a throttle body was carried successfully by the "Converge-CFD" software and the

"Cygwin" command terminal.

The Realizable k-epsilon model was used for modeling the turbulent flow in the pipe. The eddy movement created by the

turbulent flow near the wall is less in size whereas as we move away from the wall of the pipe the size of eddy grows

considerably. The turbulent flow is characterized by unsteady velocity fluctuations and highly disordered motions. 

Flow separation, recirculation, and reattachment zones were observed during simulation. The change in the curvature of the

pipe can be accounted for that. The fluid while moving through the pipe experiences the elbow as such the velocity near the

wall increases, due to which separation occurs, beyond the separation point the fluid tries to reattach itself to the boundary

layer, swirling eddies were also observed. The second separation region is encountered when the fluid strikes the surface of

the throttle body leading to an increase in the static pressure near the front portion of the throttle body. The rear end of the

throttle body experiences the separation region and reattachment region producing negative pressure.

The variation of the turbulent boundary layer thickness was also observed. The straight section of the pipe has more

boundary layer thickness, in elbows, the boundary layer thickness is less.

Y+ values were also calculated, the range of Y+ was  30-300 which captures the turbulent region.