Lab 4 - Driven Cavity Problem

Objectives

This lab explores solutions to the 2-D, incompressible Navier-Stokes CFD equations. Specifically, this lab seeks to solve well-posed, isothermal, laminar, incompressible N-S equations in the streamfunction-vorticity form. The accomplishment of this objective occurs via the following proceedure:

Problem Statement

For the case of 2-D, isothermal, 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.

Following the usual procedure, the unsteady time term is reexpressed in terms of a Taylor series expansion with a Euler-implicitness term.

The vorticity equation may now be expressed in template pseudo-code.

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

The streamfunction equation is similarly derived in the template pseudo-cod where the wall velocity condition is expressed everywhere as zero.

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

Again as expected, the non-linear nature of the equations requires solution via an iterative method. We choose the time honored Newton iterative method which then requires the calculation of the Jacobian. The procedure is as follows:

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

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 as follows:

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

As one final note, it is expected that the solution will experience instabilities due to dispersion error. A modified Taylor Weak Statement (TWSh) can provide added diffusion to help stabalize the system. For a steady-state solution, the following modifications are made to the vorticity equation:

There are three different options to handle this extra diffusion term. The first option simply keeps a time step in the system. This is undesireable since we seek a steady-state solution.

By defining a local time scale via the element length and velocity magnitude, the system can be reformed without a time dependence:

The final option involves performing the Taylor series expansion on x rather than t for the steady-state solution. This leads to:

The second method is chosen for the template pseudo-code as detailed below:

{FOMG} = {FOMG}
         +(BETA)(UM11)()(11;-0.5)(B211)(OMGA)
         +(BETA)(UM12)()(12;-0.5)(B211)(OMGA)
         +(BETA)(UM12)()(21;-0.5)(B211)(OMGA)
         +(BETA)(UM22)()(22;-0.5)(B211)(OMGA)
         +(BETA)(UM11)()(13;-0.5)(B212)(OMGA)
         +(BETA)(UM12)()(14;-0.5)(B212)(OMGA)
         +(BETA)(UM12)()(32;-0.5)(B212)(OMGA)
         +(BETA)(UM22)()(24;-0.5)(B212)(OMGA)
         +(BETA)(UM11)()(31;-0.5)(B221)(OMGA)
         +(BETA)(UM12)()(32;-0.5)(B221)(OMGA)
         +(BETA)(UM12)()(41;-0.5)(B221)(OMGA)
         +(BETA)(UM22)()(42;-0.5)(B221)(OMGA)
         +(BETA)(UM11)()(33;-0.5)(B222)(OMGA)
         +(BETA)(UM12)()(34;-0.5)(B222)(OMGA)
         +(BETA)(UM12)()(43;-0.5)(B222)(OMGA)
         +(BETA)(UM22)()(44;-0.5)(B222)(OMGA)

Discussion of Results

GWS


Uniform Mesh

Figures 1-3 show the non-uniform vorticity solution for three Reynolds numbers. It it seen, as expected, that dispersion error becomes a problem as the Reynolds number increases. This is manifested as oscillatory waves around the singularities where the moving top meets the stationary walls.


Figure 1



Figure 2



Figure 3


GWS


Non-Uniform Mesh Adaptation

One way to eliminate the dispersion error is via mesh refinement. Since there are two singularities which exhibit oscillations, it was decided to impliment a symmetrical progression ratio non-uniform mesh as shown in Figure 4. Table 1 details the results of a small optimization study for the symmetric progression ratio. The resulting data does not make much since. The highlighted values select the progression ratio combination that results in the lowest energy norm errors. These values are then used to resolve the system and mitigate for dispersion error as shown in Figures 5-7.


Figure 4


Table 1: Progression ratio optimization.



Figure 5



Figure 6



Figure 7


GWS


Flow Solutions

Here is a representative speed and pressure solution for the case of Re = 1000. The velocity and pressure gradients are simply solved via a post-processing operation.


Figure 8



Figure 9


TWS


Artificial Dispersion

Finally, the TWS is implimented in an attempt to mitigate dispersion error via the introduction of artificial diffusion. A trial and error approach led to an optimal value for beta of 0.2 as shown in Figure 10 below. When implimented alongside the symmetric progression ratio non-uniform meshing, a superior solution is obtained as seen in Figure 11. Upon inspection, it was seen that the energy norms in both the vorticity and streamfunction solutions were lower in all locations for the TWS versus the like GWS solution.


Figure 10



Figure 11

Conclusion

In conclusion, this lab has explored aspects of the streamline-vorticity Navier-Stokes GWS algorithm. The dispersion error that inevitably creeps in at higher Reynolds numbers has been deal with via a combination of non-uniform meshing as well as the addition of artificial diffusion. These combination studies have then led to an optimal TWS solution that, aside from the corner-singularities, accurately expresses the given situation of a steady-state driven cavity.