Exhaust Port

1. Introduction

The exchange of thermal energy involving a solid body and a fluid volume at the shared interfaces can be described as a conjugate heat transfer (CHT) phenomena. The CHT analysis of any given system considers the physical processes and corresponding governing equations involved for the individual subdomains separately and then correlating their solutions by matching the required conditions at the interfaces to get the distribution of temperature and heat flux in a coupled manner thereby eliminating the need of a heat transfer coefficient. A CHT analysis thus offers advantages of monitoring both fluid flow and heat transfer changes as well as real time prediction of wall heat transfer coefficient, temperature and velocity profiles across the domain. CHT analysis over recent years have aided significantly in improving the cooling performance of electronic equipments, such as design of heat sinks, to designing effective heat exchangers in waste water treatment plants to promote optimum environment for growth of microbes within anaerobic digesters to remove pollutants. One such area that has gained prominence with regards to CHT analysis is in the development of an exhaust port system.

The aim of this challenge is to perform steady-state simulations at an input velocity value of 5m/s to understand the development of the velocity and temperature profiles for the given exhaust port model and to calculate the wall/surface heat transfer coefficient on the internal solid surface. The exhaust port is essentially a construction of a piping system in a configuration appropriate with the design requirements and constraints of the concerned application to remove or guide the burnt gases away from the combustion chamber. The given exhaust port model consists of the exhaust manifold, collector, exhaust outlet pipe and the flanges for mounting in a four-into-one setup. The analysis involves verification of the computed results with a mesh independence study detailing on the mesh refinement approach based on the 'y+' value, boundary layer thickness and required number of inflation layers to capture the boundary layer flow properties. The mesh independence study would also outline cell count, wall heat transfer coefficient values and number of iterations for convergence for each mesh based on stability of behaviour of wall heat transfer coefficient values as well as some of the mesh quality attributes like element quality, skewness, orthogonality and aspect ratio. The mesh independent result would also be validated against the analytical solution of wall heat transfer coefficient for the proposed exhaust port model configuration.

2. Method

2.1 Modelling

The exhaust port setup was modelled in SpaceClaim, the default solid modeling CAD software for the Fluid Flow (Fluent) Analysis System. Since, the exhaust port setup has different pipe lengths in the exhaust manifold and that the overall contribution from all the pipes affects the quantity and quality of discharged flow gases from the exhaust outlet pipe despite having isothermal boundary conditions and uniform input flow of constant magnitude in all directions defined in a three-dimensional space, the simulations cannot be carried out for one half of the model i.e., the model size cannot be reduced for a faster computational time. Therefore, the entire model needs to be utilised for discretization and solution of the partial differential equations which will be discussed in detail regarding the setup in subsequent sections. SpaceClaim offers a feature to create a Fluid Volume using the Volume Extract option to distinguish the volume through which the fluid flow will occur against the entire model as illustrated in Fig. 1. In this case, both the solid body and the Fluid Volume will be required for discretization as wall heat transfer coefficient needs to be evaluated at the shared interface due to the nature of the simulation i.e., CHT analysis. Therefore, it is important to ensure that the nodes on both sides of the interface in the mesh match, or in other words is conformal, which also implies negligible to minimum interpolation in this region thereby improving computational time and rendering higher accuracy to the solution. For obtaining the conformal mesh, both the solid body and the Fluid Volume in the model tree were selected and moved to a new component which was renamed to 'Final Geometry'. To allow for exchange of information during the solution process owing to the fluid-solid interactions and avoid failure to create mesh during the meshing process, the 'Share' option in both 'Share Topology' under the 'Properties' window and in 'Sharing' section under the 'Workbench' tab respectively was enabled. Lastly, the original components which were no longer in use were deleted by 'right-clicking' on the model tree and choosing the 'Delete Empty Components' option.

Fig. 1: Solid Body and Fluid Volume Models, and Boundary Conditions regions

2.2 Meshing

The meshing approach adopted for the given computational domain was developed to maintain the meshing quality attributes viz., element quality, skewness, orthogonality and aspect ratio within the ranges as laid out in the Meshing training material of ANSYS outlined in the following table:

Element Quality			0 (Bad Mesh)	1 (Good Mesh)
Aspect Ratio		Triangles	1 (for an equilateral triangle)	- Best possible value
		Quadrilaterals	1 (for a square)	- Best possible value
Skewness	0-0.25 (Excellent)	0.25-0.50 (Very Good)	0.50-0.80 (Good)	0.80-0.94 (Acceptable)	0.95-0.97 (Bad)	0.98-1.00 (Unacceptable)
Orthogonal Quality	0-0.001 (Unacceptable)	0.001-0.14 (Bad)	0.15-0.20 (Acceptable)	0.20-0.69 (Good)	0.70-0.95 (Very Good)	0.95-1.00 (Excellent)

Table 1: Mesh Quality Attributes as defined in the ANSYS meshing training material

The reason for targeting the above mentioned mesh quality attributes comes from analysis across different types of geometries/models that these are few of the known ways of identifying whether the discretized mesh would develop convergence issues during the solution process and also to locate areas in the mesh where the statistics can be improved. For each of the mesh, named selections such as inlets, outlet and outer-wall-convection have been created so as to assign their respective boundary conditions in the setup stage. The shared interface between the solid and fluid domains is detected by the Fluent solver as the region of coupled thermal condition since exchange of heat transfer takes place along this boundary further highlighting the importance of conjugate heat transfer analysis. Also, depending on the name of the named selection, Fluent is able to assign inlet, outlet, and wall boundary conditions with velocity-inlet, pressure-outlet, and convection wall type boundary conditions respectively which is a useful feature to have in flow-type problems of similar nature since the user now can only add values relating to the specific boundary condition thereby making the setup process a lot faster. In Section 2.5.2, we will discuss in detail about the different meshes that were prepared, their mesh qualities and the arrival on mesh independent solution with regards to the boundary layer development.

