Week 1: Channel flow simulation using CONVERGE CFD

1. Introduction

Fluid flow within a well-bounded solid structure whose volume is completely occupied such that flow across the cross-sectional areas takes place due to pressure gradient can be characterised as internal flow. Hot and cold water flow pumped to taps in homes, transport of oil and gas distributed over long networks of piping connections can be regarded as commercial and industrial examples of internal flow.

The aim of this challenge was to perform two-dimensional, steady-state simulations, using the default physical model containing no inclusion of turbulent terms for a channel flow geometry. The analysis involved computing results for different grid sizes that would outline total cell count, pressure drop values across the domain and velocity and mass flow rate values at the inlet and outlet boundaries, number of cycles for convergence of each mesh based on the convergence tolerances as well as some of the 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:

1. Creating a new folder dedicated to running the concerned project inside a user-preferred directory.
2. 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 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.
3. Saving the case file(s) and exporting the input files into the folder created for running the concerned project which typically includes 'boundary.in', 'material_state.dat' where material_state points to either 'gas', 'liquid', 'parcel' or 'solid', 'initialize.in', 'inputs.in', 'monitor_steady_state.in', 'post.in', 'solver.in' and 'surface.dat' 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.
4. Running the simulation in a Cygwin64 terminal window.
5. 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 simulating air flow movement through a channel geometry was modelled in Converge Studio v3.0. Since, the pressure gradient across the computational domain was maintained at their initial or boundary condition values (in this case, both being the same), the simulations were carried out in a two-dimensional setup by selecting the appropriate boundary types on the different boundary zones. The extents of the channel geometry were selected as 0.1m in the direction of flow (or x-axis) and 0.01m along the remaining directional axes i.e., y and z axes while ensuring that the center of the model was positioned at the origin i.e., (0.0m, 0.0m, 0.0m) using the 'Box' shape feature when the 'Create' button was clicked on under the 'Geometry' dock. 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 were needed to be grouped together. In order to perform the grouping and creation of new boundaries, the 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 5, in this case, followed by lastly grouping the triangles that constitute a given surface and naming them accordingly, such as 'Inlet', 'Outlet', 'Top Wall' and 'Bottom Wall' (together named as TopBottomWalls), 'Front2D' and 'Back2D' as illustrated in Fig. 1. Moreover, the normal from each of the triangular surface is required to point towards the fluid domain which can be 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 can be taken by clicking on the 'Normal' tab within the 'Transform' button under the 'Geometry' dock and selecting a 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. It should be noted that the length of the channel flow model was estimated to be long enough to observe the variations in flow behaviour along the walls of the domain for different grid sizes but not necessarily leading to fully developed flow conditions as it was not the primary requirement.

Fig. 1: Fluid Volume Model and Boundary Conditions regions in two-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 [1] 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.5 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 default physical model consisting of no turbulence equations. A brief description into the working of the default physical model is provided below:

2.3.1 Default physical model

The default physical model can be described as the solution of the governing fluid equation viz., Navier-Stokes equations for a domain, in this case, discretized using the finite volume approach, without the presence of turbulence terms. The Navier-Stokes equations define the fluid flow properties such as velocity, pressure, density, temperature, viscosity, energy, etc., within a given physical domain. Additionally, the Navier-Stokes equations are developed with regards to the laws of conservation that states the net mass flow across a system is zero (Conservation of mass) and that the momentum within the same system (or control volume) is constant (Conservation of momentum). Since, energy is also a variable of assessment for the given numerical analysis, then according to the first law of thermodynamics, total heat gain of a system is equal to the sum of heat added to the system and work done on the same system.

2.3.1.1 Setup of the default physical model:

To use the default physical model in Converge Studio v3.0, we navigate to Physical Models option in the Case Setup tree and uncheck 'Turbulence modeling' to prevent inclusion of turbulence terms in the solution of governing fluid equations.

2.3.2 Boundary Conditions and Monitoring Variables

The different boundary conditions for the channel flow simulation is described in the following table:

Names	Inlet	Front2D	Back2D	TopBottomWalls	Outlet
Type	INFLOW	TWO_D	TWO_D	WALL	OUTFLOW
Pressure Boundary Condition	Specified Value (DI), Total pressure	--	--	--	Specified Value (DI)
Pressure Magnitude (Pa)	100001.0	--	--	--	100000.0
Velocity Boundary Condition	Zero normal gradient (NE)	--	--	--	Zero normal gradient (NE)
Temperature Boundary Condition	Specified Value (DI)	--	--	Specified Value (DI)	--
Temperature Backflow	--	--	--	--	Specified Value (DI)
Temperature Magnitude (K)	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	--
Surface Movement	--	--	--	Fixed, No-slip	--

