Lab 5 - Thermal Cavity, Quasi and Newton Algorithm

Objectives

This lab explores solutions to the 2-D, well-posed, laminar, incompressible Navier-Stokes CFD equations in the streamfunction-vorticity form. Specifically, the problem geometry is characterized by an enclosed cavity driven only by thermal convection. The accomplishment of this objective occurs via the following proceedure:

Problem Statement

For the case of 2-D, incompressible, laminar, unsteady flow, we can enforce the continuity equation via introduction of the vorticity and the streamfunction. In so doing, the velocities and pressure gradient "disappear" and the system reduces to two equations on the vorticity and the streamfunction. As a change from the previous lab, however, we now modify the vorticity equation with a buoyancy term and we add an energy equation. The problem geometry is an insulated cavity with a hot wall and a cold wall. Thus, the streamline is zero along the cavity boundary. The vorticity boundary condition depends upon the modified Taylor series as defined in the course lectures with the noted simplification that the velocity at the cavity boundary is everywhere zero. The temperature boundary condition highlights the opposite hot and cold walls with thermal insulation on the other two surfaces.


The streamline equation is unmodified from the previous lab.

{FPSI} =  ()()()(1122;-1)(B211)(PSI)
         +()()()(3344;-1)(B222)(PSI)
         +()()()(1324;-1)(B221)(PSI)
         +()()()(1324;-1)(B212)(PSI)
         +(-)()()(0;1)(B200)(OMGA)

Since the vorticity and temperature equations involve an unsteady time term, we use the well-developed Galerkin Weak Statement modified with the Euler implicit Taylor series as shown.

Now, the vorticity equation and pseudo template may be developed. It is noticed that this development is identical to the previous lab except for the addition of a buoyancy term.

{FOMG} = ()()()(0;1)(B200)(-OMGA)

         +(REI)()()(1122;-1)(B211)(OMGA)
         +(REI)()()(3344;-1)(B222)(OMGA)
         +(REI)()()(1324;-1)(B221)(OMGA)
         +(REI)()()(1324;-1)(B212)(OMGA)

         +( )()(PSI)(23;-1)(B3102)(OMGA)
         +(-)()(PSI)(14;-1)(B3102)(OMGA)
         +( )()(PSI)(14;-1)(B3201)(OMGA)
         +(-)()(PSI)(23;-1)(B3201)(OMGA)

         +(-,GRSH,RE2I)()()(1;0)(B201)(TEMP)
         +(-,GRSH,RE2I)()()(3;0)(B202)(TEMP)

The temperature equation development closely follows that of the vorticity equation.

{FTEMP} = ()()()(0;1)(B200)(-TEMP)

          +(PEI)()()(1122;-1)(B211)(TEMP)
          +(PEI)()()(3344;-1)(B222)(TEMP)
          +(PEI)()()(1324;-1)(B221)(TEMP)
          +(PEI)()()(1324;-1)(B212)(TEMP)

          +( )()(PSI)(23;-1)(B3102)(TEMP)
          +(-)()(PSI)(14;-1)(B3102)(TEMP)
          +( )()(PSI)(14;-1)(B3201)(TEMP)
          +(-)()(PSI)(23;-1)(B3201)(TEMP)

With the vorticity, streamline, and temperature equations defined in their aPSE pseudo template code, we now move to the solution algorithm. The nonlinearities seen in the problem are best handled via the tried-and-tested Newton algorithm with an added twist. The "twist" comes in our handling of the Jacobian matrix. Instead of calculating a 3x3 matrix at every iteration, we will instead take advantage of the two zeros and uncouple the temperature from the vorticity and streamfunction. In this way the computer algorithm can update the temperature via its own quasi-Newton algorithm. The new temperature iteration is then applied to the next iteration of the vorticity and streamfunction equations as detailed in the previous lab.

{J OMG/OMG} = ()()()(0;1)(B200)()
              +(REI)()()(1122;-1)(B211)()
              +(REI)()()(3344;-1)(B222)()
              +(REI)()()(1324;-1)(B221)()
              +(REI)()()(1324;-1)(B212)()
              +( )()(PSI)(23;-1)(B3102)()
              +(-)()(PSI)(14;-1)(B3102)()
              +( )()(PSI)(14;-1)(B3201)()
              +(-)()(PSI)(23;-1)(B3201)()

{J OMG/PSI} = ( )()(OMGA)(23;-1)(B3201)()
             +(-)()(OMGA)(14;-1)(B3201)()
             +( )()(OMGA)(14;-1)(B3102)()
             +(-)()(OMGA)(23;-1)(B3102)()