2.3 CFD Simulation Setup, Boundary Conditions, Monitoring Points and Solution

The simulations on all meshes was performed using the Reynolds Averaged Navier-Stokes (RANS) models, in this case, the k-`omega` turbulence models. The use of k-`omega` turbulence models and more specifically the k-`omega` SST turbulence model is attributed to the fact that the Reynolds number, based on the given flow conditions, lies in the turbulence flow region and for calculating heat transfer on the wall, the boundary layer requires to be resolved to the level of the viscous sublayer for which the appropriate `y^+` selected, in this case 1 at the first cell height, will only yield satisfactory results based on this turbulence modelling approach. Additionally, all simulations are being run using the double-precision feature to accurately measure the prediction quantities (wall heat transfer coefficient) and as a general measure to resolve large variations in pressure differences. A brief description into the working of the k-`omega` turbulence model is provided below:

2.3.1 k-`omega` turbulence model

The k-`omega` turbulence model consists of two transport equations to measure the convection and diffusion properties of turbulent energy viz., the turbulent kinetic energy (k) term, similar to the one in the k-`epsilon` turbulence models, to evaluate the energy in turbulence and `omega` to determine the turbulence scale. The `omega` equation in the standard k-`omega` turbulence model has a major downside in the form of a high sensitivity of the solution dependent on the freestream values of k and `omega` outside the shear layer. This drawback has been addressed particularly well by the k-`omega` SST turbulence model wherein it combines terms of the `omega` and `epsilon` equation. Within the freestream, the k-`omega` SST turbulence model makes use of the `epsilon` equation which is where the `omega` equation is most sensitive to while the `omega` equation is employed in the rest of the region. This switching behaviour of the k-`omega` SST turbulence model makes it particularly effective in dealing with adverse pressure gradient and separating flow related applications. Also, when compared with the k-`epsilon` turbulence model, the effects of turbulence in stagnation and strong acceleration zones are relatively smaller. The SST model in the k-`omega` SST turbulence model deals with the calculation of the principal shear stress in adverse pressure gradient boundary layers.

2.3.1.1 Setup of the k-`omega` SST turbulence model:

The k-`omega` model is available in the Model section of Viscous Model feature and when we choose the SST option in the type, we now set the k-`omega` model' to k-`omega` SST turbulence model. Production limiter in the options section is an option selected by default and essentially limits the generation of excessive turbulence energy around stagnation points.

2.3.2 Boundary Conditions and Monitoring Points

The different boundary conditions for the simulation of flow over a cylinder is described in the following table:

Named Selections	inlets	wall-fff-final_geometry_fluid-volume-fff-final_geometry_solid-volume	outer-wall-convection	wall-fff-final_geometry_solid-volume	outlet
Type	Velocity-inlet	wall	wall	wall	pressure-outlet
Velocity Specification Method	Magnitude, Normal to Boundary	--	--	--	--
Velocity Magnitude (m/s)	5	--	--	--	--
Turbulence Specification Method	Intensity and Viscosity Ratio	--	--	--	--
Turbulent Intensity (%)	5	--	--	--	--
Turbulent Viscosity Ratio	10	--	--	--	--
Backflow Turbulent Intensity (%)	--	--	--	--	5
Backflow Turbulent Viscosity Ratio	--	--	--	--	10
Temperature (K)	700	--	--	--	--
Backflow Total Temperature (K)	--	--	--	--	300
Thermal Conditions	--	Coupled	Convection	Heat Flux	--
Heat Transfer Coefficient (`W/(m^2K)`)	--	--	20	--	--
Free Stream Temperature (K)	--	--	300	--	--
Heat Flux (`W/m^2`)	--	--	--	0	--
Gauge Pressure (Pa)	--	--	--	--	0
Wall Motion	--	--	--	--	--
Shear Condition	--	--	--	--	--

Table 2: Boundary Conditions for conjugate heat transfer analysis in the given Exhaust port model

For the given problem, the material properties flowing through the fluid control volume is retained as air by default with Density value of 1.225 `(kg)/m^3`, Viscosity (or dynamic viscosity) value of `1.7894*10^(-5) (kg)/(ms)`, Specific Heat of 1006.43 `J/(kgK)` and Thermal Conductivity of 0.0242 `W/(mK)`. Similarly, for the solid body, the material selected is aluminium with Density value of 2719 `(kg)/m^3`, Specific Heat of 871 `J/(kgK)` and Thermal Conductivity of 202.4 `W/(mK)`. Both the materials air and aluminium can be verified by heading over to the "fff-final_geometry_fluid-volume" and "fff-final_geometry_solid-volume" zone of "Cell Zones" respectively under the "Zones" section of the "Physics" tab and click on "Edit" button for viewing the different attributes of the same.

The choice of Turbulence specification method chosen for this simulation is Intensity and Viscosity Ratio wherein the Intensity (or Turbulent Intensity) is defined as the ratio of root-mean-square (RMS) value of velocity fluctuations to average flow velocity while Viscosity Ratio (or Turbulent Viscosity Ratio) is defined as the ratio of the turbulent viscosity to the laminar viscosity. Turbulent intensity is set to a value of 5% (representing medium intensity) and the Turbulent Viscosity Ratio has its default value at 10 since this value was not directly available from any experimental findings and in that case, for an internal flow simulation, this value usually lies in a typically good scale or range as is recommended in the lecture series by ANSYS on Turbulence Modeling.

During the solution process, we are interested in monitoring the wall heat transfer coefficient plot along the interface between the solid body and the Fluid Volume to observe stability in output of the predicted value. To create the required monitoring plot, we head over to the 'Definitions' drop-down menu in the 'Reports' section of the 'Solution' tab and choose 'Area-Weighted Average' option in the 'Surface Report' followed by assigning the name as wall heat transfer coefficient, choosing 'Wall Fluxes' as the Field Variable and 'Surface Heat Transfer Coef.' as the attribute, selecting 'wall-fff-final_geometry_fluid-volume-fff-final_geometry_solid-volume' in the 'Wall Zones' list, checking the 'Report File' and 'Report Plot' boxes before clicking on the 'OK' button. Fig. 1 outlines the named selections where the boundary conditions would be employed as well as contour plots to demonstrate the velocity and temperature distribution discussed in Section 2.3.3 Solution and Section 2.4 Post-Processing respectively.

