Performing Conjugate Heat Transfer (CHT) Simulation Using Super-Cycling in Converge CFD Software & Post-processing It In ParaView

 Abstract:- The work is focused on setting up a flow simulation through an aluminum pipe to visualize the conjugate heat transfer mechanism between the pipe walls & the fluid flowing inside the pipe. We are interested in performing a grid-dependency test & also study the effects of super-cycling on simulation time. Pictorial representation of the geometry is shown below in Fig 1.1. The geometry is designed in CAD software & imported to Converge as an STL file which is then edited using the various geometric editing tools available in Converge Studio. Newtonian fluid Air is considered as the working fluent. We will also be making use of boundary flagging techniques for assigning each side of the geometry to a particular boundary & these boundaries will later be assigned to a volumetric region. Next, we set up the simulation parameters using a steady-state solver. Once the setup is complete, we export these inputs files generated by Converge CFD & run it using CYGWIN to generate 3D post output files. These 3D output files are first converted into either Paraview vtk inline binary format or EnSite format for Paraview to read & generate the flow field simulation. The paper will also cover the explanation of conjugate heat transfer & Super cycling. The simulation will be run for 3 different time intervals for super-cycling stage & 3 different mesh grid size of dx=dy=dz=0.004m ; dx=dy=dz=0.003 & dx=dy=dz=0.002 will be used. We will also be comparing the results obtained from running the simulation obtained.

                 

Fig 1.1:- Geometry of the aluminum pipe duct

 

Theory                                                                                                                                                                                                                                                                                                 

                    The term conjugate heat transfer (CHT) is used to describe processes that involve variations of temperature within solids and fluids, due to thermal interaction between the solids and fluids. The exchange of thermal energy between the two physical bodies is called the study of Heat Transfer, the rate of transferred heat is directly proportional to the temperature difference between the bodies. A typical example is the heating or cooling of a solid object by the flow of air in which it is immersed, and some other example includes conduction through solids, free and forced convection in the gases/fluids, and thermal radiation.

Conjugate heat transfer corresponds with the combination of heat transfer in solids and heat transfer in fluids. In solids, conduction often dominates whereas, in fluids, convection usually dominates. Efficiently combining heat transfer in fluids and solids is the key to designing effective coolers, heaters, or heat exchangers. Forced convection is the most common way to achieve a high heat transfer rate.

Heat transfer in solids and heat transfer in fluids are combined in the majority of applications. This is because fluids flow around solids or between solid walls and because solids are usually immersed in a fluid.

Modes of heat transfer

• Conduction: diffusion of heat due to temperature gradients. A measure of the amount of conduction for a given gradient is heat conductivity.
• Convection: when the heat is carried away by moving fluid. The flow can either be caused by external influences, forced convection; or by buoyancy forces, natural convection. Convective heat transfer is tightly coupled to the fluid flow solution.
• Radiation: transfer of energy by electromagnetic waves between surfaces with different temperatures, separated by a medium that is at least partially transparent to the (infrared) radiation. Radiation is especially important at high temperatures, e.g. during combustion processes, but can also have a measurable effect at room temperatures.

An accurate description of heat transfer modes, material properties, flow regimes, and geometrical configurations enables the analysis of temperature fields and heat transfer.

Typical design problems involve the determination of:

• Overall heat transfer coefficient, e.g. for a car radiator.
• The highest (or lowest) temperature in a system, e.g. in a gas turbine, chemical reaction vessels, food ovens.
• Temperature distribution (related to thermal stress), e.g. in the walls of a spacecraft.
• Temperature response in time-dependent heating/cooling problems.

There are several challenges when solving the fluid and solid domains together. The largest challenge to overcome is the difference in timescales associated with the fluid and solid. In general, the heat transfer in the solid occurs more slowly than the convective and diffusive timescales dictated by the fluid flow in the cylinder. This disparity in timescales is problematic for CFD simulation because each cycle can be computationally expensive. CONVERGE offers an alternative approach called super-cycling for overcoming the time-scale disparity.