{J PSI/OMG} = (-)()()(0;1)(B200)()

{J PSI/PSI} =  ()()()(1122;-1)(B211)()
              +()()()(3344;-1)(B222)()
              +()()()(1324;-1)(B221)()
              +()()()(1324;-1)(B212)()
{J TEMP/TEMP} = ()()()(0;1)(B200)()

                +(PEI)()()(1122;-1)(B211)()
                +(PEI)()()(3344;-1)(B222)()
                +(PEI)()()(1324;-1)(B221)()
                +(PEI)()()(1324;-1)(B212)()

                +( )()(PSI)(23;-1)(B3102)()
                +(-)()(PSI)(14;-1)(B3102)()
                +( )()(PSI)(14;-1)(B3201)()
                +(-)()(PSI)(23;-1)(B3201)()

Once a solution is obtained, we must perform a post-processing operation to revert the solution back into the two velocity components and the pressure gradient. The vorticity equation gives a solution for the two velocity components as follows:

{U1} = (-)()()(2;0)(B201)(OMGA)
      +(-)()()(4;0)(B202)(OMGA)

{U2} = ()()()(1;0)(B201)(OMGA)
      +()()()(3;0)(B202)(OMGA)

The streamfunction equation also gives a solution for the two velocity components:

{UPSI} = (-)()()(2;0)(B201)(PSI)
        +(-)()()(4;0)(B202)(PSI)
 
{VPSI} = ()()()(1;0)(B201)(PSI)
        +()()()(3;0)(B202)(PSI)

Finally, we may solve for the pressure gradient in the solution, remembering that the new temperature dependence will add to the code developed in the previous lab.

{PRES} =  ()(U1)()(11;-1)(B211)(U1)
         +()(U2)()(12;-1)(B211)(U1)
         +()(U1)()(21;-1)(B211)(U2)
         +()(U2)()(22;-1)(B211)(U2)
         +()(U1)()(13;-1)(B212)(U1)
         +()(U2)()(14;-1)(B212)(U1)
         +()(U1)()(23;-1)(B212)(U2)
         +()(U2)()(24;-1)(B212)(U2)
         +()(U1)()(31;-1)(B221)(U1)
         +()(U2)()(32;-1)(B221)(U1)
         +()(U1)()(41;-1)(B221)(U2)
         +()(U2)()(42;-1)(B221)(U2)
         +()(U1)()(33;-1)(B222)(U1)
         +()(U2)()(34;-1)(B222)(U1)
         +()(U1)()(43;-1)(B222)(U2)
         +()(U2)()(44;-1)(B222)(U2)

         +(-,GRSH,RE2I)()()(2;0)(B210)(TEMP)
         +(-,GRSH,RE2I)()()(4;0)(B220)(TEMP)

Discussion of Results

The first order of business in a discussion of the solution is the identification of the properties defining the problem. Table 1 below lists the constants used to define the problem. Regarding the units, any question marks indicate an assumption where English units are used.

Table 1: Thermal cavity problem constants.
 =====================================================
 NAME            VALUE                UNITS          
_____________________________________________________
 GRAVCN         32.17398              ft/s^2         
 TMIN            0                    F              
 TREF           67.4515               F              
 TABS          527.0515               F              
 ALC             1.0                  ft             
 VISCK           1.591337E-4          ?              
 CONDCT          4.175889E-6          Btu/s-ft-F (?) 
 SPHEAT          0.237                Btu/lbm-F  (?) 
 GSCNST          4.972E4              Btu/lbm-R  (?) 
 PREF            2.1168E3             lb/ft^2    (?) 
=====================================================
	

The Prandtl number is an important parameter used to characterize the thermal nature of the flow. After hunting around in the "hooks" files, the following code was discovered for the calculation of the Prandtl number. The resulting value calculated for every solution in this lab exercise was Pr = 0.679914.

RREF = PREF/(GSCNST*TABS)
VISC = VISCK*RREF
PRANDL  = VISC*SPHEAT*GRAVCN/CONDCT

In addition, the following code supplies the definition for the calculation of the Grashoff an Rayleigh numbers.

BETA2 = 1/(TABS + DELTEM/2)
GRASH = BETA2*GRAVCN*DELTEM*ALC**3/VISCK**2
RAYLE = GRASH*PRANDL

With the Rayleigh number properly defined, we can now proceed to calculate the problem variables. With the constants defined in Table 1, a given Rayleigh number allows for solution of the resulting maximum temperature, TMAX. With TMAX defined, we can then calculate the intrinsic velocity reference scale using a work balance as shown below.