Table 2: Boundary Conditions for channel flow simulation 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 and critical pressure 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. Since, the given flow analysis is for laminar conditions, the values of 'Turbulent Prandtl number' and 'Turbulent Schmidt number' are left unchanged at their default values of 0.9 and 0.78 respectively as they would play no role in energy transport equation for the former and species transport equation for the latter.

2.3.3 Determination of maximum possible flow velocity across the channel flow domain and hydrodynamic entrance length

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 major loss across the channel domain which is given by the following relation:

`H_L(MAJOR) = f*(L/D)*(V_1^2/(2*g))`; where (Eq. 4)

`L` is the length of the pipe determined as 0.1m,

`D`, in this case `D_h` as the cross-section is noncircular, is the hydraulic diameter evaluated as 0.01m i.e., (`D_h = (4*A_c)/(p) = (4*a*b)/(2*(a+b))`), and

`f` is the friction factor, in this case assuming laminar flow conditions for a non-circular channel medium with a=b as `f = 56.92/(Re_(D_h))` or `f = 56.92*mu/(V_1*D_h*rho)` as velocity `V_1` is unknown.

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 ignoring the potential energy or `rho*g*h` terms as fluid flow is horizontally aligned, we get

`P_(TOTAL) - 0.5*rho*V_1^2 + 0.5*rho*V_1^2 = P_(bound) + 0.5*rho*V_2^2 + (56.92*mu)/(V_1*D_h*rho)*(L/D_h)*(V_1^2/(2*g))` (Eq. 5)

Velocity at the inlet boundary (`V_1`) was related to the velocity at the outlet boundary (`V_2`) from the continuity equation which for a horizontal and incompressible fluid flow states that the flow rate is conserved across the domain and mathematically can be expressed as:

`A_1*V_1 = A_2*V_2`; where (Eq. 6)

`A_1` and `A_2` are the cross-sectional areas at inlet and outlet boundaries respectively which in this case are equal

Therefore, `V_1` and `V_2` can be assumed to be equal in the given scenario such that simplifying (Eq. 5) using known values of `P_(TOTAL)` as 100001 Pa, `P_(bound)` as 100000 Pa, `g` as 9.81 `m/s^2`, density (`rho`) and dynamic viscosity (`mu`) of air at 300K (or 26.85°C) as 1.1766 `(kg)/(m^3)` and `1.857*10^-5 (kg)/(m.s)` respectively [2], `V_1` at inlet boundary was determined as 1.2654m/s which was also the same as `V_2` at the outlet boundary.

The reason for determining the maximum possible flow velocity across the given channel flow domain was to estimate the Mach number (or local Mach number) so as to decide on 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 channel flow model, to speed of sound in the medium or the speed of sound in air. As the maximum possible flow velocity for the given channel flow model was evaluated as 1.2654m/s and speed of sound in air is approximately equal to 347.336 `m/s` for the same temperature conditions but 10% relative humidity [3], the Mach number resulted in an approximate value of 0.003643 which is much lesser than the typical range of 0.2 to 0.3 for which compressibility effects could be considered.

Reynolds number (based on non-circular cross-section inlet diameter):

`Re_(D_h) = (U*D_h*rho)/mu` (Eq. 7)

Substituting the calculated value of inlet velocity `U` as 1.2654m/s into (Eq. 7) with the values of other variables in the same equation mentioned previously, we get value of Reynolds number based on the hydraulic diameter as 801.7607.

Entrance region determined by the hydrodynamic entrance length is typically considered as the distance from the entrance of the channel to the position where the wall shear stress assumes about two percent of its fully developed value and is also the point wherein the value of velocity reaches its maximum value at the center of the pipe, in this case at the axis of symmetry, due to the development of a boundary layer. The entrance length which in the case of laminar flow, in non-dimensional form, can be expressed mathematically as shown below:

`(L_h)/D_h cong 0.05*Re` (Eq. 8)

Substituting the determined values of hydraulic diameter and Reynolds number as determined in (Eq. 4) and (Eq. 7) into (Eq. 8), the hydrodynamic entrance length was evaluated to a value of 0.4m from the inlet of the channel medium in the direction of flow. This implied that for the given length of the domain, the position of the hydrodynamic entrance length was situated beyond the extents defined for the same model and as such performing validation for fully developed flow conditions would not be possible in this case.

