Week 4.1: Project Steady state simulation of flow over a throttle body
1. Introduction
A huge piping network employed in different applications such as thermal power generation, petroleum or transport of water to commercial and residential complexes often involve different pipe fittings, valves and elbows due to various constraints in otherwise supplying typically fluid mediums across straight channels. The additional components also end up being responsible for controlling and regulating the supply and provide great assistance in conducting repair works in the damaged portions, however, on the downside and especially due to the disruption caused to an otherwise regular flow, losses are also incurred such as flow separation and intense mixing which impacts the energy requirements to conduct the flow in the first place.
The aim of this challenge was to perform three-dimensional, steady-state simulations using Realizable k-`epsilon` (RKE) turbulence model for an elbow-shaped pipe domain consisting of a throttle body. The analysis involved computing results for different grid sizes that would outline total cell count, velocity magnitude profiles, mass flow rate and pressure drop values across the domain, mesh verification and comparison against analytical formulations, identification of separation regions as well as mesh quality attributes like skewness, non-orthogonality, and aspect ratio.
2. Method
The typical approach adopted to performing numerical simulations in CONVERGE CFD involves the following steps:
- Creating a new folder dedicated to running the concerned project inside a user-preferred directory.
- Setting up the relevant case file(s) in Converge Studio v3.0 that involves creation of the fluid volume domain, setting up the required application type, choosing material and its state (gas, liquid, parcel or solid) for which simulation is to be performed along with the requisite properties, defining simulation parameters such as solver types, monitor variables, run time variables, and convergence tolerances, defining region(s), boundary conditions and initial conditions, specifying the physical models whose equations are intended to be solved, determining and allocating the base grid size values as well as any other refinements, and lastly selecting the variables for post-processing at the required time intervals (in terms of cycles for steady-state) along with any hidden parameters or UDF selection as part of advanced parameters.
- Saving the case file(s) and exporting the input files into the folder created for running the concerned project which typically includes 'boundary.in' for the various boundary conditions, 'material_state.dat' where material_state points to either 'gas', 'liquid', 'parcel' or 'solid', 'initialize.in' specifying initial conditions, 'inputs.in' defining the various input conditions and file writing intervals, 'monitor_steady_state.in' to assign monitor variables for the steady state simulation, 'turbulence.in' containing information about the turbulence model selected and its associated constants, 'post.in' for the different post-processing variables, 'solver.in' describing the solver chosen and its equations, and 'surface.dat' containing datasets about the surface while other input files with respect to the state of the material, thermal properties and species may also be part of the exporting process and should be ensured that they are selected for the simulation to run.
- Running the simulation in a Cygwin64 terminal window.
- Post-converting 3D output data by navigating to the output folder created within the same project folder during the simulation run or from the 'Post-Processing 3D' button in Converge Studio v3.0.
The simulation run and post-converting processes were executed using Linux-type commands in a Cygwin64 terminal window as shown below:
i) Change directory to the project containing the relevant case file(s):
cd "PATH_TO_CONCERNED_PROJECT_FOLDER"
ii) Confirm the contents of the project for the exported input files:
ls
iii) Performing simulation run with the converge executable application on multiple processors (in this case, 4) using mpiexec executable file in the background (by adding ampersand operator) along with posting output into a logfile:
mpiexec.exe -n 4 converge.exe restricted </dev/null> logfile &
iv) Checking for simulation status with regards to completion or errors:
tail logfile
v) Post-converting the 3D output data files using the 'post_convert_30_msmpi_64' utility on multiple processors (in this case, 4) using mpiexec executable file from within the output folder of the project containing the exported input files:
cd output
mpiexec.exe -n 4 post_convert_30_msmpi_64
2.1 Modelling
The computational domain for steady state simulation of air flow movement through an elbow-shaped pipe over a throttle body was modelled in Converge Studio v3.0. The length of the domain was provided as 0.2091378m in the direction of flow from the inlet to the outlet and 0.0171742m as the inlet pipe diameter. The basic element entity in CONVERGE is denoted by a triangle which would be visible by clicking on any of the surfaces of the recently drawn model and in order to demarcate boundaries, the corresponding triangles need to be grouped together. The grouping and creation of new boundaries process began with clicking on the 'Boundary' button under the 'Geometry' dock which was followed by clicking on the '+' symbol from which the pop-out message would read as 'Create a new boundary.' leading to the 'Boundary definition' window to open in which the option 'Create multiple boundaries' was checked with number of boundaries as 3, in this case, followed by lastly grouping the triangles that constitute a given surface and naming them accordingly, such as Inlet, Outlet, Elbow Pipe Walls and Throttle Body as illustrated in Fig. 1. Moreover, the normal from each of the triangular surface is required to point towards the fluid domain which was ensured by first visualising the same using the 'Normal toggle' button accessible in the 'View Toolbar' and depending on the direction along which the normal is pointing, a decision was taken by clicking on the 'Normal' tab within the 'Transform' button under the 'Geometry' dock and selecting the suitable option from one of the three 'Transform Normal Options' followed by clicking on the 'Apply' button. At this stage and for the given numerical analysis, the modelling requirements were achieved, however, to validate if the model was clean i.e., it contains no intersecting edges, overlapping triangles or any other form of geometrical constraints, a check was performed by clicking on the 'Find' button in the 'Diagnosis' dock to observe for any anomalous 'Problem Triangles' and rectify the same.
Fig. 1: Fluid Volume Model and Boundary Conditions regions in three-dimensional space
2.2 Meshing
The meshing approach adopted for the given computational domain was developed to maintain the meshing quality attributes viz., aspect ratio, skewness, and non-orthogonality as outlined in the CONVERGE 3.0 Manual provided in the following table:
|
Aspect Ratio |
1 (Ideal Value) |
||
|
Skewness |
0° (Ideal Value) |
< 70° (Acceptable) |
`>=` 90° (Likely to cause simulation failure) |
|
Non-Orthogonality |
0° (Ideal Value) |
< 70° (Acceptable) |
`>=` 90° (Likely to cause simulation failure) |
Table 1: Mesh Quality Attributes as defined in the CONVERGE 3.0 Manual [1]
The reason for pursuing the above mentioned mesh quality attributes, as specified in Table 1, can be ascertained on the fact that whether the discretized domain was likely to create convergence issues during the solution process and also to identify areas in the mesh where the statistics can be improved based on the underlying physics. In Section 2.3.4 Mesh Refinement approach, detailed discussion about the different meshes that were prepared and their mesh qualities has been highlighted.
2.3 CFD Simulation Setup, Boundary Conditions, Monitoring Variables, Determination of maximum possible flow velocity and Solution
The simulations on all meshes was performed using the Reynolds Averaged Navier-Stokes (RANS) models, in this case, the k-`epsilon` turbulence models. The use of k-`epsilon` turbulence models and more specifically the RKE turbulence model is attributed to the fact that the Reynolds number, based on the given flow conditions, lies in the turbulence flow region as well as representing high Reynolds number attached flows for which the solution variables, like the velocity profile, was modelled for the log-law layer using the standard wall function approach such that the first grid cell would lie in the range of 30<`y^+`<300 in order to yield satisfactory results. A brief description into the working of the RKE turbulence model is provided below:
2.3.1 RKE turbulence model
The RKE turbulence model proposed by Shih [2] is considered as an upgrade to the standard k-`epsilon` turbulence model while sharing the similarities of consisting two transport equations, one for the turbulent kinetic energy (k) to evaluate the energy in turbulence and the other for the turbulent dissipation (`epsilon`) to calculate the rate of dissipation of turbulent kinetic energy (k), however notable differences are present as mentioned below:
- The turbulent viscosity term is formulated differently in the RKE turbulence model.
- The dissipation rate `epsilon` in the transport equation has been modified and developed from the exact equation for transport of mean-square vorticity fluctuation.
The RKE turbulence model is able to satisfy some of the mathematical constraints of the Reynolds's stresses such as the positive values of normal stresses when the eddy viscosity definition is used in the Boussinesq relationship and Schwarz inequality for shear stresses even when the main strain rate is large which otherwise would be negative. These formulations in the RKE turbulence model have made them more popular to use over the other k-`epsilon` turbulence models especially when flow features include the presence of strong streamline curvature, vortices, rotation, boundary layers under the influence of strong and adverse pressure gradients, separation and recirculation.
2.3.1.1 Setup of the RKE turbulence model:
To use the RKE turbulence model in Converge Studio v3.0, we navigate to Physical Models option in the Case Setup tree followed by checking the 'Turbulence modeling' box to enable inclusion of turbulence terms in the solution of governing fluid equations and then selecting the Reynolds Averaged Navier-Stokes (RANS) item in the 'Turbulence modeling' window available under the Physical Models option with 'Turbulence model' drop-down menu set to 'Realizable k-`epsilon`' and 'Near Wall treatment' retained at its default selection of 'Standard wall function' while all other variables/parameters and constants remained unaltered with respect to their original values.
2.3.2 Boundary Conditions and Monitoring Variables
The different boundary conditions for simulation of air flow movement through an elbow-shaped pipe over a throttle body is described in the following table:
|
Names |
Inlet |
Throttle_Body_wall |
Elbow_pipe_walls |
Outlet |
|
Type |
INFLOW |
WALL |
WALL |
OUTFLOW |
|
Pressure Boundary Condition |
Specified Value (DI), Total pressure |
-- |
-- |
Specified Value (DI) |
|
Pressure Magnitude (Pa) |
150000 |
-- |
-- |
100000 |
|
Velocity Boundary Condition |
Zero normal gradient (NE) |
-- |
-- |
Zero normal gradient (NE) |
|
Temperature Boundary Condition |
Specified Value (DI) |
Specified Value (DI) |
Specified Value (DI) |
-- |
|
Temperature Backflow |
-- |
-- |
-- |
Specified Value (DI) |
|
Temperature Magnitude (K) |
300 |
300 |
300 |
300 |
|
Species Boundary Condition |
Specified Value (DI) |
-- |
-- |
Specified Value (DI) |
|
Species Name: `O_2` |
0.23 |
-- |
-- |
0.23 |
|
Species Name: `N_2` |
0.77 |
-- |
-- |
0.77 |
|
Wall motion type |
-- |
Stationary |
Stationary |
-- |
|
Surface Movement |
-- |
Fixed, Law of wall |
Fixed, Law of wall |
-- |
|
Turbulent Kinetic Energy (tke) |
Intensity, 0.033 fraction [3] |
Zero normal gradient (NE) |
Zero normal gradient (NE) |
-- |
|
Turbulent Dissipation (eps) |
Length scale, 0.001202m [4] |
Wall model |
Wall model |
-- |
|
Turbulent Kinetic Energy (tke) Backflow |
-- |
-- |
-- |
Intensity, 0.033 fraction [3] |
|
Turbulent Dissipation (eps) Backflow |
-- |
-- |
-- |
Length scale,0.001202m [4] |
|
Near wall treatment |
-- |
Standard wall function |
Standard wall function |
-- |
|
Law of wall roughness parameters; Absolute roughness |
-- |
0.0m |
0.0m |
-- |
|
Law of wall roughness parameters; Roughness constant |
-- |
0.5 |
0.5 |
-- |
Table 2: Boundary Conditions for the elbow-shaped pipe domain consisting of a throttle body defined using the 'Boundary' sub-option of the 'Boundary Conditions' under the 'Case Setup' tree
For the given problem, the simulation was chosen as 'Gas simulation' enabled in the 'Materials' option in the 'Case Setup' tree with 'Predefined mixtures' set to 'Air' and 'Species (not in mech.dat)' option enabled to utilise the mass fraction values of `O_2` and `N_2` as they represent the major constituting species in air by mass fraction. Since, 'Gas simulation' option was checked in the 'Materials' option, an equation of state was required to be defined that formulates a coupling relationship between the different state variables such as density and physical conditions such as pressure and temperature to depict the state of matter. In the given numerical analysis, the Redlich-Kwong (RK) cubic equation of state was set to in the 'Equation of state' drop-down menu with its default values of critical temperature (`T_c`) and critical pressure (`p_c`) retained at '133.0 K' and '3770000.0 Pa' respectively as it is a simpler form of representation in comparison to the Redlich-Kwong-Soave and Peng-Robinson set of equations but more accurate than the ideal gas law which works best for gases at low pressures, however, at higher pressures tends to report huge deviations in its predictions as it assumes that the molecular volume is very small and that neither attraction or repulsion exists between the molecules which is not true in the case of real gases. Additionally, the real gas properties would behave as a function of temperature which was again set to as the default option in the 'Gas simulation' window. The values of 'Turbulent Prandtl number' and 'Turbulent Schmidt number' were left unchanged at their default values of 0.9 and 0.78 respectively.
During the solution process, a monitor variable was created to observe the behaviour of flow-field variables viz., velocity, density, pressure and mass flow rate, by selecting 'AVG_VELOCITY' from the list of 'Monitor Variables' available by default with 'Target Type' and 'Target' set to 'Boundary' and 'Outlet' respectively and file name' as 'mass_avg_flow.out' from the 'File Name' drop-down menu. The settings of the remaining parameters for the reporting of the monitor plots included entering values of 'Minimum number of cycles to be executed for the steady-state solver' and 'Sample size, cyc' such that the given solution would run for at least 5000 cycles regardless of whether a state of convergence would have been achieved prior to it and 1000 number of samples to decide if the concerned variable has reached steady-state condition with 'Monitor start delay, cyc' having a value of 0 before clicking on the 'Validate' and 'OK' buttons.
2.3.3 Determination of flow velocity across the elbow-shaped pipe domain consisting of a throttle body, Reynolds number, Boundary layer thickness and First Cell Height
The given numerical problem was defined based on pressure boundary conditions and in order to employ the type of gas flow solver as incompressible or compressible, it was decided to conduct a preliminary assessment of the flow field leading to approximated values of velocity across the domain. Assuming an incompressible flow, the relation between static pressure and total pressure at the inlet boundary using INFLOW boundary condition can be derived as shown below [1]:
`P_1(static) = P_(TOTAL) - 0.5*rho*U_i^2 - rho*g*h`; where (Eq. 1)
`P_(static)` is the static pressure,
`P_(TOTAL)` is the total pressure, in this case, the pressure magnitude specified at the inlet boundary using INFLOW boundary condition,
`rho` is Density of fluid,
`U_i` is the velocity at the inlet boundary,
`g` is acceleration due to gravity, and
`h` is the elevation above a reference point or plane
Similarly, the relation between static pressure and total pressure at the outlet boundary using OUTFLOW boundary condition can be expressed as follows [1]:
`P_2(static) = (P_(bound) + predist*P_(fluid))/(1 + predist)`; where (Eq. 2)
`P_(bound)` is the boundary pressure, in this case, the pressure magnitude specified at the outlet boundary using OUTFLOW boundary condition,
`predist` is weighting factor to dampen reflecting pressure waves within the domain and subsequently improve convergence rates for steady-state simulations, and
`P_(fluid)` is the boundary cell pressure,
For the given numerical problem, the `predist` parameter is unchanged from its default value of 0 thereby implying that the static pressure is the same as the outlet boundary pressure designated using OUTFLOW boundary condition. Since, the flow was assumed to be incompressible across the domain, steady flow energy equation was utilised which in simplified form can be mathematically written as:
`P_1 + 0.5*rho*V_1^2 + rho*g*h_1 = P_2 + 0.5*rho*V_2^2 + rho*g*h_2 + H_L`; where (Eq. 3)
`P_1` and `P_2` are the pressure energy terms in terms of static pressure at inlet and outlet boundaries respectively,
`V_1` and `V_2` are velocities in the kinetic energy per unit volume terms (`0.5*rho*V^2`) at inlet and outlet boundaries respectively,
`h_1` and `h_2` are elevations above a reference point or plane in the potential energy per unit volume terms (`rho*g*h`) at inlet and outlet boundaries respectively, and
`H_L` is total head loss, in this case the sum of major loss across the elbow-shaped pipe domain and minor losses such as the elbow bend, throttle body (similar to a butterfly valve) and throttle body as an obstruction wherein major loss through the elbow bend is given by:
`H_L(MAJOR) = f*(L/D)*(V_1^2/(2*g))`; where (Eq. 4)
`L` is the length of the pipe determined as 0.2091378m,
`D` is pipe diameter evaluated as 0.0171742m, and
`f` is the friction factor, in this case assuming laminar flow conditions for a circular pipe as `f = 64/(Re_D)` or `f = 64*mu/(V_1*D*rho)` as velocity `V_1` is unknown. Similarly, the minor losses are expressed using the following relations:
- Minor loss for the elbow bend:-
`H_L(elbow_(bend)) = K_(eb)*(V_1^2/(2*g))`; where (Eq. 5)
`K_(eb)` is the loss coefficient for the elbow bend which for the given configuration of smooth and rounded 90° bend equates to a value of 0.29 [5].
- Minor loss for the throttle body as a rotating or butterfly valve:-
`H_L(TB) = K_(tb)*(V_1^2/(2*g))`; where (Eq. 6)
`K_(tb)` is the loss coefficient for the throttle body whose function is similar to that of a rotating or butterfly valve and in the given configuration, the position of the throttle body allows for uninterrupted movement of fluid thereby indicating a maximum loss coefficient value of 0.25 [5].
- Minor loss for the throttle body as an obstruction within the elbow domain:-
`H_L(TB_(obs)) = K_(tb_(obs))*(V_1^2/(2*g))`; where (Eq. 7)
`K_(tb_(obs))` is the loss coefficient for the throttle body as an obstruction within the elbow domain and is expressed as the ratio of the area of the pipe to the product of the coefficient of contraction and difference of the areas of the pipe and obstruction i.e., `((A/(C_c*(A-a)))-1)^2`. For the open position of the throttle body, area of obstruction was determined to be 1.51204E-05`m^2`, coefficient of contraction as the ratio of the difference of the areas of the pipe and obstruction to the area of the pipe as 0.934728917 and consequently, the loss coefficient for the throttle body as an obstruction within the elbow domain as 0.02089.
Therefore, total head loss can be expressed as:
`H_L(TOTAL) = H_L(MAJOR) + H_L(elbow_(bend)) + H_L(TB) + H_L(TB_(obs))`
OR
`H_L(TOTAL) = {[(64*mu)/(V*D*rho)*(L/D)] + 0.29 + 0.25 + 0.020890035}*V^2/(2*g)` (Eq. 8)
Substituting `P_1` as static pressure in velocity `U_i` term from (Eq. 1), which is equivalent to `V_1` term of (Eq. 3), `P_2` as static pressure from (Eq. 2) which is equivalent to the outlet boundary pressure designated using OUTFLOW boundary condition and `H_L(TOTAL)` from (Eq. 8) as `H_L` into (Eq. 3), we get
`P_(TOTAL) - 0.5*rho*V_1^2 + 0.5*rho*V_1^2 + rho*g*h_1 - rho*g*h_1 = P_(bound) + 0.5*rho*V_2^2 + rho*g*h_2 + {[(64*mu)/(V_1*D*rho)*(L/D)] + 0.29 + 0.25 + 0.020890035}*V_1^2/(2*g)`
OR
`P_(TOTAL) = P_(bound) + 0.5*rho*V_2^2 + rho*g*h_2 + {[(64*mu)/(V_1*D*rho)*(L/D)] + 0.29 + 0.25 + 0.020890035}*V_1^2/(2*g)` (Eq. 9)
Applying the continuity equation for a system of one-inlet, one-outlet, open throttle-body position and steady and incompressible flow, as shown below, it was observed that as there is negligible change in cross-sectional area of flow through the elbow bend, the inlet and exit velocities were assumed to be approximately the same such that `V_1 = V_2 = V`. Simplifying (Eq. 9) using known values of `P_(TOTAL)` as 150000 Pa, `P_(bound)` as 100000 Pa, `g` as 9.81 `m/s^2`, `h_2` as 0.09773m, `rho` and `mu` of air at 1atm or 101235 Pa and 300K (or 26.85°C) as 1.1766 `(kg)/(m^3)` and `1.857*10^-5 (kg)/(m.s)` respectively [6], `V` throughout the domain was determined as 284.6636m/s. The reason for determining the flow velocities across the flow domain was to estimate the Mach number (or local Mach number) and accordingly justify the use of compressible or incompressible flow equations with respect to the given numerical problem. Mach number, by definition, is a dimensionless quantity which represents the ratio of local flow velocity, in this case local air flow across the elbow-shaped domain, to speed of sound in the medium or the speed of sound in air. As the flow velocity through the given elbow-shaped pipe domain with the throttle body in fully opened position was evaluated as 284.6636m/s and speed of sound in air is approximately equal to 347.336m/s for the same temperature conditions but 10% relative humidity [7], the Mach number resulted in an approximate value of 0.81956 which is much higher than the typical range of 0.2 to 0.3 for which compressibility effects are required to be considered. Therefore, the above justification of incompressible flow does not hold true and subsequently a new approach with regards to a compressible flow was taken such that the relation between static pressure and total pressure at the inlet boundary using INFLOW boundary condition was now derived as shown below [1]:
`P_1(static) = P_(TOTAL)*(1+(gamma-1)/(2)*(U_i^2)/(gamma*R*T))^-((gamma)/(gamma-1)) - rho*g*h`; where (Eq. 10)
`gamma` is the ratio of specific heats, equal to 1.4 [8],
`R` is the ideal gas constant, equal to 287.05 J/kg-K [9], and
`T` is temperature of fluid
As proved from initial assessment, the assumption that flow was incompressible across the domain was incorrect, the energy equation now had to be rewritten for compressible flow while considering it to be isothermal as temperature remained constant as shown below:
`(P_1)/(rho_1*g)log_eP_1 + V_1^2/(2*g) + z_1 = (P_2)/(rho_2*g)log_eP_2 + V_2^2/(2*g) + z_2 + H_L`; where (Eq. 11)
`z_1` and `z_2` are elevations above a reference point or plane, similar to `h_1` and `h_2`, at inlet and outlet boundaries respectively,
Continuity equation for a compressible flow indicates that mass flow rate across the domain is conserved and can be mathematically written as:
`rho_1*A_1*V_1 = rho_2*A_2*V_2`; where (Eq. 12)
`rho_1` and `rho_2` are the density values at inlet and outlet boundaries respectively
Rewriting the above expression in terms of velocities at the inlet and outlet boundaries i.e., (`V_1`) and (`V_2`), we get
`V_2/V_1 = (rho_1*A_1)/(rho_2*A_2)` (Eq. 13)
Substituting the values of `A_1` and `A_2`, in this case of equal magnitude, we get
`V_2/V_1 = (rho_1/rho_2)` (Eq. 14)
Also, as mentioned before that the process is isothermal, the relation between pressure and density can be expressed as shown below:
`P_1/rho_1 = P_2/rho_2` (Eq. 15)
On a closer look at (Eq. 2), (Eq. 10), (Eq. 11), (Eq. 14) and (Eq. 15), it was noticed that there are 5 equations and 6 unknowns and therefore, it was decided to allocate a guess pressure value of 125000Pa for static pressure at the inlet boundary which not only adhered to being less than the total pressure at the same boundary but also corresponded to a value exactly between the total pressure at the inlet and the static pressure at the outlet resulting in an approximated value of velocity as 179.537m/s. However, as total pressure was unknown at the outlet boundary, value of density at the inlet boundary was required to be determined using the cubic equation of state as expressed below [1]:
`P = (R*T)/(v-b) - (a)/(v^2+u*b*v+w*b^2)` (Eq. 16)
Since, the selected cubic equation of state was Redlich-Kwong (RK), the values of parameters `u`, `w`, `b` and `a` were 1, 0, `beta_(rk)*v_c` and `alpha_(rk)*(R^2*T_c^2.5)/(p_c)`. Moreover, the values of `v_c` in `b` and `a` were `(R*T_c)/p_c`, `beta_(rk)` in `b` was 0.08664 and `alpha_(rk)` in `a` was 0.42748. The notation `v_c` refers to the critical volume. Substituting, the values of ideal gas constant (R), critical temperature (`T_c`), critical pressure (`p_c`), coefficients `alpha_(rk)` and `beta_(rk)` in parameters `a` and `b`, and ideal gas constant (R) and temperature (`T`) in the first term on the right hand side of (Eq. 16), we get its simplified version at total pressure of 150000 Pa as shown below:
`150000 = (86115)/(v-0.0008774) - (1.906)/(v^2+0.0008774*v)` (Eq. 17)
It should be noted that the terms ideal gas constant (R) and critical pressure (`p_c`) in parameter `a` were multiplied by a factor of 1e-3 for dimensional consistency. Solving the cubic equation above (Eq. 17), we get the value of term `v` or the molar volume as 0.573704`m^3/(kg)`. Inverting the value obtained for the molar volume, the density of air at given conditions of temperature and total pressure using the Redlich-Kwong (RK) cubic equation of state was determined as 1.7391`(kg)/m^3`. Utilising this value of density at the inlet, and dynamic viscosity of air at 1atm or 101235 Pa and 300K (or 26.85°C) as `1.857*10^-5 (kg)/(m.s)` [6], a corresponding Reynolds number, based on inlet pipe diameter, was estimated as 288764.0548 (see Eq. 18) indicating the flow was in the turbulent regime as it is well above the starting limit for turbulent internal flow of 4000.
Reynolds number (based on circular cross-section inlet diameter):
`Re_D = (U*D*rho)/mu`; where (Eq. 18)
`mu` is the dynamic viscosity of fluid
Substituting the value of density and more specifically at the inlet boundary in (Eq. 14) and (Eq. 15) and velocity at the inlet boundary in (Eq. 14), expressions for velocity and pressure at the outlet boundary are formulated in terms of density at the same boundary, i.e.,
`V_2 = (312.2328/rho_2)` (Eq. 19)
`P_2 = 86251.5094*rho_2` (Eq. 20)
Since, the flow is no longer incompressible and laminar, a new relationship or value to define the friction factor was required for the compressible and turbulent flow regime which in this case was evaluated by taking reference from the Moody diagram [10] for the curve representing a smooth pipe, as was for the given case, indicating a value of about 0.016. Recalculating the total head loss from (Eq. 8) in terms of the new approximated value of velocity at the inlet and major head loss while retaining the minor losses as expressed in (Eq. 5), (Eq. 6) and (Eq. 7), we get
`H_L(TOTAL) = H_L(MAJOR) + H_L(elbow_(bend)) + H_L(TB) + H_L(TB_(obs))`
OR
`H_L(TOTAL) = (f*(L/D) + 0.29 + 0.25 + 0.020890035)*(V_1^2/(2*g)) = 1241.517m` (Eq. 21)
Bringing all pressure energy terms and kinetic energy terms on the left and right hand sides of (Eq. 11) respectively, the equation can be rewritten as:
`(P_1)/(rho_1*g)log_eP_1 - (P_2)/(rho_2*g)log_eP_2 = V_2^2/(2*g) - V_1^2/(2*g) + 0.09773 + 1241.517` (Eq. 22)
Since, the process is isothermal, implementing (Eq. 15) into (Eq. 21), we get:
`(P_1)/(rho_1*g)log_eP_1 - (P_1)/(rho_1*g)log_eP_2 = V_2^2/(2*g) - V_1^2/(2*g) + 0.09773 + 1241.517`
OR
`(P_1)/(rho_1*g)[log_e(P_1/P_2)] = (V_2^2/(2*g)) - (V_1^2/(2*g)) + 0.09773 + 1241.517` (Eq. 23)
Substituting the values of terms `P_1` (total pressure), `rho_1` and `V_1` in their numerical form, and `V_2` and `P_2` (total pressure) in terms of `rho_2` as seen in (Eq. 19) and (Eq. 20) above into (Eq. 23), the value of `rho_2` obtained was 1.309`(kg)/m^3`. Reutilising the value of `rho_2` in (Eq. 19) and (Eq. 20), the values of `V_2` and `P_2` (total pressure) evaluated are 238.528m/s and 112903.2258Pa respectively. It can be observed that the values of velocity and pressure (static and/or total, wherever applicable) that were initially guessed and then calculated at the inlet and outlet boundaries indicate a realistic solution thereby implying the flow to be physically possible to travel across the domain. It should be noted, as has been reiterated before, that these calculations do not represent converged state values but is an initial assessment of the flow field in order to establish whether the flow is incompressible or compressible and subsequently deduce approximated values of Reynolds number, boundary layer thickness and first cell height.
2.3.4 Mesh Refinement approach
In this numerical simulation, a baseline mesh for the given computational domain was created with default element size of 0.04m and without fixed embedding, illustrated in Fig. 2, to get an initial assessment of the flow field. Thereafter, successive mesh refinements were carried out by utilising values for different mesh sizes in conjunction with the empirically estimated turbulent boundary layer thickness, as seen in Table 4 and Table 5, and then corrected to adjust according to the approximated value of `EMBED SCALE`, recorded in Table 6, that were then part of each individual case setup and was entered into corresponding columns representing the 'Base grid size' section of the 'Base grid' sub-option in the 'Grid Control' settings under the 'Case Setup' tree to compare the results obtained for the velocity magnitude profiles, mass flow rate and pressure drop across the entire domain.
Boundary layer thickness, in case of turbulent flow, is usually dependent on experiments and/or empirical relations as it difficult to predict the same because of the intense mixing of different layers involved within it. Therefore, the mathematical expression used for determining the boundary layer thickness is obtained from one-seventh-power law combined with empirical data for turbulent flow through smooth pipes and is expressed below as:
`(delta/x) = 0.38/(Re_x)^(1/5)`; where (Eq. 24)
`delta` is the boundary layer thickness, and
`Re_x` is the local Reynolds number at a point `x` in the direction of flow within the smooth pipe, in this case the elbow-shaped pipe domain
The point `x` depends on the flow characteristics developing within the elbow-shaped pipe domain such that a presence of strong inviscid core occurs from the entrance of the domain to the location wherein the boundary layers from the opposite walls merge to form a fully developed flow due to the dominating viscous effects which is also referred to as the entrance length. The entrance length, in this case, is also the point at which the boundary layer thickness will have the maximum height before which the velocity reaches free-stream velocity value. Entrance length, in the case of turbulent flow, when the domain length is greater than 10 times the inlet boundary diameter can be expressed mathematically as:
`L_e approx 10D` (Eq. 25)
Substituting the value of `D` in (Eq. 25), the entrance length for the given elbow-shaped pipe domain was evaluated to a value of 0.171742m and as such lies towards the exit of the straight section of the pipe after the bend. Similarly, utilising the value of `L_e` in place of `D` in (Eq. 18), the local Reynolds number was calculated to be 2887640.548. Therefore, the turbulent boundary layer thickness based on the entrance length definition and the Reynolds number at the same point i.e., `x = L_e`, in (Eq. 24) was determined as 0.003331m. Similarly, the fully-open throttle body represents a flat plate, assumed to be smooth, was also considered with regards to boundary layer development for which the characteristic length would be represented by the length of the throttle body in the direction of flow resulting in a Reynolds number of 275217.613. This value of Reynolds number based on the characteristic length for the throttle body indicates that the flow is in the transitional regime but the nature of the boundary layer is laminar as the estimated Reynolds number does not exceed the critical Reynolds number of `5*10^(5)` with respect to a smooth flat plate geometrical considerations such that the thickness of the boundary layer would result in a value of 0.0001532m from the resulting expression (Eq. 26) as seen below:
`(delta/L) = 4.91/(Re_L)^(1/2)` (Eq. 26)
In order to capture the boundary layer development for high Reynolds number attached flows, wall modelling approach of using a wall function was chosen such that the first grid cell was in the range of 30 < `y^+` < 300 based on the RKE turbulence model selected. To perform calculations for `y^+` value, the values of wall shear stress and frictional velocity we required to be calculated as shown below:
Wall Shear Stress:
`tau_w=(1/2)⋅C_f⋅rho⋅(U_∞)^2`; where (Eq. 27)
`tau_w` is Wall shear stress,
`C_f` is Skin Friction coefficient,
`U_∞` is the freestream velocity, in this case the calculated value of inlet boundary velocity
Local skin friction coefficient (`C_f`) based on a semi-empirical expression can be related to the Reynolds number for two cases viz., turbulent boundary layer for pipe flow (in this case, elbow shaped pipe domain) and laminar boundary layer for flat plate (i.e., throttle body), as shown below:
`C_f(turb_(eb))=0.0594/(Re_x)^(1/5)` (Eq. 28)
`C_f(laminar_(tb))=1.33/(Re_L)^(1/2)` (Eq. 29)
Frictional velocity (`U_tau`) can be computed mathematically as:
`U_tau=sqrt(tau_w/rho)` (Eq. 30)
For local Reynolds number with regards to the characteristic length of inlet diameter for the elbow-shaped pipe domain and characteristic length of throttle body as 2887640.548 and 275217.613 respectively with the previously mentioned values of dynamic viscosity, density and free-stream velocity, the skin-friction coefficients (`C_f(turb_(eb))` and `C_f(laminar_(tb))`), wall shear stress (`tau_w`) and friction velocity (`U_tau`) values were calculated using (Eq. 28), (Eq. 29), (Eq. 27) and (Eq. 30) respectively as shown in the following table:
|
Characteristic Length (L) |
Local Reynolds Number |
Skin-friction Coefficient (`C_f`) |
Wall Shear Stress (`tau_w`)[kg m^-1 s^-2] |
Friction Velocity (`U_tau`)[m s^-1] |
|
Inlet Pipe diameter (D) |
2887640.548 |
`3.0316*10^(-3)` |
84.972 |
6.989 |
|
Throttle Body diameter (d) |
275217.613 |
`2.535*10^(-3)` |
71.053 |
6.392 |
Table 3: Skin-friction coefficient, wall shear stress and friction velocity values for the governing flow conditions of the elbow-shaped pipe domain and throttle body
Calculation of cell-centroid size (`y_P`), First Cell Height (`Delta_y`) and embed scale with respect to Base Grid Size:
`y^+` is defined as the non-dimensional wall distance for a wall bounded flow that can be expressed mathematically as:
`y^+=(ρ*U_tau*y_P)/mu`; where (Eq. 31)
`y_P` is distance from the cell-centroid 'P' to the wall
Simplifying (Eq. 31) to determine `y_P`, we get;
`y_P=(y^+*mu)/(U_tau*rho)` (Eq. 32)
Since, `y^+` is computed at the cell-centroid, the first cell height (`Delta_y`) can accordingly be mathematically expressed as:
`Delta_y = 2*y_P` (Eq. 33)
The approach adopted for evaluating the distance from the cell-centroid 'P' to the wall (`y_P`) and correspondingly the first cell height (`Delta_y`) was to target a `y^+` value of 290 which based on the required transition between it and the base mesh size led to the different mesh sizes as reported in Table 4 along with the `EMBED SCALE` as is demonstrated below and tabulated in Table 5 and Table 6.
From the CONVERGE 3.0 Manual [1], the mathematical formulation relating the embed size with the base grid size based on embed scale can be given by:
`embed = (base grid)/(2^(EMBED SCALE))`; where (Eq. 34)
`base grid` is the Base Grid Size defined for every element across the domain,
`EMBED SCALE` is the scaling parameter by which the base grid size is to be refined either at the selected boundary, geometrical shape or region, and
`embed` is the scaled down or refined grid size, in this case also equal to the First Cell Height i.e., `Delta_y`,
The missing variable in (Eq. 34) is `EMBED SCALE` and therefore equating for the same variable was required with regards to the given grid sizes.
|
Mesh |
dx(m) |
dy(m) |
dz(m) |
|
Baseline |
0.004 |
0.004 |
0.004 |
|
Mesh 1 |
0.00124 |
0.00124 |
0.00124 |
|
Mesh 2 |
0.00115 |
0.00115 |
0.00115 |
|
Mesh 3 |
0.00106 |
0.00106 |
0.00106 |
Table 4: Mesh sizes representing spatial grid discretization values across each of the three-directional axes
|
Mesh |
`EMBED SCALE` |
Embed Layers |
Initial First Cell Height (`Delta_y`)[m] |
dy (`base grid`)[m] |
element size % increase (w.r.t `Delta_y`) |
Initial Target `y^+` (`y_P`) |
|
1 |
0.485 |
3.759 |
0.0008861 |
0.00124 |
~40 |
~290 |
|
2 |
0.3761 |
3.759 |
0.0008861 |
0.00115 |
~30 |
~290 |
|
3 |
0.258 |
3.759 |
0.0008861 |
0.00106 |
~20 |
~290 |
Table 5: Base Grid and its transition w.r.t initial First Cell Height (`Delta_y`), `EMBED SCALE`, Target `y^+`(`y_P`) and Embed Layers for different meshes of the elbow-shaped pipe domain
It should be noted that the values of `EMBED SCALE` are required to be a positive integer and therefore the values indicated in Table 5 above were rounded off and consequently new first cell height (`Delta_y`), target `y^+` and expansion ratio values were evaluated as is reported in Table 6 below. Similarly, the number of embed layers was initially determined based on the ratio of the calculated boundary layer thickness to the initial first cell height (`Delta_y`) and then recalculated based on the new first cell height (`Delta_y`).
|
Mesh |
`EMBED SCALE` |
Embed Layers |
New First Cell Height (`Delta_y`)[m] |
dy (`base grid`)[m] |
expansion ratio (`Delta_y` to dy) |
New Target `y^+` (`y_P`) |
|
1 |
1 |
5.372 ~ 6 |
0.00062 |
0.00124 |
~2 |
~202.9 |
|
2 |
1 |
5.793 ~ 6 |
0.000575 |
0.00115 |
~2 |
~188.176 |
|
3 |
1 |
6.28 ~ 7 |
0.00053 |
0.00106 |
~2 |
~173.45 |
Table 6: Base Grid and its transition w.r.t First Cell Height (`Delta_y`) for rounded-off `EMBED SCALE`, new Target `y^+` (`y_P`) and Embed Layers values for different meshes of the elbow-shaped pipe domain
Since, the Reynolds number calculated for the given domain using the throttle body diameter in the direction of flow as the characteristic length was found to be in the transitional regime but the boundary layer to be laminar, it was decided to refine the grid size on the throttle body using the different base grid sizes mentioned in Table 6 above leading to the `EMBED SCALE` values on it as reported in Table 7 below.
|
Mesh |
`EMBED SCALE` throttle body |
Laminar Boundary Layer Thickness (`Delta_y_(LAM)`)[m] |
dy (`base grid`)[m] |
|
1 |
3.0168 ~ 3 |
0.0001532 |
0.00124 |
|
2 |
2.9081 ~ 3 |
0.0001532 |
0.00115 |
|
3 |
2.7906 ~ 3 |
0.0001532 |
0.00106 |
Table 7: `EMBED SCALE` for the throttle body at different meshes of the elbow-shaped pipe domain
Baseline Mesh:
Fig. 2: Baseline Mesh (two-dimensional sliced view)
Mesh 1:
Fig. 3: Mesh 1 (two-dimensional sliced view)
Mesh 2:
Fig. 4: Mesh 2 (two-dimensional sliced view)
Mesh 3:
Fig. 5: Mesh 3 (two-dimensional sliced view)
Mesh and time statistics, peak velocity magnitude values near the outlet, and pressure drop across the domain for different meshes:
|
Mesh |
Cells |
Run Time [s] |
First Cell Height (`Delta_y`)[m] |
`y^+` (`y_P`) (Line 1) |
Peak Velocity Mag. [m/s] (Line 1) |
% var |
Pressure drop [Pa] |
% var |
|
Baseline |
702 | 86.26 | 0.00120904 | 1400.437 | 252.016 |
-- |
20988 |
-- |
|
1 |
207002 |
8886.7714 |
0.00060904 |
220.142 |
234.699 |
-6.87 |
23646 |
12.664 |
|
2 |
241324 |
10614.2204 |
0.00033404 |
200.698 |
235.609 |
0.386 |
23518 |
-0.541 |
|
3 |
284566 |
13859.2074 |
0.00013903 |
243.325 |
235.633 |
0.0102 |
23589 |
0.3019 |
Table 8: Mesh and time statistics, peak velocity magnitude values at Line 1, pressure drop values across the domain, `y^+` at (`y_P`) for different meshes at Line 1 and percentage variation in values of quantities of interest between subsequent meshes
Mesh quality attributes for the different meshes:
|
Mesh |
Max Aspect Ratio |
Max Skewness |
Max Non-Orthogonality |
expansion ratio (first to 2nd layer of boundary cells) |
expansion ratio (last row of boundary cells to dy) |
|
Baseline |
2.8842 | 1.17578e-38 | 67.0893 | -- | 3.308 |
|
1 |
3.2996 |
1.17578e-38 |
86.6679 |
1.018 |
~2 |
|
2 |
3.4323 |
1.17578e-38 |
87.9037 |
1.721 |
~2 |
|
3 |
4.1017 |
1.17578e-38 |
88.1113 |
3.812 |
~2 |
Table 9: Mesh quality attributes depicting the maximum values of aspect ratio, skewness and non-orthogonality for the different meshes of the elbow-shaped pipe domain
From the meshing statistics in Table 8, it can be seen that the prediction of peak velocity magnitude values at Line 1 and pressure drop values between mesh number 2 (or mesh 2) with 241324 cells and mesh 3 with 284566 cells are quite close to each other with less than 0.02% and 0.4% in variations respectively and as such a mesh independent result can be said to have been achieved with Mesh 2 configuration in terms of the computational time and resources. The %var is calculated as the ratio of the difference of the predicted values in two consecutive meshes to the predicted value of the previous mesh e.g., ((Peak Velocity Magnitude for Mesh 1 - Peak Velocity Magnitude for Mesh 2)/Peak Velocity Magnitude for Mesh 1). It can also be observed that the `y^+` values predicted for all the meshes at the first layer cell-centroid position are in the log-law region, however, the differences between the initially calculated first cell height and the same computed by CONVERGE were quite apart and growing with every mesh size and this may be attributed to the constantly changing domain structure due to the elbow bend in the direction of flow such that the cells cut at the boundary, using the cartesian cut-cell approach, attempts to fit the fixed embedding region with the required and defined number of cells before reaching the base grid size. Additionally, mesh quality attributes of all the different meshes prepared in Table 9 indicates that they lie in the range of acceptable value for skewness, as previously stated in Table 1, but deteriorates for aspect ratio and non-orthogonality that was particularly observed in the region between the throttle body and the elbow pipe walls as the change in the cell structures were more erratic and discontinuous which can be ascertained to be a consequence of the different mesh size refinements defined for the two boundaries.
2.4 Post-Processing
The extended part of the solution process using Converge Studio v3.0 included selecting the variables of interest whose quantitative and qualitative aspects would highlight the behaviour of the flow, in this case within the given fluid domain, under different circumstances and they were chosen by clicking on the 'Post variable selection' sub-option of the 'Output/Post-Processing' setting under the 'Case Setup' tree and checking the boxes 'Density', 'Mesh Quality', 'Pressure', 'Temperature' and 'Velocity' in the 'Typical Variables' section and 'Region ID' in the 'Geometry/Location' section of the 'Typical Variables/Geometry' tab, 'Turbulent Kinetic Energy', 'Turbulent Viscosity' and '`y^+`' in the 'Turbulence' section of the 'Combustion/Turbulence' tab, and adding 'Mass Fraction' under 'Species' section of the 'Species/Passive' tab for `O_2` and `N_2` species with respect to the concerned numerical analysis. Visualising the different flow variable quantities in two-dimensional space in ParaView required the creation of a slice of the given fluid domain by using the "Slice" tool in the 'Commons' toolbar typically located above the 'Pipeline Browser' window based on default layout settings and choosing 'Z Normal' under 'Plane Parameters' section to align the plane normal to the z-axis direction such that the entire domain was viewed in the direction of flow followed by unchecking the 'Show Plane' and 'Triangulate the slice' boxes for disabling the viewing of the plane borders and production of triangles in the output respectively. Plotting the velocity magnitude profiles near the outlet just after the analytically calculated fully-developed flow condition of the given domain was performed by using ParaView's in-built "plot over line" tool. The "plot over line" tool was accessed by clicking on 'Filters' present in the data analysis group, followed by entering the coordinates as seen in Table 10, and then clicking on the "Apply" button to view the said profile. In order to compare numerical simulation results of velocity magnitude profiles against its analytical counterparts, following expressions were created, one at a time, before plotting them at the coordinates as seen in Table 10.
Expressions:
r: arc_length-(D_inlet/2)
ratio_radius_center_to_pipe: 1-abs((r)/(D_inlet/2))
friction_factor_1: (100*`Re_D`)^(-1/4) [11]
n: 7 [12]
U_max_friction_factor_1: U*(1+(1.33*sqrt(friction_factor_1)))
U_max_power_law: U/(((2*(n^2))/((n+1)*((2*n)+1)))
Power_law_velocity_profile_U_max_friction_factor_1: U_max_friction_factor_1*(ratio_radius_center_to_pipe^(1/n))
Power_law_velocity_profile_U_max_power_law: U_max_power_law*(ratio_radius_center_to_pipe^(1/n))
Friction factor utilised in the above list of expressions relates to the Blasius correlations on account of ease in calculating it as it contains no term for pipe roughness making it effectively valid for only smooth pipes or in this case, a smooth elbow-shaped pipe domain. There is, however, a drawback in using the Blasius formulation for friction factor as it is known to hold true until Reynolds number of 100000, which for the given case had clearly crossed the limit. Therefore, other formulations for friction factor were also considered such as the Kármán–Prandtl resistance equation that is also supported by Nikuradse’s experiments [13] and Unified Friction Formulation by Brkić and Praks [14], and their corresponding profiles can be viewed under Section 3 Results. It was ensured that the last time-step size of the simulation run was selected to view the velocity and `u^+` vs `y^+` profile results. Additionally, velocity (U) and pressure (p) field plots were created to view the development of the flow field characteristics across the fluid domain, and glyph and tube based on stream tracers were added for identifying the high and low velocity regions. Legends, titles, and x and y labels were added for better understanding of the graphical representations.
|
Line 1 |
x |
y |
z |
|
Point 1 |
-0.00331697 |
0.01921 |
0.020174 |
|
Point 2 |
-0.00331697 |
0.0020424 |
0.020174 |
Table 10: Line Coordinates for viewing the velocity magnitude and `u^+` vs `y^+` profiles
To save an animation plot of the different field plots, the option 'Save Animation' was accessed under the 'File' menu and clicked on followed by choosing the appropriate settings under the different sections such as 'Size and Scaling', 'Coloring', 'Compression Options' and 'Animation Options' to name a few, assigning an appropriate name to the animation file, in this case by the variable whose field plot is being indicated before clicking on the 'OK' button.
2.5 Convergence Level
In the given numerical analysis, a steady-state simulation was performed using a Density-based solver, to account for the compressibility effects, of type modified Pressure Implicit with Splitting of Operator (PISO) which is an iterative pressure-velocity coupling algorithm involving a predictor step for predicting solution of the momentum equation from which a pressure equation is obtained and solved whose correction is then consequently reapplied to the momentum equation such that this process continues until the desired level of accuracy is achieved. At the end of the momentum predictor and first corrector step, remaining transport equations are solved in series. In other words, momentum and pressure are solved at the first PISO iteration which enables calculation for velocity terms required by other transport equations and at the end of each iteration in the PISO loop, convergence criterion is checked for which in the case of compressible flow requires the density correction error term to be less than the value of PISO tolerance which remains unchanged for the given numerical simulation at its default value of 1e-03. The selection of solver types for the different equations within 'Solver parameter [Steady-state]' sub-option in the 'Simulation Parameters' settings under the 'Case Setup' tree included assigning SOR to Momentum, Density, Energy, Turbulent Kinetic Energy (tke) and Turbulent dissipation rate (eps) equations, and 'CONVERGE BiCGSTAB' to the Pressure equation. The SOR solver type is a generally preferred algorithm owing to its advantage in being computationally less expensive than the different forms of biconjugate gradient stabilized (BiCGSTAB) algorithms or HYBRID (wherever applicable) except for the pressure equation unless the number of iterations for this equation (pressure) is less than 15 on average making it suitable for cases that contains a compressible flow with smaller value of time-steps. However, for the concerned numerical simulation, the default 'CONVERGE BiCGSTAB' algorithm was retained as it was difficult to estimate whether or not the pressure equation would always result in less than 15 iterations on average with subsequent mesh refinements. BiCGSTAB algorithm works best in scenarios for which simulations approach convergence at a slower rate making it more suitable for multiphase flow simulations like volume of fluid (VOF). Also, as previously mentioned in Section 2.3.2 Boundary Conditions and Monitoring Points, the behaviour of velocity, density, pressure and mass flow rate at the outlet were being monitored as the quantities of interest to observe the variation in values over successive cycles. The solver types, convergence tolerance limits, and minimum and maximum iterations defined for the equations corresponding to the RKE turbulence model equations and the monitor variable plots corresponding to the different mesh sizes are shown in the table and figures below:
|
Equations |
Solver Type |
Convergence Tolerance |
Minimum Iterations |
Maximum Iterations |
Preconditioner |
|
Momentum |
SOR |
1e-05 |
2 |
50 |
None |
|
Pressure |
CONVERGE BiCGSTAB |
1e-08 |
2 |
300 |
None |
|
Density |
SOR |
1e-04 |
2 |
20 |
None |
|
Energy |
SOR |
1e-04 |
2 |
20 |
None |
|
Turbulent Kinetic Energy (tke) |
SOR (pre-defined) |
1e-03 |
2 |
50 |
None (pre-defined) |
|
Turbulent dissipation rate (eps) |
SOR (pre-defined) |
1e-03 |
2 |
50 |
None (pre-defined) |
Table 11: Solver Types, Convergence Tolerance Values, Minimum and Maximum Iterations and Preconditioners of the different equations used in the RKE turbulence model
Fig. 6: Plots of monitor variables velocity, density, pressure and mass flow rate for steady-state simulation using RKE turbulence model for baseline mesh
Fig. 7: Plots of monitor variables velocity, density, pressure and mass flow rate for steady-state simulation using RKE turbulence model for mesh 1
Fig. 8: Plots of monitor variables velocity, density, pressure and mass flow rate for steady-state simulation using RKE turbulence model for mesh 2
Fig. 9: Plots of monitor variables velocity, density, pressure and mass flow rate for steady-state simulation using RKE turbulence model for mesh 3
From Fig. 6, Fig. 7, Fig. 8 and Fig. 9, it was noticed that the monitor variables show a highly stable output with minimum variations well before they reach the end of total number of cycles specified in the 'End time' parameter under the 'General' tab of 'Simulation Parameters' settings under the 'Case Setup' tree as 15000 indicating that a steady-state solution was obtained. Moreover, the converged state velocity at the inlet for all the different mesh refinements from Mesh 1 to Mesh 3 was recorded, on average, to be around 180.5m/s and as such resembles quite closely with the velocity value estimated from initial assessments at the same boundary to be at 179.537m/s while at the outlet, the difference was considerably larger with final state numerically simulated values at around 216m/s as opposed to the analytical estimation at the same boundary to be at 238.528m/s mentioned in Section 2.3.3 Determination of flow velocity across the elbow-shaped pipe domain consisting of a throttle body, Reynolds number, Boundary layer thickness and First Cell Height.
3. Results
3.1 Velocity magnitude and pressure gradient values using field plots
Fig. 10 and Fig. 11 below shows the velocity magnitude and pressure gradients across the elbow-shaped pipe domain consisting of a throttle body for different mesh sizes in two-dimensional space using field plots at the slice region created by orienting the plane normal to the flow direction as was specified in Section 2.4 Post-Processing. It can be observed that while the Baseline mesh was unable to predict the gradients for both velocity and pressure regions mainly in the vicinity of the throttle body, the differences in the same became indistinguishable with increase in the spatial resolution as they were consequently able to offer improved visualisation of the development of viscous effects, pressure drop and flow separation regions across the domain.
Fig. 10: Velocity Magnitude field plots across the domain for different mesh sizes in two-dimensional space
Fig. 11: Pressure field plots across the domain for different mesh sizes in two-dimensional space
3.2 Velocity magnitude values using line plots
Fig.12 below shows the velocity magnitude profile plots at Line 1 wherein it can be observed that all the refined meshes have almost identical values in the core of the flow region with nearly similar peak magnitudes, as was also reported in Table 8, but the agreement with the analytical power-law velocity profiles were not essentially reasonable and the reason may be established on the grounds that the fluid movement undergoes several changes during the course of its travel including at the elbow bend and flow over the fully opened throttle body leading to separation before reattaching further downstream wherein the turbulent boundary layer thickness is regained. Besides, the numerical results for the velocity magnitude profiles do not necessarily qualify for validation against its analytical counterparts as the location at which the comparisons have been drawn ideally represents the fully-developed flow conditions for a straight pipe with a constant cross-sectional area whereas in the given scenario the conditions indicate an elbow-bend pipe with a throttle body as an obstruction and therefore, the position of the hydrodynamic length (or point in the flow at which fully-developed flow conditions are achieved) would typically be located further downstream, in this case, beyond the given extents at which the outlet of the domain has been demarcated.
Fig. 12: Plot of comparison of velocity profiles between numerical and analytical (power-law) at 0.2068m from inlet of elbow-shaped pipe domain with a throttle body at Line 1
3.3 Comparison of law of the wall and logarithmic-law velocity profiles
Since, the gradients of the velocity field are varying by large magnitudes near the wall on account of the developing boundary layer, it was necessary to ensure that this region was properly resolved which in the given case was performed using the RKE turbulence model and a standard wall function approach such that the first grid cell was in the log-law layer or within the range of 30 < `y^+` < 300. The layer closest to the wall (or the viscous sub-layer) was suggested by Prandtl [15] to be the zone wherein the velocity profile is linear such that the relation `u^+ = y^+` is obeyed, also confirmed by Nikuradse's experiments [16], usually upto `y^+ = 5` for smooth surfaces and similarly the logarithmic law-overlap layer, proposed by Von Karman [17], is the region wherein the velocity behaves as being directly proportional to the logarithm of the distance such that `u^+ = (1/kappa)ln y^+ + B`, typically in the range of 30 < `y^+` < 300, such that values of `kappa` and 'B' were experimentally evaluated as 0.42 [18] and 5.5 for smooth wall surfaces which in Converge are also known as 'Von Karman's constant' and 'Law of the wall parameter' respectively. A third layer exists as a region between the viscous sub-layer and the log-law layer, called as the Buffer layer, where both the viscous and turbulence effects occur but neither velocity profile from the other two regions are followed or are accurate to represent the ongoing physical effects with maximum deviation between the two velocity profiles occurring at the intercept position of around `y^+ = 11`, similar to the observation in Fig. 13 below. Investigations were carried out to describe a universal turbulent velocity behaviour in an attempt to fit Nikuradse's experimental data [16] near the wall such as the three-layer Karman model [19] represented by three separate algebraic equations for the three different layers in the inner region of the velocity profile, Reichardt's [20] and Spalding's [21] law of the wall models defined with single algebraic equations for `y^+>=0` with the latter showing comparable conformity with the given simulation results in the viscous sub-layer and log-law layer. It should be realised that the velocity profiles as seen in Fig. 13 corresponding to the different layers in the inner region and the universal turbulent velocity behaviour were designed for incompressible fluid flow with values of its associated quantities belonging to the averaged and converged state conditions corresponding to last time-step size, however, in the given numerical simulation and based on Reynolds number of the inlet pipe diameter of circular cross-section, the flow was found to be compressible and turbulent and therefore, they should be at best considered to be suggestive and not the real representation, including the value mentioned in the plot at the cell-centroid of the first grid cell as opposed to its actual value mentioned in Table 7, and while different formulations do exist for the compressible law of the wall, sufficient information about the wall temperatures, wall density and wall thermal conductivity among many others were not known at the time of performing the given numerical simulation.
Fig. 13: Plot of comparison of law of the wall and logarithmic-law velocity profiles for Mesh 2 at 0.2068m from inlet of the elbow-shaped pipe domain with a throttle body and `y^+` value at distance from the cell-centroid 'P' to the wall (`y_P`)
3.4 Separation and Recirculation regions
As observed in Fig. 14 and Fig. 15 below, velocity at the inner pipe wall of the flow approaching the elbow bend was found to be of a significantly greater magnitude, due to the presence of a curvature, than the same quantity at the outer pipe wall subsequently causing separation. The fluid flow also undergoes further separation as it traverses its way over the fully opened throttle body before the transition of the high-velocity takes place from the inner pipe wall to the outer pipe wall in the downstream direction. Additionally, the region in front of the throttle body becomes a high pressure zone as a result of the fluid impinging directly on its surface with the attached flow advancing around it causing local acceleration in the crevices between the elbow pipe walls and the throttle body followed by detachment and generation of a low pressure region behind the throttle body. The region between the separated flow and ensuing reattachment was identified as recirculating flow region, in this case, located in front and behind the throttle body, and in the curved section of inner pipe walls.
Fig. 14: Elbow-shaped pipe domain consisting of a throttle body coloured using local pressure field while streamlines in the form of tubes oriented in the direction of flow as per the local velocity field magnitude for Mesh 2
Fig. 15: Animation of elbow-shaped pipe domain consisting of a throttle body coloured using local pressure field while streamlines in the form of glyphs (arrows) oriented in the direction of flow as per the local velocity field magnitude for Mesh 2
4. Discussion
Air enters at the inlet of the given cross-section of the elbow-shaped pipe domain consisting of a throttle body, containing a mixture of two major contributing gas species `O_2` and `N_2` by mass fraction, as turbulent flow. The development of eddy structures in terms of movement and intensity is considerably at lower rates near the elbow pipe walls in the straight section of the domain from the inlet because of the "no slip" condition such that the turbulent boundary layer formed experiences a steep velocity gradient close to the wall but changes at a slower rate in its core region. This sudden change in velocity magnitude values close to the wall enhances the magnitude of wall shear stress creating a more fuller shape of the fully developed velocity profile as opposed to parabolic if the same flow were to be laminar. As the wall-attached flow approaches towards the elbow bend of the given domain, velocity at the inner pipe wall begins to accelerate to account for the curvature effects leading to generation of swirling secondary flows and thereby separation. This flow separation is further supplemented by the presence of the throttle body as it acts as an obstruction leading to two separate paths for the fluid flow. The first flow path adopted is over the throttle body with the fluid particles initially striking at the front portion of the structure, producing a zone of high pressure, with the sudden change in cross-sectional area causing flow to be separated before giving rise to a boundary layer over it, in this case laminar with regards to the Reynolds number based on the throttle body diameter as the characteristic length, followed by the sudden change in cross-sectional area at the trailing edge causing the flow to decelerate such that a zone of unfavourable or adverse pressure gradient is created leading to flow separation. Therefore, two regions of recirculating flow are as such established both over the leading and trailing edges of the throttle body. The second flow path is around the throttle body with a similar impingement of fluid particles at the front of its leading edge but in this case the flow remains attached over it followed by subsequent rapid advancement towards the narrow gap between the throttle body and elbow pipe walls, again due to the curvature effects, and then separation in the trailing edge forming yet another recirculating region as it attempts to rejoin the first flow path which in this case can be considered as the major portion of streamline flow because of the area it covers over the throttle body within the confined domain. The impact on the turbulent boundary layer thickness can therefore be described to be increasing in the straight section of the model from the inlet due to the increasing viscous effects, with decrease in the length of the same boundary layer as it approaches the elbow bend before beginning to thicken again in the downstream section and because of the added effects of the pressure drop, velocity of the air heading towards the outlet increases. The recirculation zones recognised as low-pressure regions lead to shear layer formations that eventually goes on to reattach, depending on the length of the domain, as a consequence of development of the unstable and large-scale turbulent structures beyond which a state of equilibrium for the turbulent boundary layer profile is achieved with no reverse flow.
5. Conclusion
The numerical analysis for flow within an elbow-shaped pipe domain consisting of a throttle body under steady-state conditions using RKE turbulence model in Converge Studio v3.0 and executed using the Cygwin64 terminal was well understood. It was realised that based on the wall modelling approach of using a standard wall function and assigning appropriate `y^+` values, the turbulent or log-law region was well captured for all mesh sizes, although as an approximation to the incompressible law of the wall model, along with visualisation of the flow separation in the curved region of inner pipe wall and throttle body and an eventual reattachment further downstream. The assumption of an initial static pressure value at the inlet leading to determination of flow-field variables like velocity, density and pressure also stresses the importance of assigning correct boundary and initial conditions as the final state or converged values especially at the inlet were understandably close allowing for selection of an appropriate solver and empirically estimating the boundary layer profile growth. The CFD model, however, requires further verification in terms of incorporating additional information about the compressibility effects and validation as part of an extended scope wherein both the length of the inlet and outlet sections of the given domain would be sufficiently long enough to allow for fully-developed flow conditions to occur based on the calculated Reynolds number.
References
[1] CONVERGE 3.0 Manual. © 2020 Convergent Science. All rights reserved. July 15, 2020.
[2] Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z., and Zhu, J., 1995, “A New - Eddy-Viscosity Model for High Reynolds Number Turbulent Flows - Model Development and Validation.”, Computers Fluids, 24(3):227-238
[3] Basse, N.T. (2017), "Turbulence intensity and the friction factor for smooth- and rough-wall pipe flow", Fluids, vol. 2, 30.
[4] Flórez Orrego; et al. (2012). "Experimental and CFD study of a single phase cone-shaped helical coiled heat exchanger: an empirical correlation". Proceedings of ECOS 2012 – The 25th International Conference on Efficiency, Cost, Optimization, Simulation and Environmental Impact of Energy Systems, June 26–29, 2012, Perugia, Italy. ISBN 978-88-6655-322-9.
[5] SUEZ's degremont® water handbook, minor losses in the pipelines, fittings, valves... for water, www.suezwaterhandbook.com, viewed 20 June 2021, <https://www.suezwaterhandbook.com/formulas-and-tools/formulary/hydraulics/minor-losses-in-the-pipelines-fittings-valves-for-water>
[6] Engineers Edge, LLC, Viscosity of air, Dynamic and Kinematic, www.engineersedge.com, viewed 20 June 2021, <https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm#:~:text=The%20viscosity%20of%20air%20depends,2%20%2Fs%20or%2014.8%20cSt.>
[7] The Engineering ToolBox, Speed of Sound in Air, www.engineeringtoolbox.com, viewed 20 June 2021, <https://www.engineeringtoolbox.com/air-speed-sound-d_603.html>
[8] The Engineering ToolBox, Ratio of Specific Heat of Gases, www.engineeringtoolbox.com, viewed 20 June 2021, <https://www.engineeringtoolbox.com/specific-heat-ratio-d_608.html>
[9] The Engineering ToolBox, Universal and Individual Gas Constants, www.engineeringtoolbox.com, viewed 20 June 2021, <https://www.engineeringtoolbox.com/individual-universal-gas-constant-d_588.html>
[10] Moody, L.F. (1994). "Friction Factors for Pipe Flow". Transactions of the ASME. 66 (8): 671-684.
[11] BLASIUS, P. R. H. 1913. Das Aehnlichkeitsgesetz bei Reibungsvorgangen in Flüssigkeiten. Forschungsheft 131, 1-41
[12] Hinze, J. O., Turbulence, 2nd Ed., McGraw-Hill, New York, 1975
[13] Nikuradse, J. “Stroemungsgesetze in rauhen Rohren.” Ver. Dtsch. Ing. Forsch., 361, 1933
[14] D. Brkić and P. Praks, Unified friction formulation from laminar to fully rough turbulent flow, Applied Sciences (2018). 8 (11), article 2036 URL <https://www.scipedia.com/public/Brkic_2018a>
[15] L. Prandtl, "Uber die ausgebildete Turbulenz," Z. fur angew. Math. Mech., vol. 5,1925, p. 136
[16] J. Nikuradse, "Gesetzmiissigkeiten der turbulenten Stromung in glatten Rohren," Forschung a.d. Geb. Ing., No. 356, 1932.
[17] Von Karman, T. H. (1931). Mechanical similitude and turbulence. NACA Technical Memorandum 611.
[18] Kármán, Theodore von. 1930. “Mechanische Ähnlichkeit und Turbulenz.” Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch–Physikalische Klasse, 58–76.
[19] T. von Karman, "The Analogy Between Fluid Friction and Heat Transfer," Trans. ASME, vol. 61, 1939, pp. 705-710.
[20] H. Reichardt, "Vollstandige Darstellung der turbulenten Geschwindigkeitsverteilung in glatten Leitungen," Z. angew. Math. Mech., vol. 31, 1951, pp. 208-219
[21] Spalding, D.B. Heat Transfer to a Turbulent Stream from a Surface with a Step-wise Discontinuity in Wall Temperature. ASME Conference. Intern. Dev. in Heat Trans. Boulder, Colo. Part 11:439-466. 1961