UREF = (2*GSCNST*BETA2*ALC*DELTEM)^0.5

Table 2 lists the appropriate maximum temperatures and reference velocities for the given range of Rayleigh numbers.

Table 2: Thermal cavity variable calculations.
=====================================
 RAYLE       TMAX           UREF     
              (F)          (ft/s)    
_____________________________________
  10^3    6.101250E-4    8.630782E-3 
  10^4    6.101282E-3    2.729293E-2 
  10^5    6.101600E-2    8.630782E-2 
  10^6    6.104780E-1    2.729293E-1 
=====================================
	

Finally, we also wish to test the solution when using a different non-dimensional velocity scale. A thermal balance results in a reference velocity independent of the Rayleigh number.

UREF = VISCK/ALC/PRANDL = 2.340497E-4

Use of this velocity then results in a Peclet value of unity as requested by the course instructor.



Figure 1 below shows the resulting steady-state temperature solution for Ra = 103. As expected, this case with the lowest hot wall temperature results in mild, vertical, isotherms. Figure 2 shows the resulting steady-state temperature solution for Ra = 104. The higher Rayleigh number indicates a higher hot wall temperature. The higher buoyancy forces then begin to induce a rotation as seen by the vertical to horizontal isotherms. The situation is even more pronounced in Figure 3 where Ra = 105. Finally, in Figure 6 at Ra = 106 we see a sharp flowfield character. Any higher Rayleigh number would most likely trip the flow into turbulence.


Figure 1



Figure 2



Figure 3



Figure 4


In the previous 4 temperature solutions, it is seen that the temperature isotherms "bunch up" against the X = 0 and X = 1 ft walls. This immediately suggests that a non-uniform progression ratio mesh might provide a more robust solution. Figures 5-7 below detail a "best" mesh study over a range of symmetric progression ratios. It should be noted that a symmetric progression ratio of 1.1 indicates a progression ratio of 1.1 for 0 < X1 < 0.5 ft along with an equivalent progression ratio of 0.9 for 0.5 < X1 < 1 ft. The results in the streamline energy norm especially shows a clearly optimum progression ratio of 1.15 / 0.85 which is declared the "best" mesh.


Figure 5



Figure 6



Figure 7


To highlight the applicability of the non-uniform "best" mesh, Figure 8 shows the vorticity energy norm for the Ra = 105 solution. Figure 9 below then shows the vorticity energy norm for the same Rayleigh number using the "best" 1.15/0.85 progression ratio mesh. As seen, the peaks at the walls are reduced by a factor of two. Thus, the non-uniform mesh results in a more smooth energy norm pattern which indicates a more optimum solution!


Figure 8



Figure 9


Table 3 allows for some final comparisons between other solution methods. As seen, three separate solutions were obtained for the case of Ra = 105. An implimentation of the thermal balance definition for the reference velocity resulted in the same rate of solution convergence. The big surprise, however, is that the resulting global vorticity energy norm is 7 orders of magnitude greater than the work balance method. Clearly the work balance is the more optimal choice. Also, a full Newton iteration solution was obtained utilizing the full 3X3 Jacobian matrix. As seen, an equivalent solution is obtained, as seen in the vorticity norm. The solution, required less than half the number of time steps, however, and also required half as much computing time. This was a surprise. The author expected a trade-off between convergence rate and solution time. With the results as seen below, why should a quasi-Newton algorithm be used at all?

Table 2: Progression ratio optimization.
========================================================================
 PROCESS     U scale    RAYLE     STEPS   SOLUTION TIME     OMGE FINAL  
________________________________________________________________________
  Quasi         Work     10^3       15       0:10.63       2.622588E-01 
  Quasi         Work     10^4       28       0:23.15       5.247588E-01 
  Quasi         Work     10^5      147       3:25.00       1.060679E+00 
  Quasi      Thermal     10^5      147       3:26.00       5.318758E+07 
  Newton        Work     10^5       51       1:19.24       1.060679E+00 
========================================================================
	


Conclusion

In conclusion, this lab has expanded the scope of the previous lab with the addition of a thermally driven flows. The resulting solutions are much more benign as compared to the previous lab. This is attributed to the fact that the small Rayleigh numbers result in relatively small Reynolds number which indicates highly laminar flow. Also, there is no moving plate to introduce singularities along the boundary. The appropriateness of the work balance definition of the reference velocity has been seen as well as the existence of a "best" mesh that provides finer element spacing against the cavity walls where the isotherms are most dense. Finally, the robustness of the full Newton iterative method is seen which breeds the question, why didn't we use it in the first place?