2.3.3 Solution

The nature of the given CFD simulation is steady-state while the type of solver chosen is Pressure-based solver as we are dealing with a low-speed incompressible flow problem i.e., the speed of the fluid flow is typically in the subsonic range and the density of the material does not change or the change is almost negligible within the given control volume, although the Fluent Theory Guide also claims that the Pressure-based solver can now solve and operate beyond its previous limitations. The Pressure-based solver calculates the Pressure field by evaluating the pressure or pressure-correction equation that is determined from the manipulation of continuity and momentum equations. Velocity formulation is chosen to be absolute as the velocity in the given fluid control volume is not rotating. Before we run the solution, we need to initialise the given fluid volume which in this case is done using the Hybrid Initialisation option. Hybrid Initialisation as described in the Fluent Theory Guide solves a Laplace equation to create a velocity flow field that will adhere to the domain geometry and a pressure field across the same domain to allow smooth transition of high and low pressure values. A contour plot to depict the development of temperature gradient across the exhaust port model during the simulation run and recorded as a MPEG animation file was setup. The contour plot was now available to be created under the 'Graphics' section of the 'Results' tab since the solution was initialized. After a new contour plot option has been selected from the 'Contours' drop-down menu, a name can be assigned, in this case, 'temperature-contour', followed by clicking on the 'Filled' check-box, choosing 'Temperature' as the field variable with 'Static Temperature' as the attribute and selecting all the surfaces in the 'Surfaces' list. We can now click on the 'Save/Display' button. To create an animation plot of the same contour, we can head back to the 'Solution' tab and click on the 'Solution Animations' option available under the 'Create' drop-down menu of the 'Activities' section. In the 'Animation Definition' window, we assign the name of the animation as 'Temperature' and choose 'temperature-contour' which is the name of the contour plot before clicking on 'OK'. Additionally, during the entire run, we will not be checking on the convergence criteria which can be disabled by selecting the 'Residuals' option in the 'Reports' section of the 'Solution' tab and unchecking the 'continuity', 'x-velocity', 'y-velocity', 'k', `omega`, and 'energy' for k-`omega` turbulence model to monitor the solution based on the number of iterations.

2.4 Post-Processing

In order to plot the velocity and temperature profiles, streamline, and velocity, temperature, and surface heat transfer coefficient contour plots across the domain length for the given fluid flow problem at different surfaces, Line and Plane locations at different positions in the exhaust port model were drawn. Apart from the Line and Plane locations described in the tables below, expressions to represent the velocity and temperature profiles in the direction of flow were created. Legends, titles, x and y labels have been added for velocity and temperature profiles, streamlines and contour plots for a better understanding of the graphical representation and can be viewed under Section 3 Results.

Expressions:

Maximum temperature and velocity (across Fluid Volume domain):

maximum temperature = maxVal(Temperature)@fff final_geometry_fluid volume

maximum velocity = maxVal(Velocity)@fff final_geometry_fluid volume

Velocity and temperature values against axial distance in the direction of flow normalised by model inlet diameter (h = 0.17m):

y by h = (Velocity /(maximum velocity / 1 [m]) ) + (Y/0.17)

y by h temp = (Temperature /(maximum temperature ) )*1[m] + (Y/0.17)

z by h = (Velocity /(maximum velocity / 1 [m]) ) - (Z/0.17)

z by h temp = (Temperature /(maximum temperature ) )*1[m] - (Z/0.17)

y-axial distance normalised by model inlet diameter (h = 0.17m):

y by h = Y/0.17

x-axial distance normalised by model inlet diameter (h = 0.17m):

x by h = X/0.17

It is important to note that the first term on the right hand side for the 'y by h' and 'z by h' expressions to evaluate normalised velocity and temperature values is actually a dimensionless number, however CFD Post does not recognise the same unless 'maximum velocity' is not divided by 1m for velocity and similarly the resulting expression for temperature multiplied by 1m respectively, since it performs the addition of similar quantities only i.e., units on both sides of the '+' operator need to be in m.

Definition of Lines at different distances along the spanwise flow direction (Refer Fig.1 indicated by yellow coloured lines):

Line 1 at the pipe inlet:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0	1.25
Point 2	0.305	0	1.25

Table 3: Definition of Line 1

Line 2 @0.1m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.1	1.245
Point 2	0.305	0.1	1.245

Table 4: Definition of Line 2

Line 3 @0.2m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.2	1.23
Point 2	0.305	0.2	1.23

Table 5: Definition of Line 3

Line 4 @0.3m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.3	1.19
Point 2	0.305	0.3	1.19

Table 6: Definition of Line 4

Line 5 @0.4m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.4	1.11
Point 2	0.305	0.4	1.11

Table 7: Definition of Line 5

Line 6 @0.465m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.465	1
Point 2	0.305	0.465	1

Table 8: Definition of Line 6

Line 7 @0.59m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.475	0.875
Point 2	0.305	0.475	0.875

Table 9: Definition of Line 7

Line 8 @0.715m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.475	0.75
Point 2	0.305	0.475	0.75

Table 10: Definition of Line 8

Line 9 @0.84m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.475	0.625
Point 2	0.305	0.475	0.625

Table 11: Definition of Line 9

Line 10 @0.965m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.475	0.5
Point 2	0.305	0.475	0.5

Table 12: Definition of Line 10

Line 11 @1.09m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.475	0.375
Point 2	0.305	0.475	0.375

Table 13: Definition of Line 11

Line 12 @1.215m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.475	0.25
Point 2	0.305	0.475	0.25

Table 14: Definition of Line 12