2.3.4 Solution

The nature of the given CFD simulations are steady-state which was set to in the 'Solver' drop-down menu under the 'Solver' tab of the 'Run parameters' window as a sub-option of the 'Simulation Parameters' in the 'Case Setup' tree while the 'Simulation mode' and 'Gas flow solver' in the same window are set to 'Full hydrodynamic', by default, and 'Incompressible' respectively. The option 'Full Hydrodynamic' refers to a Normal Converge simulation mode [1] as opposed to its other option of 'No Hydrodynamic solver' which is more applicable to problems involving complex surface features and consequently their mesh creations and/or problems involving moving parts whose positions greatly vary with time such that errors can be dealt with at earlier stages without having to perform numerical simulations of transport equations, spray, combustion, etc., making the process a lot more efficient while 'Incompressible' indicates the gas does not exhibit compressibility effects or essentially is incompressible due to Mach number being less than the range of 0.2-0.3 for compressibility effects to be considered as explained in Section 2.3.3 Determination of maximum possible flow velocity across the channel flow domain. As a result of choosing the 'Incompressible' option for the 'Gas flow solver', the 'Navier-Stokes solver type' under the 'Navier-Stokes solver' tab was set to 'Pressure-based' in the 'Solver parameters' section of the 'Simulation Parameters' in the 'Case Setup' tree while the 'Navier-Stokes solver scheme' in the same tab was left at its default setting of 'PISO'. The reason for adopting the 'PISO' setting in the 'Navier-Stokes solver scheme' is because it is computationally cheaper over the 'SIMPLE' setting while both of them take into consideration the hyperbolic nature of the transport equations and continue to use the elliptic nature of the pressure equation so as to disseminate information across the domain at a rapid rate, however, the latter setting is more appropriate for flows that are dominated locally by viscosity, well-resolved boundary layers and is more stable over the former in transient simulations wherein a larger time step would be more feasible or in the case of steady-state simulations, a larger pseudo time-step. The remaining parameters in the subsequent tabs of the 'Solver parameters' window are unchanged with respect to their default values except for the 'Preconditioner' drop-down menu in the 'Main' tab of selected 'Pressure' equation under the 'Equations' tab which was changed to 'None' from 'SOR' to avoid increase in computational cost and more importantly was not essentially required for the given numerical analysis. Before running the solution, the given fluid volume was required to be initialised and more importantly a region describing the entire domain was to be created on which initial values could be assigned for which the 'Regions and initialization' sub-option was clicked on in the 'Initial Conditions and Events' under the 'Case Setup' tree, followed by clicking on the 'Add' button to add a new region in which 'Velocity' was assigned a component-wise value of 0m/s along all three directional axis i.e., velocity was initialised as (0.0, 0.0, 0.0) m/s, 'Temperature' was allocated as 300K, 'Pressure' as 100001.0 Pa which was the same as the pressure value defined at the Inlet surface using the 'INFLOW' boundary type, 'Air' as species was clicked on wherein the mass fraction values for `O_2` and `N_2` were reversed or in other words mass fraction of `O_2` was 0.77 while mass fraction of `N_2` was 0.23. It should be noted that the mass fraction values defined for `O_2` and `N_2` in the 'Regions and initialization' window are only initial values which would eventually be replaced by the corresponding mass fraction values of the same species as defined at the inlet and outlet boundary names mentioned in Table 1. In order for the boundary conditions to be employed on the given fluid domain, they need to point to the same, in this case, the region defined in the 'Regions and initialization' sub-option of the 'Initial Conditions and Events' under the 'Case Setup' tree which implied that it was required to revisit the 'Boundary' sub-option of the 'Boundary Conditions' and assign each boundary name to the region created. During the solution process, it was decided to monitor the average mass flow rate and average velocity values to observe stability and accordingly assess the convergence behaviour of the same. The monitor variables, in this case the steady-state monitor variables, were created by clicking on the 'Steady-state monitor' sub-option of the 'Run parameters' setting under the 'Case Setup' tree followed by creating two monitor variables, individually, with one of them assigned the name 'MASS_FLOW_RATE' and the other 'AVG_VELOCITY', both available under the 'Monitor Variables' column as a drop-down option as they are predefined followed by assigning both of them as 'Target Type' to be 'Boundary' and 'Target' as 'Outlet' followed by clicking on the 'Validate' and 'OK' buttons successively. In the same window, the 'Minimum number of cycles to be executed for the stead-state solver' was altered to a value of 5000, and Monitor start delay (cyc/s) and Sample size (cyc/s) for both monitor variables were set to 0 and 1000 respectively indicating for the former that the monitor variables would be checked for steady-state convergence immediately after the entered value while for the latter the steady-state condition would be examined based on the number or range of samples. All other parameters and their corresponding values remained unaltered. Additionally, during the entire run, the convergence tolerances for different variables of the Navier-Stokes equations were also being reported.

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 can be 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, 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 for a two-dimensional flow problem 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' box for disabling the viewing of the plane borders. Plotting the velocity profiles at the outlet of the channel flow model 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 4, and then clicking on the "Apply" button to view the said profile. Similarly, the pressure profile across the entire domain was represented using the "plot over line" tool, however, with coordinates as depicted in Table 5.

