* 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 = 10^{6} using three different Jacobian formulations:

- Tensor Product Jacobian
- Quasi-Newton Jacobian
- Full Jacobian

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

- Steady-state Temperature Field
- Steady-state Temperature Energy Norm Field
- Steady-state Pressure Field
- Steady-state Potential (Phi) Field
- Temperature Energy Norm Iterative Convergence
- Maximum Allowable Time Step Size

**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.

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)()

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)()

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 = 10^{6}.

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 ==========================================================================

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

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

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**

- Provides accurate results
- Difficult to achieve steady-state convergence
- Steady-state potential solution exhibits 2dX oscillatory behavior

**Quasi-Newton Jacobian**

- Does converge to a steady-state solution, however,
- Provides inaccurate results and should not be trusted

**Tensor Product Quasi-Newton Jacobian**

- Provides accurate results
- Fast convergence to a steady-state solution
- Also exhibits suspect 2dX oscillatory behavior in steady-state potential field

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.