Line 13 @1.34m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.475	0.125
Point 2	0.305	0.475	0.125

Table 15: Definition of Line 13

Line 14 @1.465m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.135	0.475	0
Point 2	0.305	0.475	0

Table 16: Definition of Line 14

Line 15 @1.475m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.145	0.6	0
Point 2	0.3	0.6	0

Table 17: Definition of Line 15

Line 16 @1.6m from pipe Inlet downstream:

Point	x[m]	y[m]	z[m]
Point 1	0.145	0.725	0
Point 2	0.3	0.725	0

Table 18: Definition of Line 16

Definition of Plane:

Plane 1 at centre of the exhaust outlet pipe in the flow direction (Refer Fig.1 indicated by green coloured planes):

Option = XY Plane

Plane	X[m]	Y[m]	Z[m]
Plane 1	0	0.0	0.0

Table 19: Definition of Plane 1

Streamline plot is created using Type "3D Streamline" to visualise the individual flow profiles from the inlet to the collector wherein the mixing of exhaust gases takes place before being discharged from the exhaust outlet pipe. In the streamline definition, the "Variable" option is selected as "Velocity" and the option for "Cross Periodics" is checked on along with the "Forward" option under "Direction" drop-down retained.

A contour plot that displays the temperature gradients across the entire solid domain was created. For the contour plot, the "Range" was set to "Local" which was available under the "Geometry" tab. An optimal value of 101 was set for "# of Contours" and was decided primarily based on the visibility of "Show Contour Bands". Also, velocity and temperature gradient values was created using the Plane 1 definition to observe the variations in the profiles from those noticed at the inlet with the Streamline plot.

2.5 Verification

In this numerical simulation, we have carried out verification of calculation (or mesh independence study) based on the definition of boundary layer development for internal flow in the given exhaust port model and the resolution of the boundary layer to the viscous sublayer with k-`omega` SST being the turbulence model of choice more of which has been highlighted in Section 2.5.2 Mesh Refinement approach.

2.5.1 Baseline Mesh

A baseline mesh for the given computational domain was created with default element size of 0.15m and without inflation layers, illustrated in Fig. 2. Since, the fluid flow and exchange of heat transfer along the interface between the solid body and Fluid Volume has many curvature bound regions, 'Capture Curvature' option was set to 'Yes'. For the given numerical analysis, three different schemes viz., SIMPLE, SIMPLEC and COUPLED, were compared to understand the effect on convergence and stability of wall heat transfer coefficient plot as they are preferred algorithms for solving steady-state numerical problems unlike PISO which is generally used for transient calculations with large time steps. A detailed description of each of the three schemes have been mentioned in Section 2.5.3 Convergence Level.

Fig. 2: Baseline Mesh

Mesh statistics, coefficient of drag and coefficient of lift values for Baseline Mesh:

Mesh	Element-Size(m)	Nodes	Elements	wall heat transfer coef. (SIMPLE solver) [W m^-2 K^-1]	wall heat transfer coef. (SIMPLEC solver)[W m^-2 K^-1]	wall heat transfer coef. (COUPLED solver)[W m^-2 K^-1]
Baseline	0.15	27452	137651	21.2911	21.28963	21.28916

Table 20: Mesh statistics and wall heat transfer coefficient values (SIMPLE, SIMPLEC, COUPLED) for Baseline mesh of the exhaust port setup

Mesh quality attributes for Baseline Mesh:

Mesh	Element Quality(avg.)	Aspect Ratio(avg.)	Skewness(avg.)	Orthogonal Quality(avg.)
Baseline	0.81993	Triangles-1.8974	0.25067	0.74799

Table 21: Mesh quality attributes depicting the average values of element quality, aspect ratio, skewness and orthogonal quality for Baseline mesh of the exhaust port setup

Section 2.5.2 Mesh Refinement approach

The initial results obtained from the Baseline mesh in Section 2.5.1 Baseline Mesh with regards to the convergence behaviour of the three different algorithms, SIMPLE, SIMPLEC and COUPLED as can be seen in Fig. 10, Fig. 11 and Fig. 12 respectively, it was observed that the monitor plots with all three solvers show a highly stable output well before they reach the total number of iterations with SIMPLEC and COUPLED algorithms predicting almost similar wall heat transfer coefficient values from Table 20. However, the residuals for SIMPLE and SIMPLEC schemes are converging between the '1e-01' and '1e-02' values for continuity equation and '1e-04' and '1e-05' values for energy equation (although the convergence criteria option for all equations were checked off) whereas for the COUPLED scheme the trends are closer to the convergence limits based on the residual values mentioned in Table 26 making the results obtained by the latter algorithm more reliable and accurate for it to be reused to solve the pressure-velocity equations for the four different meshes that were prepared to obtain a mesh independent solution and later to be compared against the analytical solution. The time taken to solve the steady state problem was also factored into consideration before choosing the COUPLED solver which was more than that taken by SIMPLE and SIMPLEC solvers as was the expected behaviour and also mentioned in Section 2.5.3 Convergence Level. The mesh refinement hereafter was created while adhering to the mesh quality attributes mentioned in Section 2.2 Meshing. For the inlets and outlet named selections, "Edge Sizing" was chosen to uniformly divide the meshing elements and to minimize errors due to numerical diffusion. The "Edge Sizing" type was selected as "Number of Divisions" and the value was set to 36 implying that the division was taking place at a point every 10° on the circle. Inflation layers were added to capture the flow behaviour and exchange of heat transfer accurately between the solid body and Fluid Volume also for prediction of wall heat transfer coefficient values. Lastly, both the solid body and the Fluid Volume was discretized using the "Body Sizing" feature with 'Capture Curvature' option set to 'Yes' which because of the relatively complex features led to mesh creation of predominantly tetrahedral cells.

In order to add inflation layers, the boundary layer development for flow within the exhaust port setup was required to be evaluated based on the governing conditions and flow properties which is shown below:

Reynolds number (based on inlet pipe diameter):