It was ensured that the last time-step size of the simulation run was selected to view the pressure and velocity 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. Legends, titles, and x and y labels were added for better understanding of the graphical representations.

Line	x	y	z
Point 1	0.1	-0.005	0
Point 2	0.1	0.005	0

Table 3: Line Coordinates for viewing the velocity profiles

Line	x	y	z
Point 1	0	0	0
Point 2	0.1	0	0

Table 4: Line Coordinates for viewing the pressure profile

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 'OK'.

2.5 Mesh Refinement approach

In this numerical simulation, mesh refinement was carried out by utilising values for different mesh sizes, as described in the following table (Table 5) as 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 values at the outlet and pressure drop across the entire domain.

Mesh	dx(m)	dy(m)	dz(m)
Baseline	0.001	0.001	0.004
Mesh 1	0.0002	0.0002	0.0002
Mesh 2	0.00015	0.00015	0.00015
Mesh 3	0.0001	0.0001	0.0001

Table 5: Mesh sizes representing spatial grid discretization values across each of the three-directional axes

Baseline Mesh:

Fig. 2: Baseline Mesh

Mesh 1:

Fig. 3: Mesh 1

Mesh 2:

Fig. 4: Mesh 2

Mesh 3:

Fig. 5: Mesh 3

Mesh and time statistics, peak velocity magnitude values at the outlet and pressure drop across the domain for different meshes:

Mesh	Cells	Execution Time [s]	Peak Velocity Mag. [m s^-1] at outlet	Pressure drop [Pa] across domain
Baseline	1000	26.1994	1.25257	0.66
1	25000	344.4723	1.25743	0.6
2	45420	2031.0343	1.25385	0.59
3	100000	2104.4720	1.25748	0.6

Table 6: Mesh and time statistics, peak velocity magnitude values at the outlet and pressure drop across the domain for different meshes, 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
Baseline	1.00024	1.17578e-38	1.17578e-38
1	1.00024	1.17578e-38	1.17578e-38
2	4	1.17578e-38	36.8699
3	1.00024	1.17578e-38	1.17578e-38

Table 7: Mesh quality attributes depicting the maximum values of aspect ratio, skewness and non-orthogonality for the different meshes of the channel flow model

From the meshing statistics in Table 6, it can be seen that variations in the reporting of peak velocity magnitude values at the outlet and pressure drop across the domain for different meshes are similar to the second decimal place for velocity magnitude and less than one hundredth apart in the case of pressure drop respectively implying that a more refined mesh does not have a significant impact on the expected results, however, becomes computationally more expensive as a consequence of increase in execution time. Additionally, mesh quality attributes of all the different meshes prepared in Table 7, except Mesh 2, indicates that the maximum aspect ratio, maximum skewness and maximum non-orthogonality remain unaltered due to the quality of cells across the domain continuing to remain of the same type i.e., quadrilateral in two-dimensional form, while they continue to get refined to smaller spatial discretization values. It should be noted that a reason for deviation in mesh quality attributes for Mesh 2 from the remaining meshes could not be ascertained despite the setup being consistently same across all meshes and that it had no impact in terms of prediction of values of the previously mentioned quantities.

2.7 Convergence Level

In the given numerical analysis, a steady-steady simulation was performed using a Pressure-based solver 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 incompressible flow requires the pressure 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 and Energy 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. 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.3 Boundary Conditions and Monitoring Points, the behaviour of velocity, density and pressure 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 default physical 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

Table 8: Solver Types, Convergence Tolerance Values, Minimum and Maximum Iterations and Preconditioners of the different equations used in the default physical model

Fig. 6: Plots of monitor variables velocity, density, mass flow rate and pressure for steady-state simulation using the default physical model for baseline mesh

Fig. 7: Plots of monitor variables velocity, density, mass flow rate and pressure for steady-state simulation using the default physical model for mesh 1

Fig. 8: Plots of monitor variables velocity, density, mass flow rate and pressure for steady-state simulation using the default physical model for mesh 2

Fig. 9: Plots of monitor variables velocity, density, mass flow rate and pressure for steady-state simulation using the default physical model for mesh 3

From Fig. 6, Fig. 7, Fig. 8 and Fig. 9, we can notice that the monitor variables show a highly stable output 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 has been obtained. It should be noted that for Mesh 3, the number of Maximum Iterations to solve Pressure equations using CONVERGE BiCGSTAB solver was increased from 300 to 500 in order to achieve a converged solution at the end of each cycle.

3. Results

3.1 Velocity magnitude and pressure drop values using field plots

Fig. 10 and Fig. 11 below shows the velocity and pressure gradients across the channel flow model for different mesh sizes 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 with increase in the spatial resolution, both the velocity and pressure gradients offer improved visualisation of the viscous effects and pressure drop respectively, however, in the given scenario the difference is most notable when transitioning from the baseline mesh to mesh 1 as opposed to mesh 1 to mesh 2 or mesh 2 to mesh 3.

Fig. 10: Velocity Magnitude field plots across the domain for different mesh sizes

Fig. 11: Pressure field plots across the domain for different mesh sizes

3.2 Velocity magnitude and pressure drop values using line plots

The figures below i.e., Fig.12 and Fig. 13, shows the velocity magnitude profiles at the outlet of the channel flow model and pressure drop profile across the domain. It can be observed, in general, that for the velocity profiles, a similar range of their individual peak magnitude values was achieved across the different mesh sizes with the variations in the gradients being most notable when switching from the baseline mesh to mesh 1 and the difference becoming noticeably smaller from mesh 1 to mesh 2 and mesh 2 to mesh 3. The pressure drop profile too showed a remarkable increase in its maximum value when changing from baseline mesh to mesh 1 and the changes thereafter became increasingly smaller and more difficult to visualise in the same plot.

Fig. 12: Plot of Velocity Magnitude profiles at the outlet for different mesh sizes

Fig. 13: Plot of Pressure drop profiles across the domain for different mesh sizes

4. Discussion

Air, in the form of a mixture of two major contributing gas species `O_2` and `N_2` by mass fraction, enters at the inlet of the channel flow model of given cross-section (dimension representing length perpendicular to flow direction or height in case of two-dimensional space) with an almost uniform profile before encountering the boundary walls of the domain and as such, until this point, the flow could be considered as "inviscid" or the presence of viscous effects were mostly negligible. As air continues its travel downstream from the inlet, it's movement begins to get altered especially along the boundary surfaces of the channel flow model due to the "no slip" condition and the increasing presence of viscous effects causes the velocity of air flow at the boundary walls to arrive at a value of 0. The decelerating flow adjacent to the wall surfaces leads to a change in the velocity profile from uniform to a parabolic shape (characteristic of the laminar flow) with every layer from the wall to the center of the channel flow model increasing in velocity leading to the development of a boundary layer.

5.Conclusion

The numerical analysis for flow within a channel flow model under steady-state conditions using the default physical model in Converge Studio v3.0 and executed using the Cygwin64 terminal was well understood. It was realised that based on the values of spatial discretization, an improved resolution was offered into the visualisation of the boundary layer development that was well captured along the boundary surfaces of the channel flow model as well as a more accurate computation of the velocity and pressure drop gradients with regards to the flow properties. The CFD model can also be correspondingly investigated and validated as part of an extended scope towards evaluation of fully-developed flow profiles by enlarging the domain based on the estimated entrance length value for the given flow conditions which would also be the point wherein the boundary layer thickness would presumably have reached its maximum thickness.

References

[1] CONVERGE 3.0 Manual. © 2020 Convergent Science. All rights reserved. July 15, 2020.

[2] Engineers Edge, LLC, Viscosity of air, Dynamic and Kinematic, www.engineersedge.com, viewed 13 May 2020, <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.>

[3] The Engineering ToolBox, Speed of Sound in Air, www.engineeringtoolbox.com, viewed 13 May 2020, <https://www.engineeringtoolbox.com/air-speed-sound-d_603.html>