Lab 8 - Unsteady Convection and Diffusion, GWS, TWS and FVS


For a 1-D linear unsteady convection-diffusion model problem with 40 elements, this lab attends to the following objectives:

Problem Statement

The problem at hand begins with a mathematical description of the classic convection-diffusion problem.

The well-proven Taylor series modification produces four new parameters defined as alpha, beta, gamma, and mu.

A rearrangement of the above equation shows how the four new parameters affect the time term and also supply extra diffusion.

To move towards a computational form, we define the Galerkin Weak Statement...

... and finally, to handle the time dependency, the GWSh is modified by another Taylor Series and the trapezoidal rule which in turn introduces the Euler implicitness parameter (theta).

It should be clearly evident the left-hand side of this equation represents the Jacobian matrix while the right-hand side represents the residual matrix. We thus naturally choose the Newton method to solve the equation.

The aPSE template pseudo-code for the [JAC] and {RES} matrices is seen below. The 'mu' TWS modification is thrown out for the more useful other three parameters. Finally, since the problem at hand is driven solely by convection, the diffusion is turned-off by simply assigning DIFF a value of zero.

{RES} = ()()(U1)(0;0)(A3001)(Q)
[JAC] = ()()()(0;1)(A200)()


Discussion of Results

GWS + ØTS Implimentation

Gaussian IC,    Ø = 1/2,    0.25 < C < 1.0

Since the velocity and x-step are both unity, the Courant number is simply a measure of the time stepping for the system. As seen in Figures 1-4 below, the ability of the system to propertly propagate the Gaussian IC dimishes with increasing C. This is logical since an increase in C is nothing more than an increase in the time step. Since theta is 1/2, all the time step dissipation comes from the choice of time step.

Figure 1: Ø = 0.5,    C = 0.25

Figure 2: Ø = 0.5,    C = 0.5

Figure 3: Ø = 0.5,    C = 0.75

Figure 4: Ø = 0.5,    C = 1.0

In Figure 5 below, theta is increased to 3/4 resulting in added time step dissipation. For C = 0.5, the increase in theta is seen to smooth out any trailing oscillations but at the cost of severe amplitude degredation.

Figure 5: Ø = 0.75,    C = 0.50

TWS + ØTS Implimentation

Algorithm Comparisons with Gaussian IC

Table 1 outlines the algorithms chosen for the TWS comparisons. It should be noted that Case A is simply the optimization of the Jiang Least Squares following the theory as detailed in the course lecture (Viewgraph CST.14).

Table 1: TWS Algorithm Identification.

As shown below, the Jiang Least Squares optimization actually creates an unreasonable amount of dispersion. From an 'eyeball' perspective, the Raymond/Garder and Jiang Least Squares algorithms seem to perform equally well. Case D and Case E turn out to be equivalent, resulting in perfect convection of the Gaussian. This required an increase of the time step from 0.5 to 1, however. Figure 10 impliments the Case D-E algorithms with C = 0.5 which provides a direct comparison with the Raymond/Garder and Jiang Least Squares algorithms. As seen, it no longer produces exact results, but does beat its competitors.

Figure 6: CASE A - Modified Jiang Least Squares

Figure 7: CASE B - Raymond/Garder

Figure 8: CASE C - Jiang Least Squares

Figure 9: CASE D & E - Theory-derived Optimal

Figure 10: Modified Case D & E (C = 0.5)

Figure 12 provides a direct comparison of all the algorithms for C = 1/2 at the final time step. The actual curve was theorized by taking the Gaussian IC and shifting it in the domain by 30 nodes (60 iterations * C for U = 1 and dX = 1 yields a displacement of 30 units). It is interesting to note how every algorithm actual retards the actual displacement. This perfectly jives with the theory which we developed in class that showed a lag in the phase velocity due to O(dx6) truncation error!

Figure 12: Comparison of above cases after 60 steps.

Finally, the above algorithms are compared in the energy norm for C = 1/2 at the final time step. It should be noted that we can no longer use the energy norm to judge algorithms. As seen in the figure, the energy norm comparison indicates superiority in the Case A solution. Of course this would have a lower energy norm and a more even distribution since it experiences extreme diffusion.

Figure 13: Comparison of energy norms for above cases at 60th step.

TWS + ØTS Implimentation

Sine Wave IC

Now we attempt to propagate a sine wave. The results follow the above conclusions to no surprise. Under the best C = 1/2 algorithm, we see trailing dispersion errors with a consistent magnitude. Only under the optimal situation is the sine wave perfectly propagated.

Figure 14: Modified Case D & E (C = 0.5)

Figure 15: Optimal Case D & E

TWS + ØTS Implimentation

Square Wave IC

For the square wave, the best C=1/2 algorithm produces truly horrible results. This doesn't come as a big surprise if I remember that the fourier reconstruction of a square wave involves an infinite series of periodic components while the sine wave only has one periodic component. I would theorize that a good majority of those higher order frequencies are being truncated as the wave passes through successive TWS modifications. Only the truly optimal algorithm preserves the square wave as it propagates.

Figure 16: Modified Case D & E (C = 0.5)

Figure 17: Optimal Case D & E


In conclusion, this lab clearly summarizes the theory learned in our Computational Spectral Theory lectures. It is clearly seen that the TWS terms can provide time dispersion and artificial diffusion to counter the 2dX dispersive error mechanism in the GWS. The Fourier modal analysis spectral theory learned in class has guided us to the optimal algorithm a priori as shown in the unsteady convection solutions of this lab. Thus, the complicated mathematics I endured in the theory, as demonstrated by the 20-page homework assignment, actually did pay off - as this lab verifies!