`Re = (U*D*rho)/mu`, where (Eq. 1)

`U` is the freestream velocity,

`D` is the characteristic length, in this case, the inlet pipe diameter,

`rho` is Density of fluid, and

`mu` is the dynamic viscosity of fluid

Since, the Reynolds number for the given flow conditions within the exhaust port setup by substituting the given values in (Eq.1) equates to 58189.896, it implies that the fluid flow lies in the turbulent flow regime as was observed from experimental findings for flow in a pipe. 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 as expressed below:

`(delta/x) = 0.38/(Re_x)^(1/5)`, where (Eq. 2)

`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

The point `x` depends on the flow characteristics developing within the smooth pipe. The fluid entering the pipe has as a uniform profile before it encounters the walls of the solid body and therefore until this point, the flow can be considered "inviscid" i.e., there is no presence of viscous effects. The fluid undergoes changes when moving along the walls of the pipe as a result of the "no slip" condition and due to the presence of viscous effects, the velocity of flow at the surface comes to a value of 0. The flow adjacent to the wall begins to decelerate causing a change in the velocity profile from uniform to a flatter or fuller profile (characteristic of the turbulent flow) with every layer from the wall to the center of the pipe increasing in velocity leading to the development of a boundary layer. The length from the entrance of the pipe, wherein the inviscid core is strong, to which the boundary layers from the opposite walls of the pipe merge to form a fully developed flow, within which the viscous effects dominate, is termed 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 pipe length is greater than 10 times the pipe diameter can be expressed mathematically as:

`L_e approx 10D` (Eq. 3)

Substituting the value of 'D' in (Eq. 3), the entrance length for the given exhaust port setup was evaluated to a value of 1.7. Similarly, utilising the value of `x` in place of `D` in (Eq. 1), the local Reynolds number was calculated to be 581898.96. 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. 2 was determined to be 0.0454m.

In order to capture the boundary layer development and the exchange of heat transfer between the solid body and the Fluid Volume, we need to resolve the viscous sub layer for which the first grid cell, in this case at the first cell height, needs to be at a value of `y^+` equal to 1 based on the k-`omega` SST turbulence model selected. To ensure `y^+` value of 1, we need to calculate the values of wall shear stress and frictional velocity as shown below:

Wall Shear Stress:

`tau_w=(1/2)⋅C_f⋅rho⋅(U_∞)^2`, where (Eq. 4)

`tau_w` is Wall shear stress,

`C_f` is Skin Friction coefficient,

`rho` is Density of fluid, and

`U_∞` is the freestream velocity

Skin Friction Coefficient (`C_f`) based on a semi-empirical expression can be related to the Reynolds number as shown below:

`C_f=0.0594/(Re_x)^(1/5)` (Eq. 5)

Frictional velocity (`U_tau`) can be computed mathematically as:

`U_tau=sqrt(tau_w/rho)` (Eq. 6)

`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. 7)

`y_P` is distance from the cell-centroid 'P' to the wall such that `Delta_y = 2*y_P`. In other words, `y^(+)` at `y_P` is ensured to be 0.5, which is well within the demarcated boundaries for viscous sub-layer, such that a corresponding `y^(+)` at `Delta_y` would be equal to 1.

Simplifying (Eq. 6) to determine First Cell Height (`Delta_y`), we get;

`Delta_y = 2*(y^(+)*mu)/(rho*U_tau)` (Eq. 8)

For local Reynolds number 581898.96 with dynamic viscosity value of `1.7894*10^(-5) (kg)/(ms)`, density value of 1.225 `(kg)/m^3`, free-stream velocity values of 5`m/s` and `y^+` value of 1, the skin-friction coefficient (`C_f`), wall shear stress (`tau_w`), friction velocity (`U_tau`) and eventually the first cell height (`Delta_y`) values using (Eq. 5), (Eq. 4), (Eq. 6) and (Eq. 8) respectively are calculated in the following table:

Reynolds Number	Skin-friction Coefficient (`C_f`)	Wall Shear Stress (`tau_w`)[kg m^-1 s^-2]	Friction Velocity (`U_tau`)[m s^-1]	First Cell Height (`Delta_y`)[m]
581898.96	`4.1765*10^(-3)`	0.0639	0.228	`6.4067*10^(-5)`

Table 22: Skin-friction coefficient, wall shear stress, friction velocity and first cell height values for the governing flow conditions of the exhaust port setup

Calculation of Inflation Layers based on First Cell Height (`Delta_y`) and Turbulent Boundary Layer Thickness (`delta`):

From ANSYS Meshing User's Guide, the default growth rate type for determining height of inflation layers is Geometric. By mathematical formulation, the total height at layer 'n' for geometric growth rate type can be given by:

`H=h*((1−r^n)/(1−r))`, where (Eq. 9)

H = Total Height at layer 'n' i.e., Turbulent boundary layer thickness (`delta`),

r = height ratio = 1.2,

h = Initial height = First Cell Height i.e., `Delta_y`, and

n = Layer number

The missing variable in (Eq. 9) is n, or in this case total number of inflation layers. Equating for n in (Eq. 9), we get total number of inflation layers as a round figure value of 28. Due to constraints in computational resources and limitations in maximum number of cells/nodes for solution by ANSYS Student license to 512,000, it was difficult to achieve 28 inflation layers with an appropriate transition ratio in body element size value. Therefore, the mesh refinement approach involved increase in inflation layers from a starting value of 10, considered good practice from different applications perspective, while maintaining the turbulent boundary layer thickness along with gradual decline in body element size value and accordingly four different meshes were prepared to identify the mesh independent solution, the details of which are reported in Table 23.

Mesh	No. of Inflation Layers	Turbulent Boundary Layer Thickness [m]	Body Sizing Element Size [m]	% increase in element size (w.r.t Turbulent Boundary Layer Thickness)
1	10	0.0454	0.08172	80
2	12	0.0454	0.07718	70
3	14	0.0454	0.07264	60
4	16	0.0454	0.0681	50

