Week 6: Conjugate Heat Transfer Simulation

1. AIM:

• Simulate a simple Conjugate Heat Transfer simulation of flow through pipe.
• Explain the concept of CHT and Converge® Supercycling.
• Understand and propound upon the concept of Y+ and Wall-Functions in CFD.
• Perform a grid-independency test upto the best capabilites of our computational resources.
• Study the effect of supercyling parameters.
• Post process in Paraview.

2. GEOMETRY, BOUNDARY CONDITIONS AND INITIALIZATION:

• Geometry is an aluminum pipe of 0.2m length and 30mm I.D. with 5mm thickness.
• Inlet Boundary conditions: Dirichlet velocity of 3.5 m/s with temperature of inlet air of 300K. This was specified considering an inlet reynolds number of 7000.        
• Outlet Boundary conditions: At the outlet Dirichlet pressure of 101325Pa was specified.
• Wall: The wall of pipe were given a heat influx of 10000W/m2. Thickness of pipe is about 5mm.
• Timestepping and Cycles: Simulation was run for 0.5s with time minimum time step of 1E-09 and max of 1 second. The convection, diffusion and mach CFLs were kept at standard values of 1, 2 and 50 respectively.  The solution was continously monitored and CFL values were appropriately changed to accelerate the convergence rate.
• Turbulence Modelling: RANS based K-Epsilon was used with O'Rourke & Amsden wall treatments for heat transfer explained below with different turbulent wall functions depending on situations.
• Interface: The face which serves as a boundary between solid and fluid was marked as an interface. This boundary type is special in that it has two unique sets of boundary conditions, one for each side of the INTERFACE. Without any absolute roughness specification, Converge treats solid-solid and solid-fluid interfaces as if they were in perfect contact. Converge enforces a thermal continuity across this boundary, this means it has following constraints:

Along with this, law-of-the-wall boundary conditions were applied to the interface fluid region which    means for the velocity law-of-the-wall boundary condition and simulations with a k-ε turbulence model CONVERGE uses the Launder and Spalding [Launder and Spalding, 1974] wall model. This is given by following equation:

For same law-of-the-wall boundary for temperature, O'Rourke & Amsden wall treatments for heat transfer were selected. This is modelled as:

REGION INITIALIZATION:

As a good guess the domain was initialized with value of 101325Pa with 300K temperature, this will aid in faster convergence. Velocity was initialized with 3.5m/s along the pipe to develop the flow field faster.

Turbulence Kinetic Energy was initialized at 1 `m^2/s^2`and 100 `m^2/s^3`

What is Conjugate Heat Transfer? And why do we need Supercycling?

Conjugate heat transfer (CHT) is when heat transfer occurs simultaneously within and
between fluid and solid regions. In fluid domain, the  convection is the dominant mode of heat transfer where as in solids, it is conduction. This is exactly the crux of the problem. 

Ideally a CHT heat transfer would solve both solids and fluids heat transport equations on the same timescale. This might be perfect approach to solve when we are intersted in transient problems like heating of liner after 10 seconds of running the engine. But this process will take a very long time if you desire to have a steady state temperature profile of liner to design the cooling system maybe. The solver on its time-scale will result in fast heat transfer in the fluids but very slow transport in solids which is why you will see a long time for system to reach steady state. 

In order to counter this, Converge offers a novel approach to solve such problems without impractical simulation run times of normal solver. Supercycling iterates between fully-coupled transient & steady-state solvers with following algorithm:

 Supercycling parameters for grid independence test:

What is Y+ and why should CFD engineer pay any heed to it?

As most of us are aware of, when a fluid is flowing through a pipe or over a surface of solid there is a formation of something called as boundary layer due to friction and viscosity effects. This boundary layer starts from the no slip wall (where the velocity is zero) and slowly diffferent shear layers have increasing velocity as we move away in direction normal to the wall as shown below.

Apart from this development in normal direction, the boundary layer also slowly develops more and more turbulence as it moves on the surface in parellel direction as shown. This turbulent boundary layer isn't actually entirely turbulent, it can be split into 3 portions:

1) Viscous sublayer where viscous forces prevail and the layer is laminar 

2) The buffer region where there is a transitional phase between laminar viscous layer and the exterior log law region or turbulent region.

3) The outer portion of boundary layer is where there is a turbulent formulation, it is also termed as log law region as it is modelled with a logarithmic profile.

WALL FUNCTIONS:

In order to correctly capture this boundary layer, one has to discretize the near wall region into a lot of finer cells so that each cell can accurately RESOLVE the boundary layer, this is possible in lot of cases but in most realistic CFD simulations the geometry is just too large and complex to resolve it all. In order to avoid such computation costs, one can use coarser grids near the wall and MODEL the boundary layer using something called as WALL FUNCTIONS rather than RESOLVING it with much finer grids. Fundamentally speaking, WALL FUNCTIONS are nothing but estimation of a given transported variable near the wall. These functions are formulated by various experimental studies.

Now using wall functions requires knowledge of something called as Y+ (more on this later). Each and every turbulence model has some form of wall functions assigned with it in every CFD code. In case of Converge with K-Epsilon model we have following types of flow field:

1) Standard wall function: Usual Y+ range - 30 < Y+ < 300, only holds good for first cell in log law region.

2) Scalable wall function: Usual Y+ range - Y+ > 11.225, scales the standard wall functions from log law region to buffer region.

3) Non-equilibrium: Usual Y+ range - 30 < Y+ < 300, highly recommended for adverse pressure gradient flows near boundary since it compensates for assumption of local equilibrium, when the production of the turbulent kinetic energy is equal to the rate of its distruction, is no longer valid.

For mostly all the simulations we used Zero normal gradient boundary for K and used Scalable wall functions for Epsilon which are modelled in CONVERGE as :

Apart from these wall functions, discrete wall treatments are also available to model heat transfer : O'Rourke and Amsden (Amsden, 1997), Han and Reitz (Han and Reitz, 1997), Angelberger (Angelberger et al., 1997), and GruMo-UniMORE (Berni et al., 1993).

Back to Y+:

Another problem in CFD is same code will be used to simulate number of different fluids with different properties, different geometries and velocities. So a boundary layervelocity profile to normal distance of wall for one fluid might not be same to the other. Thus there was a need to generalize it across all fluids, in short we had to divide the normal distance with different fluid and flow field properties to non-dimensionalize it in a more general range. This non-dimensionalized distance is termed as a Y+. It has to be calculated to estimate the boundary layer profile and then a decision has to be made on which region of boundary layer would we be placing the first cell in order for it to be compatible with particular turbulence model and wall function. The following graph shows 3 layers clearly

So now in our case of K-Epsilon model with given flow field properties we can calculate the first cell height with formula for each grid we use in mesh independence study. This is very imporatnt because any heat transfer from interface will be carried by the convection of flow and hence it can give a lot of variation in the experiment if this isn't implemented correctly. 

NOTE: These calculations are based on the empirical experiment of Blasius boundary layer over the flat plate, although it is pretty accurate ball park in almost all the other experiments we run we will never hit the exact value at each point when we plot the contours of the same. We will verify the same when we plot Y+ contours.

The above formula can be used to ball park the Y+ values you will be seeing near the wall but are not always 100% on the spot as we already informed they are based on Blasius boundary layer experimentations. In our case we ran initial solution and from the developed flow field plot we obtained Y+ and then re-run the solution with correct set of Wall functions. 

MESH INDEPENDENCE STUDY:

The study was run for multiple mesh sizes: 

• CASE 1: dx = 5e-3, dy=5e-3, dz=5e-3. Scaleable wall functions used

• CASE 2: dx = 4e-3, dy=4e-3, dz=4e-3. Scaleable wall functions used

• CASE 3: dx = 3e-3, dy=3e-3, dz=3e-3. Scaleable wall functions used

• CASE 4: dx = 2e-3, dy=2e-3, dz=2e-3. Scaleable wall functions used

As seen below, with refinement of mesh the temperature of fluid is pretty much packed closely, but the solid temperature is still converging. By the looks of it seems the temperature will still take a lot of grid refinement till it converges independently of grid. With available computational resources this was the best possible refined grid we could solve but in order to get further accuracy it has to be solved on a further refined grid. Also in most cases the Y+ range was good for Scaleable wall function and K-Epsilon turbulence model but for the finest grid of 0.002 and below it would be recommended to use K-OMega SST model since the Y+ values will be out of range of available wall functions for K-Epsilon model. 

 CELL COUNTS FOR GRIDS:

EFFECT OF SUPERCYCLING STAGE-TIME:

With K-epsilon model the grid size of 0.02 with scaleable wall functions was the one which was most refined grid we had so it was selected for studying the supercycling parameter changes, parameters were also varied in 3 different values from 0.01 to 0.03.

All were run and compared on the same graph to gain better understanding:

 By the look of graphical plots above and the simulation time from verbose output, it looks like with each increased supercycling interval the solution overall converges faster but with respect to the actualy flow time, it takes more time to converge as expected and shown in the graph above.

VELOCITY AND TEMPERATURE CONTOUR VIDEO: