Week 6: Conjugate Heat Transfer Simulation
1. Introduction
Fluid flow as observed in varying applications for different industries typically takes place through pipes that have a circular cross-sectional area owing to their capabilities to withstand large pressure differences across the inner and outer regions without undergoing significant deformation as well as offering more heat transfer for considerably smaller pressure drop values when compared against pipes of different flow areas. Also, the effective rise in heat transfer coefficient is further catalysed when the fluid flow is characterised to be turbulent due to the associated intense mixing involved within it causing a rapid rise in friction force.
The aim of this challenge was to perform three-dimensional, super-cycling, Conjugate Heat Transfer (CHT) simulations using Realizable k-`epsilon` (RKE) turbulence model for a pipe flow geometry. The analysis involved computing results for different grid sizes and super-cycle stage intervals that would outline total cell count, velocity, temperature and pressure drop values across the domain, comparison between predicted and computed `y^+` values, and mesh quality attributes like skewness, non-orthogonality, and aspect ratio.
2. Method
The typical approach adopted to performing numerical simulations in CONVERGE CFD involves the following steps:
- Creating a new folder dedicated to running the concerned project inside a user-preferred directory.
- Setting up the relevant case file(s) in Converge Studio v3.0 that involved creation of the solid region and fluid volume domain, setting up the required application type, choosing material and its state (gas and solid) for which simulation is to be performed along with the requisite properties, defining simulation parameters such as solver types, run time variables, and convergence tolerances, defining regions, boundary conditions and initial conditions, specifying the physical models whose equations are intended to be solved, determining and allocating the base grid size values as well as any other refinements, and lastly selecting the variables for post-processing at the required time intervals along with any hidden parameters or UDF selection as part of advanced parameters.
- Saving the case file(s) and exporting the input files into the folder created for running the concerned project which typically includes 'boundary.in' for the various boundary conditions, 'material_state.dat' where material_state points to 'gas' and 'solid', 'initialize.in' specifying initial conditions, 'inputs.in' defining the various input conditions and file writing intervals, 'turbulence.in' containing information about the turbulence model selected and its associated constants, 'post.in' for the different post-processing variables, 'solver.in' describing the solver chosen and its equations, and 'surface.dat' containing datasets about the surface while other input files with respect to the state of the material, thermal properties and species may also be part of the exporting process and should be ensured that they are selected for the simulation to run.
- Running the simulation in a Cygwin64 terminal window.
- Post-converting 3D output data by navigating to the output folder created within the same project folder during the simulation run or from the 'Post-Processing 3D' button in Converge Studio v3.0.
The simulation run and post-converting processes were executed using Linux-type commands in a Cygwin64 terminal window as shown below:
i) Change directory to the project containing the relevant case file(s):
cd "PATH_TO_CONCERNED_PROJECT_FOLDER"
ii) Confirm the contents of the project for the exported input files:
ls
iii) Performing simulation run with the converge executable application on multiple processors (in this case, 4) using mpiexec executable file in the background (by adding ampersand operator) along with posting output into a logfile:
mpiexec.exe -n 4 converge.exe restricted </dev/null> logfile &
iv) Checking for simulation status with regards to completion or errors:
tail logfile
v) Post-converting the 3D output data files using the 'post_convert_30_msmpi_64' utility on multiple processors (in this case, 4) using mpiexec executable file from within the output folder of the project containing the exported input files:
cd output
mpiexec.exe -n 4 post_convert_30_msmpi_64
2.1 Modelling
The computational domain for simulating air flow movement through a pipe model subjected to heat flux on its external walls was modelled in Converge Studio v3.0. The pipe model was constructed for an inner flow radius of 0.015m and the solid regions at a radius of 0.02m with the extent selected at 0.31m in the direction of flow (or z-axis) while ensuring that the center of the model was positioned at the origin i.e., (0.0m, 0.0m, 0.0m), where the inlet boundary would also be situated, using the 'Cylinder' 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 need 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, SolidOuterPipeWall, SolidOuterPipeEdges, and SolidFluidInterface as illustrated in Fig. 1. Moreover, the normals from each of the triangular surfaces are required to point towards the fluid volume 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 normals are 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 observed that an error with "Nonmanifold edges" would be presented in the 'Diagnosis' window of the 'Diagnosis' dock at this time and they would highlight the edges especially in the region shared by the 'Inlet', 'Interface' and 'SolidOuterPipeEdges' boundaries on one side of the domain while 'Outlet', 'Interface' and 'SolidOuterPipeEdges' boundaries on the other as three triangles, one from each of the surfaces mentioned previously on either side of the domain, are being shared by a single edge and thus contributing to the "Nonmanifold edges" in the given domain. These problem triangles, however, would disappear when relevant boundary conditions are defined, in this case of type "INTERFACE", and therefore can be ignored so as to proceed with the required 'Case Setup' for the given numerical problem.
Fig. 1: Fluid Volume and Solid pipe regions for the pipe flow domain subjected to external heat flux with regards to the CHT simulation and Boundary Conditions 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.3.3.2 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 flow velocity from Reynolds number and Solution
The simulations on all meshes was performed using the Reynolds Averaged Navier-Stokes (RANS) models, in this case, the k-`epsilon` turbulence models. The use of k-`epsilon` turbulence models and more specifically the RKE turbulence model is attributed to the fact that the Reynolds number, based on the given flow conditions, lies in the turbulence flow region as well as representing wall attached flows for which the solution variables, like the velocity profile, was modelled for the log-law layer, however, using the Scalable wall function approach so as to allow for a certain degree of relaxation to model the same in order to yield satisfactory results based on the grid sizes selected, more of which has been highlighted in Section 2.3.3.2 Mesh Refinement approach. A brief description into the working of the RKE turbulence model is provided below:
2.3.1 RKE turbulence model
The RKE turbulence model proposed by Shih [2] is considered as an upgrade to the standard k-`epsilon` turbulence model while sharing the similarities of consisting two transport equations, one for the turbulent kinetic energy (k) to evaluate the energy in turbulence and the other for the turbulent dissipation (`epsilon`) to calculate the rate of dissipation of turbulent kinetic energy (k), however notable differences are present as mentioned below:
- The turbulent viscosity term is formulated differently in the RKE turbulence model.
- The dissipation rate `epsilon` in the transport equation has been modified and developed from the exact equation for transport of mean-square vorticity fluctuation.
The RKE turbulence model is able to satisfy some of the mathematical constraints of the Reynolds's stresses such as the positive values of normal stresses when the eddy viscosity definition is used in the Boussinesq relationship and Schwarz inequality for shear stresses even when the main strain rate is large which otherwise would be negative. These formulations in the RKE turbulence model have made them more popular to use over the other k-`epsilon` turbulence models especially when flow features include the presence of strong streamline curvature, vortices, rotation, boundary layers under the influence of strong and adverse pressure gradients, separation and recirculation.
2.3.1.1 Setup of the RKE turbulence model:
To use the RKE turbulence model in Converge Studio v3.0, we navigate to Physical Models option in the Case Setup tree followed by checking the 'Turbulence modeling' box to enable inclusion of turbulence terms in the solution of governing fluid equations and then selecting the Reynolds Averaged Navier-Stokes (RANS) item in the 'Turbulence modeling' window available under the 'Physical Models' option with 'Turbulence model' drop-down menu set to 'Realizable k-`epsilon`' and 'Near Wall treatment' switched to 'Scalable wall function' while all other variables/parameters and constants remained unaltered with respect to their original values.
2.3.2 Boundary Conditions and Monitoring Variables
The different boundary conditions for the CHT simulation of a pipe flow domain subjected to external heat flux is described in the following table:
|
Names |
Inlet |
SolidOuterPipeWall |
SolidOuterPipeEdges |
SolidFluidInterface |
Outlet |
|
Type |
INFLOW |
WALL |
WALL |
INTERFACE |
OUTFLOW |
|
Pressure Boundary Condition |
Zero normal gradient (NE) |
-- |
-- |
-- |
Specified Value (DI) |
|
Pressure Magnitude (Pa) |
-- |
-- |
-- |
-- |
101325.0 |
|
Velocity Boundary Condition |
Zero normal gradient (NE) |
-- |
-- |
-- |
Zero normal gradient (NE) |
|
Velocity Magnitude (m/s) |
(u v w) `-=` (0 0 3.6826) |
-- |
-- |
-- |
-- |
|
Temperature Boundary Condition |
Specified Value (DI) |
Heat Flux | Zero normal gradient (NE) |
Specified Value (DI) |
-- |
|
Temperature Magnitude (K) |
300 |
-- |
-- |
300 |
300 |
|
Flux (`W/m^2`) |
-- |
-10000 |
-- |
-- |
-- |
|
Temperature Backflow |
-- |
-- |
-- |
-- |
Specified Value (DI) |
|
Temperature Backflow Magnitude (K) |
-- |
-- |
-- |
-- |
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 |
|
Turbulent Kinetic Energy (tke) |
Intensity, 0.03836 fraction [3] |
-- |
-- |
-- |
Intensity, 0.03836 fraction [3] |
|
Turbulent Dissipation (eps) |
Length scale, 0.00105m [4] | -- | -- | -- | Length scale, 0.00105m [4] |
|
Wall motion type |
-- |
Stationary |
Stationary |
Stationary |
-- |
|
Surface Movement |
-- |
Fixed, Slip |
Fixed, Slip |
-- |
-- |
|
Contact Resistance (`(K*m^2)/W`) |
-- |
-- |
-- |
0 |
-- |
Table 2: Boundary Conditions for CHT simulation of a pipe flow domain subjected to external heat flux defined using the 'Boundary' sub-option of the 'Boundary Conditions' under the 'Case Setup' tree
| SolidFluidInterface | Forward side | Reverse side |
| Region | Fluid Volume | Solid Pipe |
| Velocity Boundary Condition | Law of wall | Slip |
| Temperature Boundary Condition | Law of wall; Coupled | Specified Value (DI); Coupled |
| Law of wall roughness parameters; Absolute Roughness (m) | 0 | -- |
| Law of wall roughness parameters; Roughness constant | 0.5 | -- |
| Heat model | O'Rourke and Amsden | -- |
| Turbulent Kinetic Energy (tke) Boundary Condition | Zero normal gradient (NE) | -- |
| Turbulent Dissipation (eps) Boundary Condition | Wall model | -- |
| Near wall treatment | Scalable wall function | Global |
Table 3: Configuration of "INTERFACE" type Boundary Condition for CHT simulation of a pipe flow domain subjected to external heat flux 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' and 'Solid 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. In the 'Solid simulation' window, 'aluminum.dat' item was selected under the 'Solid database' window when 'Predefined solids' option was clicked on followed by clicking on the 'Validate' and 'OK' buttons to confirm the same as the choice of material for the solid pipe walls and edges. The values of 'Turbulent Prandtl number' and 'Turbulent Schmidt number' were left unchanged at their default values of 0.9 and 0.78 respectively.
In regards to defining the boundary conditions for the solid region, the application of heat flux on the 'SolidOuterPipeWall' boundary was supplied at -10000`W/m^2` with the negative sign indicating transfer of heat into the domain while the 'SolidOuterPipeEdges' boundaries was set to 'Zero normal gradient (NE)' temperature boundary conditions to reflect these surfaces as adiabatic walls such that heat transfer is effectively prevented between the specified area of the solid region and the fluid volume. The 'SolidFluidInterface' boundary on the other hand represents the zone wherein the "INTERFACE" type boundary condition is defined to designate the distinction between the two regions viz., solid pipe and fluid volume respectively, and for boundary values needed to be mentioned such that the forward side represents the surface towards which the normal vectors are pointing to, in this case the fluid volume, while the reverse side would highlight the surface that is pointing away from the normal vectors i.e., the solid pipe and as such these boundary values would be similar to the 'Inlet' and 'Outlet' boundary conditions with the 'Heat model' set to the default option of 'O'Rourke and Amsden' for the former (forward side) and 'SolidOuterPipeWall' and 'SolidOuterPipeEdges' for the latter.
2.3.3 Determination of flow velocity from Reynolds number and Mesh Refinement approach
2.3.3.1 Determination of flow velocity from Reynolds number
Reynolds number (based on circular cross-section inlet 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
Substituting the value of characteristic length as inlet pipe diameter at 0.03m, and density and dynamic viscosity of air at 1atm or 101235 Pa and 300K (or 26.85°C) as 1.1766 `(kg)/(m^3)` and `1.857*10^-5 (kg)/(m.s)` [5] respectively into (Eq. 1), a corresponding flow velocity for Reynolds number of 7000 was estimated as 3.6826m/s. The reason for determining the flow velocity at the inlet of the given pipe domain was to supply the value obtained as an input for the velocity boundary condition at the same surface.
2.3.3.2 Mesh Refinement approach
In this numerical simulation, a baseline mesh for the given computational domain was created with element size of 0.004m, illustrated in Fig. 2, to get an initial assessment of the flow field. Thereafter, successive mesh refinements were carried out by utilising values for different mesh sizes such that the grid refinement factor defined as the ratio of coarse grid spatial size to the fine grid spatial size in the Journal of Fluids Engineering Editorial Policy [6] remained greater than the recommended value of 1.3, in this case reducing the base grid size for every mesh by 20%, leading to the different mesh sizes as reported in Table 5 along with the predicted `y^+` value which were then part of each individual case setup and was entered into corresponding columns representing the 'Base grid size' section of the 'Base grid' sub-option in the 'Grid Control' settings under the 'Case Setup' tree to compare the results obtained for the velocity magnitude and temperature quantities at the outlet boundary and pressure drop values across the entire domain.
Since, the Reynolds number for the given flow conditions within the pipe domain was required to be 7000, it implied 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 and is expressed below as:
`(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
It should be noted that the above expression is also valid for thermal boundary layer thickness as its formation is not determined due to thermal diffusion but the characteristics of the turbulent flow itself and thus, is dependent on the Reynolds number and not the Prandtl number. The point `x` depends on the flow characteristics developing within the pipe domain such that a presence of strong inviscid core occurs from the entrance of the domain to the location wherein the boundary layers from the opposite walls merge to form a fully developed flow due to the dominating viscous effects which is also referred to as the entrance length. The entrance length, in this case, is also the point at which the boundary layer thickness will have the maximum height before which the velocity reaches its free-stream value. Similarly, in this case, a thermal boundary layer also develops due to the temperature difference between the fluid and the walls of the pipe, subjected to an external heat flux, and as a consequence, the fluid layers close to the wall assume thermal equilibrium earlier than the fluid layers beyond it resulting in the formation of a temperature profile that extends until it reaches the center of the pipe where it becomes a thermally fully developed region marking this zone as the thermal entrance region and the length from the inlet as the thermal entry length. Entrance length, in the case of turbulent flow, from both a hydrodynamic and thermal perspective are almost identical when the domain length is greater than 10 times the inlet boundary diameter and can be expressed mathematically as:
`L_(h,text(turbulent))/D approx 10` (Eq. 3)
The corresponding value of entrance length from the expression in (Eq. 3) yielded a value of 0.3m. Utilising this value of entrance length `L_(h,text(turbulent))` in place of `D` in (Eq. 1), the local Reynolds number was calculated to be 69999.146. Therefore, the turbulent boundary layer thickness based on the entrance length definition and the Reynolds number at the same point i.e., `x = L_(h,text(turbulent))`, in (Eq. 2) was determined as 0.0122m. In order to capture the boundary layer development for the wall-bounded flow, the initial wall modelling approach of using a Standard wall function was chosen such that the first grid cell was in the range of 30 < `y^+` < 300 based on the RKE turbulence model selected. To perform calculations for `y^+` value, the values of wall shear stress and frictional velocity we required to be calculated as shown below:
Wall Shear Stress:
`tau_w=(1/2)⋅C_f⋅rho⋅(U_∞)^2`; where (Eq. 4)
`tau_w` is Wall shear stress,
`C_f` is Skin Friction coefficient,
`U_∞` is the freestream velocity, in this case the calculated value of inlet boundary velocity
Local skin friction coefficient (`C_f`) based on a semi-empirical expression can be related to the Reynolds number 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)
For local Reynolds number of 69999.146 with the previously mentioned values of dynamic viscosity, density and free-stream velocity, the skin-friction coefficient (`C_f`), wall shear stress (`tau_w`) and friction velocity (`U_tau`) values were calculated using (Eq. 5), (Eq. 4) and (Eq. 6) respectively as shown in the following table:
|
Local Reynolds Number |
Skin-friction Coefficient (`C_f`) |
Wall Shear Stress (`tau_w`)[kg m^-1 s^-2] |
Friction Velocity (`U_tau`)[m s^-1] |
|
69999.146 |
`6.379*10^(-3)` |
0.0509 |
0.208 |
Table 4: Skin-friction coefficient, wall shear stress and friction velocity values for the governing flow conditions of the given pipe domain subjected to external heat flux on the solid walls
Calculation of cell-centroid size (`y_P`) and First Cell Height (`Delta_y`):
`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
Simplifying (Eq. 7) to determine `y_P`, we get;
`y_P=(y^+*mu)/(U_tau*rho)` (Eq. 8)
Since, `y^+` is computed at the cell-centroid, the first cell height (`Delta_y`) can accordingly be mathematically expressed as:
`Delta_y = 2*y_P` (Eq. 9)
The initial approach adopted for evaluating the distance from the cell-centroid 'P' to the wall (`y_P`) and correspondingly the first cell height (`Delta_y`) was to target a `y^+` value of above 30, however, as the baseline mesh recommended for the given setup was defined using a grid size of 0.004m, the predicted `y^+` value achieved at the cell-centroid 'P' was approximately 26.358 using (Eq. 7) to (Eq. 9). This implied that further refinement would make it challenging to achieve and/or maintain the target `y^+` value to be greater than or at 30 respectively, and therefore to accommodate for the appropriate relaxation of the log-law assumption, the use of Scalable wall function was made over the Standard wall function.
|
Mesh |
dx(m) |
dy(m) |
dz(m) |
Predicted `y^+` value at `y_P` |
|
Baseline |
0.004 |
0.004 |
0.004 |
26.358 |
|
1 |
0.0032 |
0.0032 |
0.0032 |
21.086 |
|
2 |
0.00256 |
0.00256 |
0.00256 |
16.869 |
|
3 |
0.002048 |
0.002048 |
0.002048 |
13.495 |
|
4 |
0.0016384 |
0.0016384 |
0.0016384 |
10.796 |
Table 5: Mesh sizes representing spatial grid discretization values across each of the three-directional axes and their predicted `y^+` value at `y_P`
Baseline Mesh:
Fig. 2: Baseline Mesh (two-dimensional sliced view)
Mesh 1:
Fig. 3: Mesh 1 (two-dimensional sliced view)
Mesh 2:
Fig. 4: Mesh 2 (two-dimensional sliced view)
Mesh 3:
Fig. 5: Mesh 3 (two-dimensional sliced view)
Mesh 4:
Fig. 6: Mesh 4 (two-dimensional sliced view)
2.3.4 Determination of transient parameters for the fluid flow material across the pipe domain
The setup of simulation run time and time-step for the given numerical problem with respect to the fluid material was defined such that the 'Start time' in the 'Simulation Time parameters' window under the 'Case Setup' tree was set to 0s while the 'End Time' in the same window depended on the domain length and the flow velocity, in the direction of flow, such that the fluid medium passed through it at least 10 flow through times so as to achieve steady state conditions. Therefore, 'End Time' was evaluated as:
`End Time = (L/U)*10`; where (Eq. 10)
`L` is the length of the pipe determined as 0.31m, and
`U` is the fluid flow velocity determined from (Eq. 1) as 3.6826m/s
From (Eq. 10), the value of 'End Time' obtained is approximately equal to 0.842s. Time-step was decided based on the CFL number which is expressed below:
`CFL = (U*Delta t)/(Delta x)`; where (Eq. 11)
`Delta t` is time-step,
`Delta x` represents the smallest cell size, which in this case was the value of base grid size in either spatial direction.
In order to attain temporal accuracy and numerical stability, the CFL criteria was decided to be 1 and therefore, the value of time-step determined with this limiting value for different mesh sizes have been reported in Table 6.
| Mesh | `Delta t` (s) |
| Baseline | 0.001086 |
| 1 | 0.0008689 |
| 2 | 0.00069516 |
| 3 | 0.00055613 |
| 4 | 0.00044479 |
Table 6: Time-step (`Delta t`) values for different meshes of the pipe flow domain subjected to external heat flux
The above values of time-steps reflect for a domain whose smallest cell size across it would remain at the value of `Delta x`, as reported earlier, however, with the cartesian cut-cell meshing approach, the same cannot always be achieved and therefore, the time marching solution process or 'Time-step selection' chosen was of type 'Use variable time-step algorithm'. A variable time-step allows for internally calculated time-step that in turn depends on additional inputs including minimum and maximum time-step, and maximum convection CFL, among others, such that the resultant value determined based on the corresponding limiting conditions on the inputs would represent the maximum permissible time-step for the next iteration to be performed. The initial, minimum and maximum time-steps chosen for each setup is tabulated as seen in Table 7 below wherein the initial and minimum time-steps were chosen to be a value that would be a reasonably good start for the simulation, considering the given setup, while the maximum time-steps represents the `Delta t` values as calculated in Table 6.
| Mesh | Initial time-steps (s) | Minimum time-steps (s) | Maximum time-steps (s) |
| Baseline | 1e-7 | 1e-7 | 0.001086 |
| 1 | 1e-7 | 1e-7 | 0.0008689 |
| 2 | 1e-7 | 1e-7 | 0.00069516 |
| 3 | 1e-7 | 1e-7 | 0.00055613 |
| 4 | 1e-7 | 1e-7 | 0.00044479 |
Table 7: Initial, minimum and maximum time-steps values for different meshes of the pipe flow domain subjected to external heat flux
2.3.5 CHT simulation setup using Super-Cycling method
Super-Cycling in CONVERGE can be described as a method to expedite the time-dependent CHT problem to steady-state for the solid region in a given domain. This time control method allows to overcome the rather large time-scale requirement for the solid phase heat transfer by iterating between the fully-coupled transient and steady-state solvers in the following manner:
- The fluid and solid phase equations are solved at the same time using the transient solver, however, the solid heat transfer data is not stored from the beginning of the simulation but rather at the value mentioned in the 'Begin storing super-cycle data' parameter in the 'Super-Cycle modeling' window thereby allowing the fluid flow field to develop.
- At the arrival of the time specified in the 'Begin storing super-cycle data' parameter, the values for heat transfer coefficient (HTC) and near wall temperature for the cells lying at the interface are stored with a single value corresponding to the concerned time-step.
- The above step continues until the value passed in the 'Time length for each cycle stage' is reached.
- For the time equal to the summation of 'Begin storing super-cycle data' and 'Time length for each cycle stage' parameters, the solver for fluid flow is paused and time-averaged HTC and temperature for cells at the interface of the solid and fluid phases is computed for values stored in the previous step. This is followed by performing the solid heat transfer calculation using the time-averaged HTC and temperature for cells in the fluid phase near the interface as boundary conditions. There are two approaches that can be taken to solve the solid heat transfer equations viz., steady-state and transient. In the given numerical simulation, the steady-state method was preferred which obtains a solution in the following way:
- The solid heat transfer equations are solved until the value specified in 'Sensible internal energy (SIE) tolerance' is not reached, in this case 1e-7. Solid energy transfer calculations concluding at this point represent the steady-state solid temperature.
- The process repeats again from step i to step iii which consequently results in a drastic increase in fluid temperature based on the updated values of the steady-state solid temperature until the difference between the two phases does not cause the temperature of the former phase (fluid) to increase any further corresponding to the given flow rate. Also, new values of HTC and temperature for cells at the interface of the solid and fluid phases is stored at each time-step while the same quantities are utilised from the current super-cycle stage during averaging of those values.
- The moment time equals the summation of 'Begin storing super-cycle data' and 'Time length for each cycle stage' parameters, only this time the latter parameter will be one time more than its previous state i.e., 2*'Time length for each cycle stage', step iv repeats again.
The above procedure continues unless the 'End Time' parameter in the 'Simulation Time parameters' window has arrived. The 'Total number of stages' which represents the total number of stages to average data was retained at its default value of 1 while the 'Conduction CFL number' and 'Sensible internal energy (SIE) tolerance' values were altered to 100 and 1.4 respectively as was provided for the given numerical simulation.
Mesh and time statistics, peak velocity magnitude values near the outlet, and pressure drop across the domain for different meshes:
|
Mesh |
Cells |
Run Time [s] |
First Cell Height (`Delta_y`)[m] |
`y^+` (`y_P`) (Line 1) |
Mass Flow Averaged Temperature at Outlet [K] |
% var |
|
Baseline |
7776 | 404.2582 | 0.0029361 | 10.6193 | 423.371 |
-- |
|
1 |
13648 |
1645.278 |
0.00213647 |
8.12649 |
423.559 |
0.044 |
|
2 |
27256 |
2178.5668 |
0.00213566 |
8.55368 |
423.269 |
-0.0685 |
|
3 |
48000 |
4787.0263 |
0.00264846 |
11.0068 |
422.739 |
-0.1252 |
|
4 |
92680 |
14169.5278 |
0.00182647 |
7.63943 |
423.122 |
0.0906 |
Table 8: Mesh and time statistics, mass flow averaged temperature values at outlet, `y^+` at (`y_P`) for different meshes at Line 1 and percentage variation in values of quantities of interest between subsequent meshes
Mesh quality attributes for the different meshes:
|
Mesh |
Max Aspect Ratio |
Max Skewness |
Max Non-Orthogonality |
|
Baseline |
3.26826 | 1.17578e-38 | 30.7576 |
|
1 |
3.58671 |
1.17578e-38 |
40.1001 |
|
2 |
2.16748 |
1.17578e-38 |
29.4668 |
|
3 |
1.87707 |
1.17578e-38 |
4.15696 |
|
4 |
1.91062 |
1.17578e-38 |
3.47205 |
Table 9: Mesh quality attributes depicting the maximum values of aspect ratio, skewness and non-orthogonality for the different meshes of the pipe flow domain subjected to external heat flux
From the meshing statistics in Table 8, it can be seen that the computation of mass flow averaged temperature at the outlet boundary represented the least variation for mesh number 1 (or mesh 1) with 13648 cells, however, the percentage variations (% var) thereafter have been unsteady with increase in the same reported until mesh 3 and decrease for mesh 4. This resulted in averaging the mass flow averaged temperature values at the outlet determined by CONVERGE for all the meshes and evaluating the mesh whose difference would be the least from the mean leading to mesh 2 with 27256 cells which was accordingly deemed to be the mesh independent result, also in terms of the computational time and resources. The % var was calculated as the ratio of the difference of the calculated values in two consecutive meshes to the calculated value of the previous mesh e.g., ((Mass Flow Averaged Temperature at Outlet for Mesh 1 - Mass Flow Averaged Temperature at Outlet for Mesh 2)/Mass Flow Averaged Temperature at Outlet for Mesh 1). It can also be observed that the `y^+` values estimated for all the meshes at the first layer cell-centroid position were in the buffer region justifying the use of Scalable wall function as the near wall treatment modeling approach for the given problem, however, the differences between the initially calculated first cell height and the same computed by CONVERGE were quite apart and this may be attributed to the cells cut at the boundary using the cartesian cut-cell approach. Additionally, mesh quality attributes of all the different meshes prepared in Table 9 indicates that they lie in the range of acceptable value for skewness, as previously stated in Table 1, and improves with respect to aspect ratio and non-orthogonality from mesh 1 onwards for every successive refinement as the cells assume less stretching positions along a given direction and the maximum angle made between the vector connecting two cell centers and the normal of the shared face reduces respectively. An assessment was also carried out on the effects of choosing different super-cycle stage intervals both in terms of the values computed for Mass Flow Averaged Temperature at Outlet and the time taken to attain convergence of which the former has been highlighted in this section, as seen in Table 10, indicating that the overall effect on the quantity of interest in question is almost negligible while the run time decreases in a non-linear fashion possibly due to more regularly updated averaging leading to earlier convergence. Results in Section 3 dealing with visualisation of flow-field parameters using field plots will only be referred with regards to the mesh independent solution i.e., Mesh 2.
Impact of super-cycle interval on simulation time:
|
Mesh |
Super-Cycle Interval (s) |
Cells |
Run Time [s] |
Mass Flow Averaged Temperature at Outlet [K] |
% var |
|
2 |
0.05 |
27256 | 2178.5668 | 423.269 |
-- |
|
2 |
0.03 |
27256 |
2049.8641 |
423.137 |
-0.0312 |
|
2 |
0.02 |
27256 |
2051.1134 |
423.2 |
0.01489 |
|
2 |
0.01 |
27256 |
1920.1581 |
423.157 |
-0.01016 |
Table 10: Mesh and time statistics, mass flow averaged temperature values at outlet, `y^+` at (`y_P`) for different meshes at Line 1 and percentage variation in values of quantities of interest between subsequent meshes
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, and the heat transfer between the fluid and pipe walls of which the latter is subjected to external heat flux 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' (normally enabled by default) in the 'Geometry/Location' section of the 'Typical Variables/Geometry' tab, 'Turbulent Kinetic Energy', 'Turbulent Viscosity' and '`y^+`' in the 'Turbulence' section of the 'Combustion/Turbulence' tab, and adding 'Mass Fraction' under 'Species' section of the 'Species/Passive' tab for `O_2` and `N_2` species with respect to the concerned numerical analysis. Visualising the different flow variable quantities in two-dimensional space in ParaView required the creation of a slice of the given fluid domain by using the "Slice" tool in the 'Commons' toolbar typically located above the 'Pipeline Browser' window based on default layout settings and choosing 'X Normal' under 'Plane Parameters' section to align the plane normal to the x-axis direction such that the entire domain was viewed in the direction of flow followed by unchecking the 'Show Plane', 'Triangulate the slice' and 'Merge duplicated points in the slice' boxes for disabling the viewing of the plane borders, production of triangles in the output and overlapping data points respectively. Plotting the temperature, velocity and `y^+` profiles at the outlet of the pipe 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 11, 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 12.
It was ensured that the last time-step size of the simulation run was selected to view the temperature, velocity, `y^+` and pressure profile results. Additionally, temperature, velocity and pressure field plots were created to view the development of the flow field characteristics across the given domain. Legends, titles, and x and y labels were added for better understanding of the graphical representations.
|
Line |
x |
y |
z |
|
Point 1 |
0 |
-0.015 |
0.31 |
|
Point 2 |
0 |
0.015 |
0.31 |
Table 11: Line Coordinates for viewing the temperature, velocity and `y^+` profiles
|
Line |
x |
y |
z |
|
Point 1 |
0 |
0 |
0 |
|
Point 2 |
0 |
0 |
0.31 |
Table 12: Line Coordinates for viewing the pressure profiles across the domain
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.7 Convergence Level
In the given computational analysis, a transient simulation using first-order implicit Euler time advancement scheme was performed with respect to the fluid phase to accommodate for larger time-steps while also being numerically stable along with a Density-based solver, to account for the different density values of the two phases and especially during the process of heat transfer, of type modified Pressure Implicit with Splitting of Operator (PISO) which is an iterative pressure-velocity coupling algorithm involving a predictor step for predicting solution of the momentum equation from which a pressure equation is obtained and solved whose correction is then consequently reapplied to the momentum equation such that this process continues until the desired level of accuracy is achieved. At the end of the momentum predictor and first corrector step, remaining transport equations are solved in series. In other words, momentum and pressure are solved at the first PISO iteration which enables calculation for velocity terms required by other transport equations and at the end of each iteration in the PISO loop, convergence criterion is checked for which in the case of compressible flow (retained as the default option for 'Gas flow solver' in the 'Run parameters' window under the 'Case Setup' tree as the incompressible flow option created convergence issues) requires the density correction error term to be less than the value of PISO tolerance which remains unchanged for the given numerical simulation at its default value of 1e-03. The selection of solver types for the different equations within 'Solver parameters [Transient]' sub-option in the 'Simulation Parameters' settings under the 'Case Setup' tree included assigning SOR to Momentum, Pressure, Density, Energy, Turbulent Kinetic Energy (tke) and Turbulent dissipation rate (eps) equations. 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). The solver types, convergence tolerance limits, and minimum and maximum iterations defined for the equations corresponding to the RKE turbulence model equations for the transient state simulation and the monitor variable plots 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 |
SOR |
1e-08 |
2 |
5000 |
None |
|
Density |
SOR |
1e-04 |
2 |
20 |
None |
|
Energy |
SOR |
1e-04 |
2 |
20 |
None |
|
Turbulent Kinetic Energy (tke) |
SOR (pre-defined) |
1e-03 |
2 |
50 |
None (pre-defined) |
|
Turbulent dissipation rate (eps) |
SOR (pre-defined) |
1e-03 |
2 |
50 |
None (pre-defined) |
Table 13: Solver Types, Convergence Tolerance Values, Minimum and Maximum Iterations and Preconditioners of the different equations used in the RKE turbulence model
Fig. 7: Plots of monitor variables velocity, density, static pressure and temperature for transient CHT simulation using RKE turbulence model for baseline mesh
Fig. 8: Plots of monitor variables velocity, density, static pressure and temperature for transient CHT simulation using RKE turbulence model for mesh 1
Fig. 9: Plots of monitor variables velocity, density, static pressure and temperature for transient CHT simulation using RKE turbulence model for mesh 2
Fig. 10: Plots of monitor variables velocity, density, static pressure and temperature for transient CHT simulation using RKE turbulence model for mesh 3
Fig. 11: Plots of monitor variables velocity, density, static pressure and temperature for transient CHT simulation using RKE turbulence model for mesh 4
Fig. 12: Effect of super-cycle stage interval on the behaviour of Mass Flow Average Temperature plot at the Outlet
From Fig. 7, Fig. 8, Fig. 9, Fig. 10 and Fig. 11 we can notice that the monitor variables show a highly stable output well before they reach the end time of the simulation specified in the 'End time' parameter under the 'General' tab of 'Simulation Parameters' settings under the 'Case Setup' tree as 0.842s indicating that a steady-state solution has been obtained with the presence of sporadic oscillations in the average velocity and static pressure quantities especially in the first half of the same time-charts. Also, as previously mentioned, the effect of Mass Flow Averaged Temperature at the outlet boundary was assessed by varying the super-cycle stage intervals, as seen in Fig. 12, for the same yielding a trend that with lower values of the same parameter, the convergence for the quantity of interest was achieved earlier which simultaneously resulted in lower overall run time for the given simulation although this statistic cannot be truly accounted for as the number of processes running on the concerned machine was not determined to be the same during all simulations. It should be noted that a mass flow average of temperature and density quantities was considered over the area averaged values of the same quantities as the energy content is being aggregated over the entire boundary apart from them not being spatially conserved. The average velocity and static pressure quantities, on the other hand, are essentially spatially conserved and using either area averaging or mass flow averaging approaches would yield almost identical values.
3. Results
3.1 Velocity magnitude, pressure and temperature values using field plots
Fig. 13, Fig. 14 and Fig. 15 below shows the velocity magnitude, pressure and temperature values across the entire domain for mesh 2 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 progress in time from the initial state to the final state, all the flow-field quantities appear to have achieved steady-state conditions with a clear and distinct development of thin layers of velocity and temperature profiles in their respective plots from the inlet to the outlet due to the viscous and thermal boundary effects respectively on the fluid medium from the adjacent pipe walls.
Fig. 13: Velocity Magnitude field plot across the domain for mesh 2
Fig. 14: Pressure field plots across the domain for mesh 2
Fig. 15: Temperature field plot across the domain for mesh 2
3.2 Velocity magnitude, temperature and pressure and profiles using line plots
The figures below i.e., Fig.16, Fig. 17 and Fig. 18, shows the velocity magnitude and temperature profiles at the outlet of the domain and pressure drop profiles across the same respectively. 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, especially near the wall, which also reduces with every subsequent mesh refinement when compared to its predecessor barring the exception of mesh 3 whose velocity values computed near the wall are significantly higher than the remaining meshes. The temperature profiles too demonstrated similar trends as the velocity magnitude profiles with nearly identical peak magnitude values and the gradients near the wall reducing between the subsequent meshes other than the anomalous behaviour of mesh 3. Pressure drop profiles, on the other hand, showed an increase in its maximum value, however, the changes with respect to the same became smaller for every subsequent mesh.
Fig. 16: Plot of Velocity Magnitude profiles at the outlet for different mesh sizes
Fig. 17: Plot of temperature profiles at the outlet for different mesh sizes
Fig. 18: Plot of pressure drop profiles across the domain for different mesh sizes
4. Discussion
Air enters at the inlet of the given cross-section of the pipe flow domain subjected to external heat flux, containing a mixture of two major contributing gas species `O_2` and `N_2` by mass fraction, at uniform velocity but as turbulent flow which is typically identified by chaotic and unorganised mixing of fluid layers to constitute the formation of eddies that leads to acceleration in the transfer of momentum and energy and consequently higher values of friction, heat and mass transfer coefficients. The movement and intensity of eddies, however, decreases near the boundary walls as it travels downstream due to the boundary condition imposed on it and effectively is confined away from it such that the turbulent boundary layer formed experiences a steep velocity gradient close to the wall but changes at a slower rate in its core region such that the velocity towards the center of the pipe increases in order to maintain the mass flow rate. This abrupt change in velocity values close to the wall tends to create wall shear stress of greater magnitudes than laminar flow through the same domain with a more fuller shape of the fully developed velocity profile appearing for the turbulent flow. The other simultaneous effects occurring on the fluid flow is the development of a thermal boundary layer along the pipe due to differences in temperature between the flowing medium and the walls of the pipe, which is supplied energy in the form of external heat flux, causing convection heat transfer and similar to the velocity boundary layer (referred to as the turbulent boundary layer), the thickness of the thermal boundary layer continues to grow until it nearly reaches the center of the pipe and fills the entire region. Beyond the fully developed flow regions from both a hydrodynamic and thermal perspective, the velocity profile does not vary as such and neither does the shape of the temperature profile when constant surface heat flux is applied on the walls of the pipe, however, both the surface temperature and mean fluid temperature increases in an almost linear fashion in the direction of flow. Also, as the fluid flowing through the pipe domain is gaseous in nature, the effect of heating the same due to the external heat flux supplied on the walls of the pipe leads to an increase in the kinetic energy of the molecules with frequency in collisions amongst them increasing considerably causing rapid rise in the momentum transfer and thereby also the viscosity due to increased resistance in the flow with density accordingly decreasing.
5.Conclusion
The numerical analysis for CHT simulation of a pipe flow model using super-cycling with RKE turbulence model based on the governing flow conditions in Converge Studio v3.0 and executed using the Cygwin64 terminal was well understood. It was realised that using the super-cycling approach greatly improved the time-dependent CHT problem in terms of resolving the heat transfer characteristics by modeling the solid region to steady-state rather than a time-scale approach as well as choosing an appropriate value of super-cycle stage interval can greatly enhance the convergence behaviour of flow-field quantities such as temperature and essentially reduce the computational time. Moreover, based on the values of spatial discretization, an improved resolution was offered into the visualisation of the boundary layer development, from both a hydrodynamic and thermal perspective, that was well captured along the boundary surfaces of the pipe flow model and a more accurate computation of the velocity, temperature and pressure drop gradients with regards to the flow properties.
References
[1] CONVERGE 3.0 Manual. © 2020 Convergent Science. All rights reserved. July 15, 2020.
[2] Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z., and Zhu, J., 1995, “A New - Eddy-Viscosity Model for High Reynolds Number Turbulent Flows - Model Development and Validation.”, Computers Fluids, 24(3):227-238
[3] Basse, N.T. (2017), "Turbulence intensity and the friction factor for smooth- and rough-wall pipe flow", Fluids, vol. 2, 30.
[4] Flórez Orrego; et al. (2012). "Experimental and CFD study of a single phase cone-shaped helical coiled heat exchanger: an empirical correlation". Proceedings of ECOS 2012 – The 25th International Conference on Efficiency, Cost, Optimization, Simulation and Environmental Impact of Energy Systems, June 26–29, 2012, Perugia, Italy. ISBN 978-88-6655-322-9.
[5] Engineers Edge, LLC, Viscosity of air, Dynamic and Kinematic, www.engineersedge.com, viewed 04 August 2021, <https://www.engineersedge.com/physics/viscosity_of_air_dynamic_and_kinematic_14483.htm#:~:text=The%20viscosity%20of%20air%20depends,2%20%2Fs%20or%2014.8%20cSt.>
[6] Celik, I. , Ghia. U., Roache, P., Freitas, C., Coleman, H., Raad, P., "Procedure for Estimating and Reporting of Uncertainty Due to Discretization in CFD Applications," American Society of Mechanical Engineers Journal of Fluids Engineering, Vol. 130, July 2008.
[7] The Engineering ToolBox, Ratio of Specific Heat of Gases, www.engineeringtoolbox.com, viewed 04 August 2021, <https://www.engineeringtoolbox.com/specific-heat-ratio-d_608.html>
[8] The Engineering ToolBox, Universal and Individual Gas Constants, www.engineeringtoolbox.com, viewed 04 August 2021, <https://www.engineeringtoolbox.com/individual-universal-gas-constant-d_588.html>