Table 23: Body Sizing Element Size and Mesh Element Size for different meshes of the exhaust port setup with regards to Turbulent Boundary Layer Thickness

Mesh 1:

Fig. 6: Mesh 1

Mesh 2:

Fig. 7: Mesh 2

Mesh 3:

Fig. 8: Mesh 3

Mesh 4:

Fig. 9: Mesh 4

Mesh statistics and wall heat transfer coefficient values for the four different meshes:

Mesh	Nodes	Elements	wall heat transfer coefficient [W m^-2 K^-1]	% var
1	145781	331903	21.292	--
2	171396	382411	21.305	+0.00061
3	196539	431714	21.101	-0.00957
4	221353	480811	21.243	-0.00673

Table 24: Mesh statistics and wall heat transfer coefficient values for different meshes of the exhaust port setup

Mesh quality attributes for the four different meshes:

Mesh	Element Quality(avg.)	Aspect Ratio(avg.)	Skewness(avg.)	Orthogonal Quality(avg.)
1	0.3482	Triangles and Quadrilaterals - 198.83	0.24755	0.56544
2	0.30639	Triangles and Quadrilaterals - 244.78	0.23985	0.55761
3	0.27264	Triangles and Quadrilaterals - 299.2	0.23284	0.55004
4	0.24573	Triangles and Quadrilaterals - 370.06	0.22741	0.54084

Table 25: Mesh quality attributes depicting the average values of element quality, aspect ratio, skewness and orthogonal quality for the different meshes of the exhaust port setup

From the meshing statistics in Table 24, despite mesh number 1 (or mesh 1) with 10 inflation layers and Body Sizing element size of 0.08172m representing minimum percentage variation (%var) of the wall heat transfer coefficient values, it was decided that mesh number 3 (or mesh 3) with 14 inflation layers and Body Sizing element size of 0.07264m represents a mesh independent result as the percentage variation in the predicted value from mesh number 2 (or mesh 2) and also from mesh 3 to mesh 4 lies in the same range or order with more number of cells to compute the required solution. 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., ((wall heat transfer coefficient for Mesh 2 - wall heat transfer coefficient for Mesh 1)/wall heat transfer coefficient for Mesh 1). Interestingly, the quality attributes of all four meshes in Table 25, except for Skewness, represent a general trend of decline especially when compared against Baseline mesh in Table 21 as well as no significant difference in prediction of the wall heat transfer coefficient value. This can be attributed to the fact that the 'Capture Curvature' option, which is a type of 'Size Function', determines the Edge and Face sizes based on 'Curvature Normal Angle' and as such subdivides the mesh in curved regions until each individual element will span with the maximum value of this angle, in this case the default value of 18°. In other words, regions of less curvature prominence will result in larger elements while high curvature and within small gaps the elements will be smaller. Inflation layers, on the other hand, make use of high aspect ratio cells, more specifically tetrahedral cells and are also referred to as a layer of prism elements due to the stacking behaviour, in order to capture boundary layer flow properties and wall heat transfer coefficient. Since, inflation layers have been defined with the 'Total Thickness' option which offers a consistent layer of prism elements while 'Capture Curvature' option varies the size of the meshing elements depending on the model design, the transition between the inflation layers to the layer outside it is inconsistent causing large interpolation of values especially in regions moving from high curvature to less curvature. It is also important to note that despite the "Body Sizing" feature being implemented and element size was reduced in subsequent mesh refinement within the 512,000 cell count limit set by ANSYS Student license, the large interpolation of values still persists as an attempt to adopt "best fit" of meshing elements increases the prevailing inconsistencies. Results in Section 3 will only be referred with regards to the mesh independent solution i.e., Mesh 3.

2.5.3 Convergence Level

As mentioned earlier in Section 2.5.1 Baseline Mesh, three different schemes i.e., SIMPLE, SIMPLEC and COUPLED were compared with regards to convergence study. SIMPLE scheme is iterative in nature, robust (memory efficient) and is quicker per iteration than other pressure-based solvers. SIMPLEC scheme, similar to the SIMPLE scheme, allows increased under relaxation to the pressure correction term that improves in convergence times especially when employed in complicated flows with turbulence and/or supplementary physical models and as such offers greater diversity in terms of computational problems for implementation over SIMPLE scheme. While SIMPLE and SIMPLEC schemes are segregated type of schemes i.e., pressure and velocity terms are updated sequentially, the COUPLED scheme on the other hand updates the pressure and velocity terms simultaneously. This allows the COUPLED scheme to be more efficient in its single phase implementation especially for steady-state flows offering superior performance over the segregated algorithms (SIMPLE and SIMPLEC) thereby requiring fewer iterations to converge. The selection of schemes for the different terms under spatial discretization includes assigning Gradient to Least Squares Cell Based scheme, pressure to second order scheme, momentum to second order upwind scheme, Turbulent Kinetic Energy to first order upwind scheme and Turbulent Dissipation Rate to first order upwind scheme of which the Turbulent Kinetic Energy and Turbulent Dissipation Rate. Least Squares Cell Based scheme is the default gradient method chosen by ANSYS Fluent which is used to compute the change in values across the centroid of two cells along a vector path and is a computationally less expensive scheme offering comparable accuracy to the node-based gradient scheme, both of which are superior to the cell-based gradient scheme. The second order scheme for the pressure makes use of the central differencing scheme and is generally considered to be more accurate than the Standard and Linear schemes as well as robust in the case of skewed meshes. A second order upwind scheme is chosen for the momentum and energy terms that makes use of the quantities at the cell face values to compute a cell-centered solution at the cell centroid leading to a second order accurate solution. The turbulence terms viz., Turbulent Kinetic Energy and Turbulent Dissipation Rate are set at first order upwind scheme by default which makes use of quantities at the cell faces as the cell-average value indicating that it is first-order accurate. The reason to have the turbulence terms discretized using the first-order upwind scheme is based on observations from simulations across different applications when using a higher or even second order upwind scheme leads the computational cost to become more expensive (more compute time) solution while also demonstrating stability issues in the form of repeating oscillations of values and not necessarily guaranteeing always a more accurate solution. Convergence criteria option was checked off, however, the residual values which measures the local imbalance of conserved variable in each control volume such as continuity, x-velocity, y-velocity, k, `omega` and energy were monitored throughout the solution process for the defined number of iterations to allow the actual fluid behaviour to be developed i.e., the formation of a fully developed flow profile and consistent rate of heat transfer along the interface between the solid body and the Fluid Volume. The residual limits defined for the turbulence model equations and the residual plots corresponding to the different meshes are shown in the table and figures below:

	k-`omega` SST
continuity	1e-03
x-velocity	1e-03
y-velocity	1e-03
k	1e-03
omega	1e-03
energy	1e-06

Table 26: Residual Values of the different equations used in the k-`omega` SST turbulence model

Fig. 10: Scaled Residuals and wall heat transfer coefficient plots using the k-`omega` SST turbulence model for baseline mesh with SIMPLE solver

Fig. 11: Scaled Residuals and wall heat transfer coefficient plots using the k-`omega` SST turbulence model for baseline mesh with SIMPLEC solver

Fig. 12: Scaled Residuals and wall heat transfer coefficient plots using the k-`omega` SST turbulence model for baseline mesh with COUPLED solver

Fig. 13: Scaled Residuals and wall heat transfer coefficient plots using the k-`omega` SST turbulence model for mesh 1 with COUPLED solver

Fig. 14: Scaled Residuals and wall heat transfer coefficient plots using the k-`omega` SST turbulence model for mesh 2 with COUPLED solver

Fig. 15: Scaled Residuals and wall heat transfer coefficient plots using the k-`omega` SST turbulence model for mesh 3 with COUPLED solver

Fig. 16: Scaled Residuals and wall heat transfer coefficient plots using the k-`omega` SST turbulence model for mesh 4 with COUPLED solver

From Fig. 13, Fig. 14, Fig. 15 and Fig. 16, we can notice that the continuity, k and `omega` residuals do not meet with the '1e-03' default limit (although not being considered as the convergence criteria option was checked off) and the wall heat transfer coefficient plots also no longer show a highly stable output as was observed with the Baseline mesh for COUPLED solver in Fig. 12 for the same number of iterations. This relates back to the reason specified in Section 2.5.2 Mesh Refinement approach, wherein large interpolation persists between the inflation layers and the remaining mesh between the high curvature and less curvature regions leading to erratic variation in calculation of values.

3.Results

In Section 3.1, numerical simulation results for wall heat transfer coefficient are compared against analytical solution while in Section 3.2, velocity, temperature and surface heat transfer coefficient contours, and streamline plots to demonstrate the development of velocity and temperature profiles with regards to the mesh independent solution are discussed.

3.1. Validation of numerical results against analytical solution for wall heat transfer coefficient

There are three different situations through which exchange of heat takes place in the given exhaust port model. The first scenario begins with heat being conveyed from the exhaust gas to the pipe walls through forced convection during fluid flow from the inlets of the manifold to the exhaust port outlet. The convected heat then travels from the internal surface of the manifold to the external surface by conduction through the solid pipe wall which is then convected to the ambient, in this case, by natural convection as it is assumed that both the exhaust port model and the air surrounding it are in stationary conditions. In this numerical study, the variable of interest is the wall heat transfer coefficient on the internal solid surface and therefore, based on the boundary conditions available with regards to the given exhaust port model (as seen in Table 2), the heat transfer properties in the first scenario will be evaluated as follows:

Scenario 1: Forced Convection - Exchange of heat from exhaust gases to the pipe internal surface

Analytical Solution:

As the Reynolds number calculated based on the governing flow conditions and pipe inlet diameter (being the characteristic length) suggests the flow behaviour is turbulent in nature, the Dittus-Bölter heat transfer coefficient correlation can be used to calculate the Nusselt number as a property of the heat transfer attributes which is mathematically expressed as:

`Nu_D=0.023*Re_D^(4/5)*Pr^n`, where (Eq. 10)

`D` refers to the inlet pipe diameter belonging to the Fluid Volume domain,

`Pr` is the Prandtl number = `(mu*C_p)/k` with `mu` referring to the fluid viscosity, `C_p` being the specific heat capacity of the fluid at constant pressure and `k` implying thermal conductivity of the fluid = 0.744, and

`n` = constant whose value indicates whether the fluid is being heated or cooled = 0.3 (in this case as the walls are at a cooler temperature than the fluid)

The above equation is also valid for the given boundary conditions based on the limitations defined that are as follows:

`0.6<= Pr (=0.744)<=160`

`Re_D (=58189.896) > 10,000`

`(L/D)>10`

Substituting the values of Reynolds number with respect to the pipe inlet diameter, Prandtl number and constant value of `n` as 0.3, we get the Nusselt number value as 136.483. Nusselt number, as a property of heat transfer, can be defined as the non-dimensional ratio of convective heat transfer to conductive heat transfer across the same boundary. Mathematically, Nusselt number is evaluated as:

`Nu=(h*L)/k`, where (Eq. 11)

`h` is the convective heat transfer coefficient for the given flow,

`L` refers to the characteristic length, in this case, pipe inlet diameter `D`, and

`k` is the thermal conductivity of the fluid

Simplifying (Eq. 11) to determine the convective heat transfer coefficient (`h`), we get

`h=(Nu_D*k)/D` (Eq. 12)

Solving for the convective heat transfer coefficient (`h`) in (Eq. 12) with values for the Nusselt number based on pipe inlet diameter as the characteristic length and thermal conductivity of air 0.0242 `W/(mK)`, the value obtained is 19.4287 `W/(m^2K)`. It should be noted that the Dittus-Bölter heat transfer coefficient correlation was observed and found to be accurate in the range of `+- 15%` and works best when the pipe surface is considered smooth as was assumed in this case.