Super cycling is a method used by CONVERGE Studio in the case of conjugate heat transfer problems with a solid and liquid region. The basic idea of super-cycling is that the solver for the fluid domain is paused until the solver for the solid domain converges. This pausing is done in intervals that can be set by the user and will be explained in the case setup.

The Reynolds number (Re) is an important dimensionless quantity in fluid mechanics used to help predict flow patterns in different fluid flow situations. At low Reynolds numbers, flows tend to be dominated by laminar (sheet-like) flow, while at high Reynolds numbers turbulence results from differences in the fluid's speed and direction, which may sometimes intersect or even move counter to the overall direction of the flow eddy currents. These eddy currents begin to churn the flow, using up energy in the process, which for liquids increases the chances of cavitation. The Reynolds number has wide applications, ranging from liquid flow in a pipe to the passage of air over an aircraft wing. It is used to predict the transition from laminar & turbulent flow and is used in the scaling of similar but different-sized flow situations, such as between an aircraft model in a wind tunnel and the full-size version. Reynolds number can be defined as the ratio of internal forces to viscous forces within a fluid which is subjected to relative internal movement due to different fluid velocities. A region where these forces change behavior is known as a boundary layer, such as the bounding surface in the interior of a pipe. A similar effect is created by the introduction of a stream of high-velocity fluid into a low-velocity fluid, such as the hot gases emitted from a flame in the air. This relative movement generates fluid friction, which is a factor in developing turbulent flow. Counteracting this effect is the viscosity of the fluid, which tends to inhibit turbulence. The Reynolds number quantifies the relative importance of these two types of forces for given flow conditions and is a guide to when the turbulent flow will occur in a particular situation. Mathematical representation of Reynolds's number.

Where: Re = ρuD ÷ μ

• ρ is the density of the fluid (kg/m³)
• u is the flow velocity (m/s)
• D is the diameter of the tube (m)
• μ is the dynamic viscosity of the fluid (Pa·s or N·s/m² or kg/(m·s))
• ν is the kinematic viscosity of the fluid (m²/s). 

For flow in a pipe or tube, the Reynolds number is generally defined as

Where:    Re = ρuD_{H ÷ μ  =   uDH ÷ v  =  QDH ÷ vA  =  W ÷ DHμ} 

 

D_H is the diameter of the pipe (the inside diameter of the pipe is circular) (m),

Q is the volumetric flow rate  (m³/s),

A is the pipe's cross-sectional area (m²),

u is the mean velocity of the fluid (m/s),

 

μ (mu) is the dynamic viscosity of the fluid (Pa·s = N·s/m² = kg/(m·s)),                                                                

ν (nu) is the kinematic viscosity (ν = μ/ρ) (m²/s),

ρ is the density of the fluid (kg/m³),

W is the mass flow rate of the fluid (kg/s).                 

                                                                                                                 

 

In boundary layer flow over a flat plate, experiments confirm that, after a certain length of flow, a laminar boundary layer will become unstable and turbulent. This instability occurs across different scales and with different fluids, usually when Re_x ≈ 5×10⁵, where x is the distance from the leading edge of the flat plate, and the flow velocity is the freestream velocity of the fluid outside the boundary layer.

For flow in a pipe of diameter D, experimental observations show that for "fully developed" flow, laminar flow occurs when Re_D < 2300 and turbulent flow occurs when Re_D > 2900. At the lower end of this range, a continuous turbulent-flow will form, but only at a very long distance from the inlet of the pipe. The flow in between will begin to transition from laminar to turbulent and then back to laminar at irregular intervals, called intermittent flow. This is due to the different speeds and conditions of the fluid in different areas of the pipe's cross-section, depending on other factors such as pipe roughness and flow uniformity. Laminar flow tends to dominate in the fast-moving center of the pipe while slower-moving turbulent flow dominates near the wall. As the Reynolds number increases, the continuous turbulent flow moves closer to the inlet, and the intermittency in between increases, until the flow becomes fully turbulent at Re_D > 2900. This result is generalized to non-circular channels using the diameter, allowing a transition Reynolds number to be calculated for other shapes of a channel.

 

