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 O2">O2

and the N2">N2

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 O2">O2

and N2">N2

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-ε">k−ε

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

Ps1+12ρv12+ρgh=Ptotal">Ps1+12ρv21+ρgh=Ptotal

where Ps1">Ps1

is the static pressure at the inlet.

ρv122">ρv212

is the dynamic pressure at the inlet.

ρgh">ρ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">1.5e+05Pa

.

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

Equation 2

Ps1=Ptotal-12ρv12">Ps1=Ptotal−12ρv21

Similarly at the outlet of the pipe,

Equation 3 

Ps2+12ρv22=Ptotal">Ps2+12ρv22=Ptotal

Ps2">Ps2

is the static pressure at the outlet`

ρ">ρ

is the density of air

v2">v2

 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

Ps1+ρv122+ρgh1+hpump=Ps2+ρv222+ρgh2+hturbine+HL">Ps1+ρv212+ρgh1+hpump=Ps2+ρv222+ρgh2+hturbine+HL

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.

HL">HL

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

HLtotal=HLmajor+HLminor">HLtotal=HLmajor+HLminor

HLmajor=fLDv122g">HLmajor=fLDv212g

HLminor=KLv122g">HLminor=KLv212g

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

The minor loss due to elbow is 

Equation 6

HLelbow=KLelbowv122g">HLelbow=KLelbowv212g

For 90∘elbow">90∘elbow

 KLelbow">KLelbow is 0.3">0.3

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

Equation 7

HLthrottle=KLthrottlev122g">HLthrottle=KLthrottlev212g

For throttle fully open KLthrottle=0.05">KLthrottle=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=(vc-v)22g">h=(vc−v)22g

Equation 9

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

Applying continuity equation in section C-2

Equation 10

Acvc=Av">Acvc=Av

Equation 11

vc=AvCc(A-a)">vc=AvCc(A−a)

Putting the value of vc">vc

in Equation 8

Equation 12

h=v22g[ACc(A-a)-1]2">h=v22g[ACc(A−a)−1]2

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

h=v122g[ACc(A-a)-1]2">h=v212g[ACc(A−a)−1]2

The area of obstruction when the throttle valve is fully open is 1.6063e-5m2">1.6063e−5m2

The coefficient of contraction is calculated as 0.93066">0.93066

The coefficient loss due to obstruction is 0.02388">0.02388

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

Equation 13

HLtotal=(fLD+0.3+0.4+0.02388)∗v122g">HLtotal=(fLD+0.3+0.4+0.02388)∗v212g

Length of pipe =1.7732396565229E-01m">=1.7732396565229E−01m

Diameter of pipe =0.0171742m">=0.0171742m

Friction factor for the circular pipe is given as f=64Re">f=64Re

Reynolds number for the flow is given as Re=ρ∗v1∗Dμ">Re=ρ∗v1∗Dμ

Rewriting Equation 13 

HLtotal=((64∗μρ∗v1∗D)LD+0.3+0.4+0.02388)∗v122g">HLtotal=((64∗μρ∗v1∗D)LD+0.3+0.4+0.02388)∗v212g

Rewriting Equation 4 

Ptotal=Ps2+ρgh2+ρv222+((64∗μρ∗v1∗D)LD+0.3+0.4+0.02388)∗v122g">Ptotal=Ps2+ρgh2+ρv222+((64∗μρ∗v1∗D)LD+0.3+0.4+0.02388)∗v212g

Equation 14

150000=100000+1.177∗9.81∗0.119998+1.177v222+((64∗1.846e-051.177∗v1∗0.017142)0.1773230.0171742+0.3+0.4+0.02388)∗v122∗9.81">150000=100000+1.177∗9.81∗0.119998+1.177v222+((64∗1.846e−051.177∗v1∗0.017142)0.1773230.0171742+0.3+0.4+0.02388)∗v212∗9.81

where ρ">ρ

at 300K and 1atm pressure =1.177kgm3">=1.177kgm3

Dynamic viscosity of the air at 300K 1.846e-05kgm-sec">1.846e−05kgm−sec

Ptotal=150000Pa">Ptotal=150000Pa

Ps1=100000Pa">Ps1=100000Pa

As such there are two unknowns and one equation which forced us to assume that let v1=v2">v1=v2

Solving the above equation will give the two values of v1">v1

0.6253v12+0.0204v1-49998.61">0.6253v21+0.0204v1−49998.61

Therefore v1=282.755,-282.787">v1=282.755,−282.787

. Rejecting the negative value of v1">v1

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=vc">M=vc

where v">v

is the velocity of the fluid in msec-1">msec−1

c">c

is the speed of the sound in msec-1">msec−1 which is 346.3msec-1">346.3msec−1at 300K">300K

M=282.755346.3">M=282.755346.3

M=0.816">M=0.816

While all flows are compressible, the flows are usually treated as incompressible when (M<0.3">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=RTV-b-aT0.5V(V+b)">P=RTV−b−aT0.5V(V+b)

 where a=0.42748∗R2Tc2.5Pc">a=0.42748∗R2T2.5cPc

b=0.08664∗RTcPc">b=0.08664∗RTcPc

R=287.05JKgK">R=287.05JKgK

Critical Temperature Tc">Tc=132.63K">132.63K

Critical Pressure Pc">Pc=3.7858MPa">=3.7858MPa 

Absolute temperature = 300K">300K

Therefore a=1.88626">a=1.88626

 b=0.00087128471">b=0.00087128471

 Putting the values of a and b in Equation 16

150000=86115V-0.000871-1.8862617.32V(V+0.000871)">150000=86115V−0.000871−1.8862617.32V(V+0.000871)

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

2598000V3-1491511.8V2-1299.19V+0.00164293">2598000V3−1491511.8V2−1299.19V+0.00164293

Therefore V=0.57497,-0.00087,0">V=0.57497,−0.00087,0

Thus the molar volume of the air at the inlet calculated from the Redlich Kwong equation is 0.57497m3kg">0.57497m3kg

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 ρ1=1molarvolume">ρ1=1molarvolume

=10.57497">=10.57497=1.739kgm3">=1.739kgm3

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

Pstatic=Ptotal{1+γ-12∗M2}-(γ-1γ)">Pstatic=Ptotal{1+γ−12∗M2}−(γ−1γ)

where Mach number, M=va">M=va

and a=γRT">a=√γRT

Rewriting the Equation 9

P1static=Ptotal{1+γ-12∗v2γRT}-(γγ-1)">P1static=Ptotal{1+γ−12∗v2γRT}−(γγ−1)

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.

(140000150000)-13.5={1+1.4-12∗V21.4⋅287.05⋅300}">(140000150000)−13.5={1+1.4−12∗V21.4⋅287.05⋅300}

 Velocity at the inlet =109.54 msec">msec

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

Therefore HL=446.07m">HL=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 Re=ρ∗V∗Dμ">Re=ρ∗V∗Dμ

where ρ">ρ

at 300K and 1atm pressure =1.177kgm3">=1.177kgm3

D">D

is the inlet diameter of the pipe, which is 0.0171742m">0.0171742m

 Dynamic viscosity of the air at 300K 1.846e-05kgm-sec">1.846e−05kgm−sec

Therefore Re=119954.848474">Re=119954.848474

For an internal flow if Re≥4000">Re≥4000

then the flow is turbulent.

 Calculating Hydrodynamic Entrance Length Le">Le

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

Le=10∗D">Le=10∗D

Le=0.171742m">Le=0.171742m

Using the value of Le">Le

to calculate the local Reynolds number which is 1199548.48474">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 (δ">δ) calculated from one-seventh power law combined with empirical data for turbulent flow through smooth pipes is given by 

Equation 20 

δx=0.38Rex0.2">δx=0.38Re0.2x

Therefore δ=0.00397060377m">δ=0.00397060377m

Calculating Skin friction Coefficient

Equation 21

Cf=0.0599796510.2">Cf=0.0599796510.2

Cf=0.003589">Cf=0.003589

Calculating Wall Shear Stress

Equation 22

τw=0.5∗ρ∗U2∗Cf">τw=0.5∗ρ∗U2∗Cf

τw=25.3434Nm-2">τw=25.3434Nm−2

Calculating Frictional Velocity 

Equation 23

Ut=τwρ">Ut=√τwρ

Ut=4.6402msec-1">Ut=4.6402msec−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 Rex,cr=5∗105">Rex,cr=5∗105

and much lesser than the Rex,transition=50∗105">Rex,transition=50∗105

, 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

δx=4.91Rex">δx=4.91√Rex

δ=0.000235m">δ=0.000235m

Calculating Skin friction Coefficient for a laminar flat plate 

Equation 26

Cflam=0.664Rex">Cflam=0.664√Rex

Cflam=0.001985">Cflam=0.001985

Calculating Wall Shear Stress for a laminar flat plate 

Equation 27

τw=0.5∗ρ∗U2∗Cf">τw=0.5∗ρ∗U2∗Cf

τw=14.01Nm-2">τw=14.01Nm−2

Calculating Frictional Velocity for a laminar flat plate 

Equation 28

Ut=τwρ">Ut=√τwρ

Ut=3.45msec-1">Ut=3.45msec−1

Summarizing all the values

Calculating first cell height using fixed embedment from Converge-CFD

Sample Calculation

Equation 29

Embed=base-grid-size2Embed-scale">Embed=base−grid−size2Embed−scale

Keeping Embed scale to 1, for element size 4e-3">4e−3

, the value of Embed is `

Embed=0.0042">Embed=0.0042

Embed=0.002m">Embed=0.002m

First cell height =0.002m">=0.002m

While dealing with the control volume, y+">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=2yp">△y=2yp

Therefore,

yp=0.0022">yp=0.0022

yp=0.001m">yp=0.001m

y+=ρ∗Ut∗ypμ">y+=ρ∗Ut∗ypμ

For base grid size 4mm">4mm

y+=295.856">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.