Wall heat transfer coefficient (Numerical Analysis vs Analytical Solution):

Wall heat transfer coefficient [W m^-2 K^-1] - Numerical Analysis	Wall heat transfer coefficient [W m^-2 K^-1] - Analytical Solution	Absolute Difference (%)
21.418	19.4287	10.239

Table 27: Comparison of Wall heat transfer coefficient values between Numerical Analysis and Analytical for mesh independent solution of the exhaust port model

As seen in Table 27, the numerically computed value overpredicts the wall heat transfer coefficient by 10.239% when compared against the analytical solution for the same. It should be noted however that the numerically computed value represents average value while the analytical solution for the wall heat transfer coefficient generally holds good in straight sections of the pipe and more so with fully developed flow profile throughout. There are other heat transfer coefficient correlations that may be used, but given that the temperature on the inner surface of the pipe walls was unknown, the Dittus-Bölter heat transfer coefficient correlation and the prediction from numerical analysis gives a reasonably good understanding of the underlying heat exchange mechanisms especially more so from a quantitative perspective.

3.2 Temperature and velocity profiles, velocity, temperature and surface heat transfer coefficient contour plots, and streamline profile

Fig. 17, Fig. 18, Fig. 19, Fig. 20, Fig. 21 and Fig. 22 shows the velocity and temperature gradients across the exhaust port model using velocity and temperature contour plots at Plane 1, temperature contour plot across the computational domain, streamlines plot, surface heat transfer coefficient contour plot and velocity and temperature profiles at different line locations respectively. The streamlines plot, and temperature and velocity profiles along the y-axis flow indicate smooth profiles at the inlet, however, as the flow progresses downstream, the velocity and temperature profiles begin to develop boundary layer thickness which in the case of the velocity profile, because of the high Reynolds number, is turbulent in nature. In the collector zone, where the exhaust gases from all inlets meet, we begin to see a region of intense mixing that alters the velocity and temperature profiles as was observed for flow progressing in the z-axis direction. Due to the difference in the diametrical areas between the collector zone and the exhaust outlet pipe coupled with the sharp bend or curvature in the latter, a pressure drop is created thereby causing the velocity to accelerate in the exhaust pipe region as also seen in the velocity contour plot. This is also identified by the flow separating just around the corner due to the curvature leading to a lowering of pressure levels.

For Mesh 3:

Fig. 17: Velocity and Temperature contour plots for Mesh 3 at Plane 1

Fig. 18: Streamlines Profile within the Fluid Volume domain

Fig. 19: Temperature contour plot for Mesh 3 across the computational domain

Fig. 20: Surface Heat Transfer Coefficient contour plot for Mesh 3 across the computational domain at the internal solid surface

Fig. 21: Velocity and Temperature Profiles (normalised using model inlet diameter) along y-axis direction for Mesh 3

Fig. 22: Velocity and Temperature Profiles (normalised using model inlet diameter) along z-axis direction for Mesh 3

4. Discussion

As mentioned previously in Section 2.5.2 Mesh Refinement approach, the fluid enters the pipe with a uniform profile before undergoing change when it moves past the walls of the solid body and that the "no slip" condition combined with the viscous effects causes the velocity of flow at the surface to arrive at the value of 0. As the given flow is turbulent in nature, the velocity profile undergoes transition from uniform to a flatter or fuller profile (also seen in Fig. 20 and Fig. 21). Since, the flow is not fully developed at the pipe inlets, and the thickness of the boundary layers keeps increasing in the direction of flow due to increasing viscous effects, a point is reached wherein the boundary layers from the opposite walls of the pipe merge to form a fully developed flow, which is called the entrance length. The presence of a turbulent boundary layer is also characterised by high heat transfer rates as the presence of swirling regions, also called as eddies, lead to randomised and rapid transfer of mass, momentum and energy due to the intense mixing between the layers of the flow. Moreover, in the case of the given exhaust port model, the turbulent boundary layer thickness decreases in length as it approaches the bending regions and this causes increased heat transfer between the exhaust gases and the walls of the pipe. The process of the turbulent boundary layer beginning to thicken and then decrease in length is observed again in the exhaust pipe outlet region wherein the added effects of the pressure drop increases the velocity of the outgoing gases leading to subsequent increase in the heat transfer rate across the solid domain which is also indicated by rise in temperature values at the pipe walls. Since, the pipe walls are convecting heat to the stationary air surrounding the exhaust port model which is also dependent on the convective heat transfer coefficient, the temperatures along the pipe walls are controlled and are mostly higher in value towards the exhaust outlet pipe as the entire heat from the model is conveyed out from a single pipe.

5. Conclusion

The numerical analysis for CHT phenomena with respect to the given exhaust port model under steady state conditions using the k-`omega` SST turbulence model in ANSYS Fluent 2020 R1 was well understood. It was realised that based on the boundary layer development along the pipe walls over the exhaust port model and use of an appropriate viscous model, a CFD model was able to compute the velocity and temperature gradients as well as highlight changes in flow properties. Additionally, the wall heat transfer coefficient computed with the CFD model demonstrated the prediction capabilities of this mathematical model when validated against the analytical solution. However, the accuracy of the k-`omega` SST turbulence model itself depends on the accuracy of the analytical solution which, in this case, works well based on a set of limiting conditions within a particular range for both dimensional and flow properties, flowing conditions in straight and/or other geometrical patterns, friction coefficient of pipes, physics bound constraints, etc. Although, the baseline mesh was refined in order to capture the boundary layer thickness, the wall heat transfer coefficient predicted with the refined mesh was in the same range as that of the baseline mesh and the reason for the same can be ascribed to the fact that the boundary layer thickness was invariably small throughout as it kept growing and declining in length due to the varying geometrical changes to make any substantial influence on the wall heat transfer coefficient except closer to the exhaust outlet pipe region.