Wall functions are useful in telling the solver how to approach the solution near a wall. To understand this we need to understand the boundary layer & Y+.In physics and fluid mechanics, a boundary layer is the layer of fluid in the immediate vicinity of a bounding surface where the effects of viscosity are significant.

The behavior of the flow near the wall is a complicated phenomenon and to distinguish the different regions near the wall the concept of the wall has been formulated. Thus is a dimensionless quantity and is the distance from the wall measured in terms of viscous lengths. One of the reasons for the need of is to distinguish different regions near the wall or in the viscous region, however how exactly it helps in turbulence modeling or in general CFD modeling need to be well understood. If we intend to resolve the effects near the wall i.e., in the viscous sublayer then the size of the mesh size should be small and dense enough near the wall so that almost all the effects are captured. But in some cases, if the wall effects are negligible then there is an option of including semi-empirical formulae to bridge between the viscosity affected region and fully turbulent region, and in this case, the mesh need not be dense or small near the wall i.e., the coarse mesh would work.

Considering the first case i.e., near-wall modeling it is well-known that the mesh size should be small enough, however then the question follows is how small? Thus, here comes the concept of, and based on the value the first cell height can be calculated. The near-wall region is meshed using the calculated first cell height value with gradual growth in the mesh so that the effects are captured and avoiding overall heavy mesh count.

In a flow bounded by a wall, different scales and physical processes are dominant in the inner portion near the wall, and the outer portion approaching the free stream. These layers are typically known as the inner and outer layers. Considering the flow over a smooth flat plate the boundary layer can be distinguished into two types namely laminar boundary layer and turbulent boundary layer. Since we are dealing with the turbulent boundary layer let us not get into the laminar boundary layer. The typical boundary layer structure over a flat plate is shown below in Fig 1.2. In between the laminar and turbulent boundary layer there lies a transition region. Typically for flow over a flat plate, the transition usually occurs around.

The image attached below shows the boundary layer development. We first start in the laminar region, then transition to the turbulent region. Near the wall, due to the viscous forces, the velocity is dragged to zero and as the boundary layer develops, the velocity will reach free-stream velocity. This is termed as the no-slip condition. We can use the no-slip condition if we know the flow is completely laminar. This would mean that we are resolving the flow in the laminar region whatever the grid size. This, however, doesn’t hold if the flow is turbulent. Turbulent flows have a smaller laminar region near the wall and if we are to resolve the gradients at such a small region, we’d need a very small grid size. This is where the y+ value comes handy. The y+ value is a non-dimensional term that can be used to understand how coarse or fine our grid is. This can then be used to determine if wall functions need to be used or not. From the image attached below, we need to either be in the viscous sub-layer or the turbulent region. If we were in the viscous sub-layer, that would mean a y+ of less than 10 and if we were in the turbulent region, it would mean a y+ of greater than 30. This is very significant in our making the decision of using or not using a wall function.

 

   

Fig 1.2:- Boundary layer Structure.

From the above figure it can be seen that in turbulent boundary layer region flow near the wall has been analyzed in terms of three layers:

• The inner layer, or sublayer, where viscous shear dominates
• The outer layer, or defect layer, where large scale turbulent eddy shear dominates
• The overlapping layer, or log layer, where velocity profiles exhibit a logarithmic variation.

where "U0" is the free stream velocity and "ρ" is the boundary layer thickness, "y" is the vertical distance measured from the wall.

 

The significance of y+:

