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

Aim: Steady-state simulation of flow over a throttle body

Objective: 

Set up and run the steady-state simulation for flow over the throttle body.

Post-process the results and show pressure and velocity contours.

Show the mesh (i.e surface with edges)

Show the plots for pressure, velocity, mass flow rate, and total cell count.

Create an animation in which its throttle movement should be visible.

Introduction

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_s1+frac{1}{2}rhov_1^2+rhogh=P_(t otal)`

where `P_(s1)` is the static pressure at the inlet.

`frac{rhov_1^2}{2}` is the dynamic pressure at the inlet.

`rhogh` 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_(t otal)-frac{1}{2}rhov_1^2`

Similarly at the outlet of the pipe,

Equation 3 

`P_(s2)+frac{1}{2}rhov_2^2=P_(t otal)`

`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{rhov_1^2}{2}+rhogh_1+h_(pump)=P_(s2)+frac{rhov_2^2}{2}+rhogh_2+h_(turb i n e)+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_(Lt otal)=H_(Lmajo r)+H_(Lmi no r)`

`H_(Lmajo r)=ffrac{L}{D}frac{v_1^2}{<a href="https://skill-lync.com/electronics-engineering-courses/introduction-cellular-communications-2g-3g" target="_blank">2g</a>}`

`H_(Lmi no r)=K_Lfrac{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^@elbow` `K_(Lelbow)` is `0.3`

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

Equation 7

`H_(Lthrot tl e)=K_(Lthrot tl e)frac{v_1^2}{2g}`

For throttle fully open `K_(Lthr o t t l e)=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{(v_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_cv_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_(Lt otal)=(ffrac{L}{D}+0.3+0.4+0.02388)**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**v_1**D}{mu}`

Rewriting Equation 13 

`H_(Lt otal)=((frac{64**mu}{rho**v_1**D})frac{L}{D}+0.3+0.4+0.02388)**frac{v_1^2}{2g}`

Rewriting Equation 4 

`P_(t otal)=P_(s2)+rhogh_2+frac{rhov_2^2}{2}+((frac{64**mu}{rho**v_1**D})frac{L}{D}+0.3+0.4+0.02388)**frac{v_1^2}{2g}`

Equation 14

`150000=100000+1.177**9.81**0.119998+frac{1.177v_2^2}{2}+((frac{64**1.846e-05}{1.177**v_1**0.017142})frac{0.177323}{0.0171742}+0.3+0.4+0.02388)**frac{v_1^2}{2**9.81}`

where `rho` at 300K and 1atm pressure `=1.177frac{kg}{m^3}` 

Dynamic viscosity of the air at 300K `1.846e-05frac{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 `msec^-1`

`c` is the speed of the sound in `msec^-1` which is `346.3msec^-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 (`Mlt0.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.5V(V+b)}`

where `a=0.42748**frac{R^2T_c^2.5}{P_c}`

`b=0.08664**frac{RT_c}{P_c}`

`R=287.05frac{J}{KgK}`

Critical Temperature `T_c`=`132.63K`

Critical Pressure `P_c``=3.7858 MPa`

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.57497frac{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}{molar volume}``=frac{1}{0.57497}``=1.739frac{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_(t otal){1+frac{gamma-1}{2}**M^2}^-(frac{gamma-1}{gamma})`

where Mach number, `M=frac{v}{a}` and `a=sqrt(gammaRT)`

Rewriting the Equation 9

`P_(1static)=P_(t otal){1+frac{gamma-1}{2}**frac{v^2}{gammaRT}}^-(frac{gamma}{gamma-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.

`(frac{140000}{150000})^(frac{-1}{3.5})={1+frac{1.4-1}{2}**frac{V^2}{1.4*287.05*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**V**D}{mu}`

where `rho` at 300K and 1atm pressure `=1.177frac{kg}{m^3}` 

`D` is the inlet diameter of the pipe, which is `0.0171742m`

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

Therefore `R_e=119954.848474`

For an internal flow if `R_egt=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**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}{Re_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**rho**U^2**C_f`

`tau_w=25.3434Nm^-2`

Calculating Frictional Velocity 

Equation 23

`U_t=sqrt(frac{tau_w}{rho})`

`U_t=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 `R_(ex,cr)=5**10^5` and much lesser than the `R_(ex,transition)=50**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(Re_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**rho**U^2**C_f`

`tau_w=14.01Nm^-2`

Calculating Frictional Velocity for a laminar flat plate 

Equation 28

`U_t=sqrt(frac{tau_w}{rho})`

`U_t=3.45msec^-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-scal e)`

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 

`triangley=2y_p`

Therefore,

`y_p=frac{0.002}{2}`

`y_p=0.001m`

`y+=frac{rho**U_t**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 

Animation file

Velocity

Pressure

Mass flow rate

 

Line plots

Velocity 

Pressure

Mass flow rate 

Total cell count

Mesh-1

Contour plots

Velocity

Pressure

Mass flow rate 

Animation file

Velocity

Pressure

Mass flow rate

 

 

Line plots

Velocity 

Pressure

Mass flow rate 

Total cell count

 

Mesh-2

Contour plots

Velocity

Pressure

Mass flow rate 

Animation file

Velocity

Pressure

Mass flow rate

Line plots

Velocity

Pressure

Mass flow rate

Total cell count

 

Mesh-3

Contour plots

Velocity

Pressure

Mass flow rate 

Animation file

Velocity

Pressure

Mass flow rate

 

Line plots

Velocity

Pressure

Mass flow rate 

Total cell count

Mesh-4

Contour plots

Velocity

Pressure

Mass flow rate 

Animation file

Velocity

Pressure

Mass flow rate

 

Line plots

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.  

 

 

 

Source 

Generalized Cubic Equation of state [https://www.sciencedirect.com/topics/chemistry/cubic-equation-of-state]

Engineering Toolbox [https://www.engineeringtoolbox.com/]

Redlich Kwong Equation of state [https://en.wikipedia.org/wiki/Redlich%E2%80%93Kwong_equation_of_state]

Air-Thermophysical properties[https://www.engineeringtoolbox.com/air-properties-d_156.html#:~:text=Critical%20temperature%3A%20132.63%20K%20%3D%20%2D,3%20%3D%2018.89%20lbm%2Fft]

Embedding in Converge [https://skill-lync.com/knowledgebase/embedding-in-converge-2]

Throttle [https://en.wikipedia.org/wiki/Throttle]

Fluid Mechanics Fundamentals and Applications- Yunus  A. Cengel, John M. Cimbala