Lab 6 - Continuity Constraint Algorithm

Objectives

This lab explores solutions to the 2-D, well-posed, laminar, incompressible Navier-Stokes CFD equations in the pressure projection form. The problem geometry and results compare directly to the previous lab which modeled the problem in 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 execution of the thermal cavity benchmark with a 32x32 element "best" mesh for Ra = 106 using three different Jacobian formulations:

The solutions to the thermal cavity benchmark via these three different Jacobian formulations will be compared using the following methods:

Problem Statement

The pressure projection formulation of the INS equations has been well-developed in the course lectures. For this problem statement, we will instead focus on the three important formulations of the Jacobian.

Full-Newton Jacobian


The Full-Newton Jacobian formulation requires the most terms and fully incoporates all four variables (U, V, TEMP, and PHI). Only three zeros occur: 2 since the temperature and velocity potential equations are not coupled and 1 since convection does not affect the horizontal velocity component. The resulting aPSE template pseudo-code is given below.

{J U1/U1} =  ()()()(;1)(B200)()
            +()()(U1+U2)(1020;0)(B3001)()
            +()()(U1+U2)(3040;0)(B3002)()
            +()()(U1)(1;0)(B3100)()
            +()()(U1)(3;0)(B3200)()
            +(REI)()()(1122;-1)(B211)()
            +(REI)()()(1324;-1)(B212)()
            +(REI)()()(1324;-1)(B221)()
            +(REI)()()(3344;-1)(B222)()

{J U1/U2} =  ()()(U1)(2;0)(B3100)()
            +()()(U1)(4;0)(B3200)()

{J U1/PHI} =  ()()()(1;0)(B201)()
             +()()()(3;0)(B202)()

{J U2/U1} =  ()()(U2)(1;0)(B3100)()
            +()()(U2)(3;0)(B3200)()

{J U2/U2} = ()()()(;1)(B200)()
           +()()(U1+U2)(1020;0)(B3001)()
           +()()(U1+U2)(3040;0)(B3002)()
           +()()(U2)(2;0)(B3100)()
           +()()(U2)(4;0)(B3200)()
           +(REI)()()(1122;-1)(B211)()
           +(REI)()()(1324;-1)(B212)()
           +(REI)()()(1324;-1)(B221)()
           +(REI)()()(3344;-1)(B222)()

{J U2/TEMP} = (-,GRSH,RE2I)()()(;1)(B200)()

{J U2/PHI} = ()()()(2;0)(B201)()
            +()()()(4;0)(B202)()

{J TEMP/U1} = ()()(TEMP)(1;0)(B3100)()
             +()()(TEMP)(3;0)(B3200)()

{J TEMP/U2} = ()()(TEMP)(2;0)(B3100)()
             +()()(TEMP)(4;0)(B3200)()

{J TEMP/TEMP} = ()()()(0;1)(B200)()
               +()()(U1+U2)(1020;0)(B3001)()
               +()()(U1+U2)(3040;0)(B3002)()
               +(PEI)()()(1122;-1)(B3011)()
               +(PEI)()()(3344;-1)(B3022)()
               +(PEI)()()(1324;-1)(B3021)()
               +(PEI)()()(1324;-1)(B3012)()
               +()(NUSL)()(0;1)(A200)()

{J PHI/U1} = ()()(EPMN)(1;0)(B3001)()
            +()()(EPMN)(3;0)(B3002)()
            +(-)()()(1;0)(A200)()

{J PHI/U2} = ()()(EPMN)(2;0)(B3001)()
            +()()(EPMN)(4;0)(B3002)()
            +()()()(2;0)(A200)()

{J PHI/PHI} = ()()(EPMN)(1122;-1)(B3011)()
             +()()(EPMN)(3344;-1)(B3022)()
             +()()(EPMN)(1324;-1)(B3021)()
             +()()(EPMN)(1324;-1)(B3012)()

Quasi-Newton Jacobian


The Quasi-Newton Jacobian uncouples the velocity potential from the other three variables. The resulting aPSE template pseudo-code is given below.

{J U1/U1} = ()()()(;1)(B200)()
           +()()(U1+U2)(1020;0)(B3001)()
           +()()(U1+U2)(3040;0)(B3002)()
           +()()(U1)(1;0)(B3100)()
           +()()(U1)(3;0)(B3200)()
           +(REI)()()(1122;-1)(B211)()
           +(REI)()()(1324;-1)(B212)()
           +(REI)()()(1324;-1)(B221)()
           +(REI)()()(3344;-1)(B222)()

{J U1/U2} = ()()(U1)(2;0)(B3100)()
           +()()(U1)(4;0)(B3200)()

{J U2/U1} = ()()(U2)(1;0)(B3100)()
           +()()(U2)(3;0)(B3200)()

{J U2/U2} = ()()()(;1)(B200)()
           +()()(U1+U2)(1020;0)(B3001)()
           +()()(U1+U2)(3040;0)(B3002)()
           +()()(U2)(2;0)(B3100)()
           +()()(U2)(4;0)(B3200)()
           +(REI)()()(1122;-1)(B211)()
           +(REI)()()(1324;-1)(B212)()
           +(REI)()()(1324;-1)(B221)()
           +(REI)()()(3344;-1)(B222)()

{J U2/TEMP} = (-,GRSH,RE2I)()()(;1)(B200)()

{J TEMP/U1} = ()()(TEMP)(1;0)(B3100)()
             +()()(TEMP)(3;0)(B3200)()

{J TEMP/U2} = ()()(TEMP)(2;0)(B3100)()
             +()()(TEMP)(4;0)(B3200)()

{J TEMP/TEMP} = ()()()(0;1)(B200)()
               +()()(U1+U2)(1020;0)(B3001)()
               +()()(U1+U2)(3040;0)(B3002)()
               +(PEI)()()(1122;-1)(B211)()
               +(PEI)()()(3344;-1)(B222)()
               +(PEI)()()(1324;-1)(B221)()
               +(PEI)()()(1324;-1)(B212)()
               +()(NUSL)()(0;1)(A200)()

Once the Quasi-Newton Jacobian is formulated and calculated, the uncoupled Phi Jacobian must then be calculated.

{J PHI/PHI} = ()()()(1122;-1)(B211)()
             +()()()(3344;-1)(B222)()
             +()()()(1324;-1)(B221)()
             +()()()(1324;-1)(B212)()

Tensor Product Quasi-Newton Jacobian


The Tensor Product Quasi-Newton Jacobian further simplifies the Quasi-Newton Jacobian. By recognizing that the quasi-Jacobian has the generic form of a convection-diffusion equation, we can separate the horizontal and vertical components with minimal resulting error. The resulting aPSE template pseudo-code is given below.

[J1] is in the x-direction:

{J U1/U1} = ()()()(;0)(A200)(DETJ)
           +()()(U1+U2)(1020;0)(A3001)()
           +()()(U1)(1;0)(A3100)()
           +(REI)()()(1122;-1)(A211)()

{J U1/U2} = ()()(U1)(2;0)(A3100)()

{J U2/U1} = ()()(U2)(1;0)(A3100)()

{J U2/U2} = ()()()(;0)(A200)(DETJ)
           +()()(U1+U2)(1020;0)(A3001)()
           +()()(U2)(2;0)(A3100)()
           +(REI)()()(1122;-1)(A211)()

{J U2/TEMP} = (-,GRSH,RE2I)()()(;0)(A200)(DETJ)

{J TEMP/U1} = ()()(TEMP)(1;0)(A3100)()

{J TEMP/U2} = ()()(TEMP)(2;0)(A3100)()

{J TEMP/TEMP} = ()()()(;0)(A200)(DETJ)
               +()()(U1+U2)(1020;0)(A3001)()
               +(PEI)()()(1122;-1)(A211)()

... and [J2] is in the y-direction:

{J U1/U1} = ()()()(;0)(A200)(DETJ)
           +()()(U1+U2)(3040;0)(A3001)()
           +()()(U1)(3;0)(A3100)()
           +(REI)()()(3344;-1)(A211)()

{J U1/U2} = ()()(U1)(4;0)(A3100)()

{J U2/U1} = ()()(U2)(3;0)(A3100)()

{J U2/U2} = ()()()(;0)(A200)(DETJ)
           +()()(U1+U2)(3040;0)(A3001)()
           +()()(U2)(4;0)(A3100)()
           +(REI)()()(3344;-1)(A211)()

{J U2/TEMP} = (RE2I)()()(;0)(A200)(DETJ)

{J TEMP/U1} = ()()(TEMP)(3;0)(A3100)()

{J TEMP/U2} = ()()(TEMP)(4;0)(A3100)()

{J TEMP/TEMP} = ()()()(;0)(A200)(DETJ)
               +()()(U1+U2)(3040;0)(A3001)()
               +(PEI)()()(3344;-1)(A211)()

Discussion of Results

Table 1 below lists the driving properties of the thermal cavity system. As seen, this problem matches that of the last lab where Ra = 106.

Table 1: Thermal cavity problem constants.
===========================
 NAME            VALUE     
___________________________
 UREF          1.59134E-02 
 TREF          6.74515E+01 
 PREF          2.11680E+03 
 PSIREF        1.59134E-02 
 OMREF         1.59134E-02 
 RREF          2.33985E-03 
 REYNO         1.00000E+02 
 GRASH         1.47077E+06 
 PECLET        6.79914E+01 
 PRANDL        6.79914E-01 
 RAYLE         1.00000E+06 
 BETA2         1.89625E-03 
===========================
	

Table 2 lists the integration factors used in the solution of the problem for the three Jacobian formulations. As seen, the Tensor Product Jacobian formulation allowed for the greatest timestep and greatest allowable change in Q. The Full-Newton Jacobian formulation was the most "finicky" and required the lowest timestep and allowable change in Q. It should be noted that the Full-Newton Jacobian formulation never converged although both previously mentioned variables were reduced by an order of magnitude. The Quasi-Newton Jacobian formulation was in the middle of the pack.

Table 2: Thermal cavity integration factors.
===========================================================================
JACOBIAN                               |    FULL        QUASI       TENSOR 
_______________________________________|___________________________________ 
 INITIAL_TIME                          |    0.          0.          0.     
 FINAL_TIME                            |    6.E6        6.E6        6.E6   
 PROBLEM_CONVERGENCE_CRITERIA          |    1.0E-4      1.0E-4      1.0E-4 
 MAXIMUM_CHANGE_IN_Q_(DQ)              |    0.05        0.05        0.2    
 INITIAL_TIME_STEP                     |    -.001       -.001       -.001  
 TIME_STEP_MULTIPLIER                  |    1.05        1.05        1.05   
 MAXIMUM_TIME_STEP                     |    0.2         0.2         0.5    
 MAXIMUM_NUMBER_OF_STEPS               | 5000        5000        5000      
 MAXIMUM_NUMBER_OF_ITERATIONS_PER_STEP |    8           8           4      
 ITERATION_CONVERGENCE_CRITERIA        |    1.0E-4      1.0E-4      1.0E-4 
 THETA_IMPLICITNESS_FACTOR             |    0.5         0.5         0.5    
===========================================================================
	

Table 3 shows the numerical differences in the final solution between the three Jacobian formulations. It would seem that the Tensor Product Jacobian formulation is the best since it delivers energy norms that are comparable with the Full-Newton Jacobian but in half the time. The Quasi-Newton Jacobian formulation looks to be the worst. Of course, all of this assumes that the Full Newton Jacobian solution is correct even if it did not converge. We will attend to that assumption shortly. It should be noted that all solutions utilize the "best" mesh discovered in the previous lab, where a nonuniorm progression ratio of 1.15/0.85 is used along the x-axis to compress elements against the left and right cavity walls.

Table 3: Thermal cavity solution characteristcs.
==========================================================================
 JACOBIAN        |     FULL NEWTON       QUASI-NEWTON             TENS0R  
_________________|________________________________________________________
 CONVERGENCE ?   |            NO                YES                YES    
 ITERATIONS      |          1720               3118               1046    
 FINAL TIME      |       287.340 s          583.926 s          461.984 s  
 COMPUTER TIME   |    0:59:57.63 s       1:25:00.00 s       0:25:24.85 s  
                 |                                                        
 U1E             |    3.355000E+01       2.351424E+02       3.473812E+01  
 U2E             |    3.279100E+02       2.052279E+02       3.265892E+02  
 PRSE            |    3.309500E+03       3.910968E+03       3.316177E+03  
 TMPE            |    6.444500E-02       1.007483E-01       6.431813E-02  
 PHIE            |        N/A            1.213663E+00       8.026520E-10  
==========================================================================
	

Full-Newton Jacobian


"The 500 Pound Sumo Wrestler"

Like a 500 pound sumo wrestler, the Full Newton Jacobian formulation was bulky and hard to handle! The author was unable to achieve problem convergence despite repeated attempts. This is quite puzzling, however, when one refers to the solution results as detailed in the following figures. Figure 1 seems to show very good convergence in the temperature norm. Perhaps one of the other variables was having difficulty with convergence? The answer is no. Upon inspection, the U1, U2, and Phi energy norms all show very good convergence with time. The author is at a loss to explain why aPSE could not coverge upon an answer.

Figures 2 and 3 give the final temperature and temperature energy norm solutions which compare favorably to the results seen in Lab 5. Figure 4 shows a very ugly final Phi solution. The oscillatory nature of the solution speaks to possible 2dX waves which might need artificial diffusion in the form of a TWS algorithm. In spite of this, however, the final pressure solution looks clean as seen in Figure 5.


Figure 1



Figure 2 - See a video of the temperature field solutions HERE.



Figure 3



Figure 4



Figure 5


Quasi-Newton Jacobian


"The Sumo Wrestler's Ugly Boy"

Like the ugly boy of the 500 pound sumo wrestler, the quasi-newton Jacobian formulation is easier to handle but doesn't look very good! Figure 6 shows good convergence in the temperature energy norms. Figure 7 shows a final temperature solution, however, that does not match that of the Full-Newton Jacobian formulation or the results of the previous lab. For these reasons, it is highly suspect. The final temperature energy norm field as seen in Figure 8 is much higher than before. The steady-state phi field seen in Figure 9 is even more confusing in that the oscillations have disappeared while very high potentials appear. Finally, the steady-state pressure field seen in Figure 10 also does not match that of the Full-Newton Jacobian case and is also suspect.


Figure 6



Figure 7 - See a video of the temperature field solutions HERE.



Figure 8



Figure 9



Figure 10


Tensor Product Quasi-Newton Jacobian


"The Sumo Wrestler's Cute Daughter"

The Tensor Product Jacobian formulation would be most akin to the sumo wrestler's cute daughter! The solution resembles that of its parent, the Full Newton Jacobian, with the added benefit that the code is much more petite and easier on the eyes. The results seen in Figures 11-15 very closely match Figures 1-5 from the Full-Newton Jacobian case. The same nasty oscillatory-like behavior is seen in the steady-state potential field, but this Jacobian implimentation still converged to a solution without any apparent problems.


Figure 11



Figure 12 - See a video of the temperature field solutions HERE.



Figure 13



Figure 14



Figure 15

Conclusion

In conclusion, the pressure projection method has been implimented to solve the 2-D, well-posed, laminar, incompressible Navier-Stokes CFD equations driven by thermal convection. These results compare directly to the solutions from the previous lab which solved the exact same problem formulation via a streamfunction-vorticity implimentation. Specifically, implimentation via the three different Jacobian formulations has resulting in the following conclusions:

Full Jacobian

Quasi-Newton Jacobian

Tensor Product Quasi-Newton Jacobian

As can be seen, the Tensor Product Jacobian formulation comes out the clear winner. This is very puzzling, however, since it is a further simplification of the Quasi-Newton Jacobian. One would think that if the Quasi-Newton Jacobian formulation delivers incorrect results, then the Tensor Product simplification would also fail to provide good results. The author does not have a good explanation for this apparent contradiction.

Finally, the pressure projection method is seen to produce equivalent results to the streamfunction-vorticity method while opening the door to the possibility of driving the system with an initial pressure field. Finally, the similarity of the results from both labs give the author more confidence in the efficacy of these solutions for the INS equations.