The y+ value denotes wherein the y+ vs u+ curve we are with respect to the cell size. It is undesirable to have a cell size that leads to a y+ of between 10 to 30 since the modeling of flows transitioning from laminar to turbulent is not easy and is still being researched. Due to this fact, it is better to have a y+ of less than 10 meaning we'd be at the laminar region or have a y+ of greater than 30 where we can go with a wall function and still achieve results. The y+ value can be used as a marker based on which we can decide on the base grid size and make a decision based on the computing power as to whether the grid needs to be refined further or coarsened. when considering complex geometry, refinement in the boundary layer may be required to ensure the desired y+ value is achieved. In such cases, re-mesh has to be done or else mesh adaption techniques have to be used to achieve the required value across the entire model.
Thus we have learned that the wall function approach and y+ value required is determined by the flow behavior and the turbulence model is used. If we have an attached flow, then generally we can use a Wall Function approach, which means a larger initial value, smaller overall mesh count, and faster run times. If one expects to have flow separation and knows that the accurate prediction of the separation point will be having an impact on the result, then he would be advised to resolve the boundary layer all the way to the wall with a finer mesh. Unfortunately, as the y+ value is dependent on the local fluid velocity which varies across the wall significantly for most industrial flow applications, it is not possible to know the exact y+ prior to running an initial simulation. Hence, it is important for one to get into the habit of checking the y+ values as part of his normal post-processing so that one can make sure to fall in the valid range for the flow physics and turbulence model selection.

 

 

GOVERNING EQUATIONS                                                                                                                                                                                                                                                                    

 

                                             The solution for a CFD simulation is obtained by solving the Navier-Stokes (NS) equations. The NS equations comprise of 5 equations, namely the continuity, momentum (3 equations), and energy equations. The continuity and momentum equations are solved for all flow problems with the energy equation being optional and can be used only when heat transfer is taking place to save computation time.

 

An important property of Converge is that any geometry created or exported is assumed to be made of triangles. In boundary flagging, we group these triangles to a particular ‘Boundary’ & assign these boundaries to a ‘Volumetric region’. Defining these boundary conditions is a fundamental step in any CFD simulation as this helps solve the NS equations. Converge creates and exports ‘input’ file of these complex governing equations which are then solved by CYGWIN for running the simulation. For our purpose, we will be assuming that there is no heat transfer and that the temperature is constant, leaving us with 4 equations now, the continuity, and 3 momentum equations, one for each axis direction. 

                                                                                                                     

                                                                                                                                         1. Geometry

The geometry of the cylindrical duct was designed using the geometry creating & editing tools available in the Converge studio. The dimensions of the larger cylinder channel will be of a radius of 0.02 & a smaller cylinder inside the large cylinder will be of radius 0.015 with a length of 0.2 in the z-direction which would add as "thickness" to the pipe. Once the geometry is designed, we perform our first ‘diagnostic’ check for probable errors. We will be getting several intersection errors. We fix this by using the ‘patch’ tool to close the front & back faces of the cylindrical pipe as shown in Fig 1.1. Once the geometry is fixed, we run a second ‘diagnostic’ run & find a new error which is the non-manifold error. This is caused due to the fact that we are running a conjugate heat transfer that requires an "interface", in Fig 1.2 we can observe the circle with a red boundary representing the inlet face of the pipe which has all its edge triangles shared by 3 other sides of a triangle. This is a common error when running a CHT problem. The interface is represented in Fig 1.3 with Green boundary. Next, we proceed to transform the ‘Normal’s’ inside the geometry.

 

Fig 1.1:- Cylindrical Pipe with closed Inlet & Outlet.

   

Fig 1.2:- A single side of a triangle is shared by 3 triangles which in turn creates the "Non-manifold Problems" error for every Conjugate Heat Transfer problem.

 

                                                         

                                              Fig 1.3:- The "green" region of the cylindrical pipe represents the "interface" between solid & fluid interactions.

 Every geometry will have a normal vector which is perpendicular to the geometry, for this problem we use the “Normal Toggle” option to check for the direction of the normal’s. If the normal’s are pointing outside of the geometry, then it’s essential to transform these normals to point inside the geometry, where the fluid flow will occur as shown in Fig 1.4.

 

     

                                                                                                        Fig 1.4:- Transforming Normals.

 

In Converge, all the surface information is stored in the form of 3 entities vertices, edges & triangles. Therefore, any geometry created or exported from other CAD software into Converge are converted into triangles & that is the fundamental entity in Converge, where every part of the geometry is assumed to be a triangle. The process of grouping these triangles into boundaries is known as boundary flagging. Each boundary is assigned with a distinct ID as shown in (Fig 1.5). For our geometry, we have 5 distinct surfaces Inlet, Outlet, Solid outer wall, Solid thickness & Interface.

                                                                                                                          Fig 1.5:-  Boundary flagging of the geometry.

 

                                                                                                                                             2. Case Setup

  Having edited the geometry, the next stage of the process is to set up the simulation parameters. This step involves specifying certain properties to capture the physics of the problem perfectly, such as materials used, simulation time parameters, solver parameters, initial & boundary conditions to solve the Navier Stokes equation which are PDE’s, defining body forces, type of fluid flow, species involved & grid size of the mesh.

We start setting up the time-based general flow simulation by first defining the Materials involved in the simulation. Since we are simulation fluid (air) flow inside an aluminum channel we choose ‘gas simulation' & 'solid simulation' with our choice of fluid for this project is ‘air’, hence we choose that as our pre-defined mixture. For our solid simulation, our predefined material would be 'Al'. For gas simulation parameters we will be using the Redlich-Kwong equation of state with a critical temperature of 133K & critical pressure of 3770000pa. For global transport parameters, we use the default values of ‘turbulent Prandtl number’ and ‘Schmidt number’ which is 0.90 & 0.78 respectively. Since the fluid flow is air our ‘species’ will be a mixture of Oxygen (O2) & Nitrogen (N2) with a chemical composition of 23% & 77% respectively.

 The next parameter for our case setup is the Simulation Parameters, which include the ‘run parameters’, ‘simulation time parameters’ & the ‘solver parameters’. For performing the Conjugate Heat Transfer problem inside a cylindrical duct, we will be running compressible gas flow using 'Transient' solver at full hydrodynamics simulation mode, since our geometry is simple & has no moving part & we also require to solve the NS equations. In addition to this, we will be using ‘density-based’ PISO Navier strokes solver. As for simulation time parameters, we will be running the simulation for 0.5 seconds with an initial & minimum time step as (1e-7).

As earlier, after designing the geometry, we grouped the triangles and assigned them to a particular boundary ID. Similarly, in Boundary Condition, we group the 5 boundaries & assign each of them to a volumetric region. To do so we add two volumetric regions from Initial Conditions & Events one for solid region & other for the fluid region & assign the volumetric region (Region 0) which is the fluid region to the 'Inlet' & 'Outlet' boundaries. The second region i.e (Region 1) which is the solid region was assigned t 'Solid outer wall' & 'solid thickness'. These volumetric regions are required to set up the initial conditions for solving PDE. Initial conditions are assigned at the volume whereas boundary conditions are given at boundaries. The Inlet & Outlet boundaries are inflow & outflow boundaries respectively where the air will be flowing inside the duct at speeds of 5 m/s in the z-axis with temperature boundary condition specified as 300 K. At the outlet, static pressure will be of 101325 pa  & velocity boundary condition will be following Neumann boundary conditions with a specified backflow temperature of 300k with Air being the species if backflow occurs. Solid outer wall & solid thickness are stationary following the 'slip' Velocity Boundary Condition. While the solid outer wall was added with an additional heat flux of -10000 (-ve because heat is flowing inside the boundary system). The Solid Thickness boundary is a stationary wall with fixed surface movement following the ‘slip’ conditions & a temperature boundary condition of ‘zero normal gradients’ which is Neumann boundary condition. For our Interface boundary condition, we split the boundary into 2 regions forward & reverse direction. For Forward direction, the Normal's will be pointing towards the fluid domain & will be following the law of the wall velocity boundary condition & Temperature boundary condition with Turbulent kinetic energy of zero normal gradients (zng). The Reverse direction will be assigned for the solid region. For Physical Models, we add Reynolds Averaged Navier Stokes (RANS) with RNG K-ε as our turbulence model with Von Karman Constant of 0.42 & 5.5 as our law of wall parameter. In addition, we need to use the Super cycling Method to speed up the entire simulation. We use a steady-state solver with variable start time & time length. We will be running 3 cases with different mesh grid size for 3 different supercycle time interval of [0.01, 0.02, 0.03]. The final parameter that we need to setup is defining the Mesh that is the Base grid size. For our case, we would be running the simulation for 3 different mesh grid sizes of dx=dy=dz=0.04, dx=dy=dz=0.03 & dx=dy=dz=0.02. 

 

                                                                                                                                 3. Post Processing

The function of Converge studio is to set up the simulation & then create several input files that are then exported to a particular folder. To run these input files we use CYGWIN, a command-line interface that reads these inputs files & solves the complex PDE’s of governing equation to generate several output files. These output files are then post-converted into 3D output files which are readable files for ParaView.

 

                                                                                                                                         4. Results

 

                                                   Case 1 : Base mesh grid size of dx=dy=dz=0.004 & Super Cycling stage interval of 0.01, 0.02, 0.03 

 

                                                  

                                    A. Pressure, Temperature & Velocity Contours For Super Cycling Stage Interval of 0.01

 

 

                     

                                                                                                                                                                                                                 

                                                                          

 

 

                                                            B. Pressure, Temperature & Velocity Contours For Super Cycling Stage Interval of 0.02           

 

                                                    

 

 

 

                                                                                                                                                                                                                

                                                         

 

 

                                                                       C. Pressure, Temperature & Velocity Contours For Super Cycling Stage Interval of 0.03.  

 

 

                                                                                                                                                                                                                

                                                           

 

                                             

 

 

 

                                                Case 2 : Base mesh grid size of dx=dy=dz=0.003 & Super Cycling stage interval of 0.01, 0.02, 0.03 

 

 

                                         A. Pressure, Temperature & Velocity Contours For Super Cycling Stage Interval of 0.01

 

   

   

                                                                                                                                                                                                             

                                                      

 

                                                         B. Pressure, Temperature & Velocity Contours For Super Cycling Stage Interval of 0.02

 

                                     

 

                                                                                                                                                                                                                        

                                                              

 

     

 

 

 

   

                                                             C. Pressure, Temperature & Velocity Contours For Super Cycling Stage Interval of 0.03

 

                                                                                                                                                                                                        

                                                  

 

 

       

 

                                                      Case 3 : Base mesh grid size of dx=dy=dz=0.002 & Super Cycling stage interval of 0.01, 0.02, 0.03 

 

 

 

 

                                                A. Pressure, Temperature & Velocity Contours For Super Cycling Stage Interval of 0.01

 

 

                  

                                                                                                                                                                                                                                    

                                                                              

                                                                    

 

 

                                                             B. Pressure, Temperature & Velocity Contours For Super Cycling Stage Interval of 0.02.

 

 

 

                                                                                                                                                                                                                       

                                                                 

                                                                            

 

 

                                                               C. Pressure, Temperature & Velocity Contours For Super Cycling Stage Interval of 0.03.

 

 

 

 

                                                                                                                                                                                                           

                                                      

 

                                                                                                                                          

 

 

 

 

                                                                                                            Observations

Successfully performed Conjugate Heat Transfer (CHT) by subjecting an aluminum pipe to a heat flux near the outer wall with air flowing inside the duct through the inlet. For case1, with a base mesh grid size of dx=dy=dz=0.004m, the total number of cells generated for the entire geometry is roughly 5000 cells, for case2, with a base mesh grid of dx=dy=dz=0.003m, the total number of cells generated is ~10,000 & for the final case i.e. case3, with a base mesh grid size of dx=dy=dz=0.002m, the total number of cells generated for the geometry is ~35,000, which indicates that the number of mesh cells generated in the geometry is highest when grid cell size is lowest. The smaller the mesh grid size the larger number of cells would be generated. From Temperature contour animation we observed that for all 3 cases the temperature increases near the walls of the cylinder duct, this is due to the fact that we applied "heat flux" near the walls of the geometry. We also observe that some layers of the fluid region near the walls also getting heated due to the heat transfer process. The temperature plots at the outlets show that for every 0.01 sec, 0.02 sec & 0.03 sec there is a sharp rise in the mean temperature of the solid domain for all 3 base grid cases at the outlet until it reaches a steady-state which indicates the super-cycling effect.