UiB : MatNat : Informatikk : Undervisning : Kurs : SC200 : SC200 Thierry Matthey


SC200 - Report 2 - Thierry Matthey, Inga Berre

Rayleigh-Benard instability

[Problem description] [Introduction] [The Model] [Numerical Methods] [Conclusion]


Introduction:

We want to study the stability of a fluid situated between two plane and horizontal plates maintained at two different temperatures, see Fig 1. Each plate is maintained at an uniform temperature. We will try to evaluate the conditions and the nature of the convective movement of the fluid.

Fig 1

In the case that the upper plate has a temperature higher then the temperature of the lower plate, the thermal stratification is stable.

If we consider the opposite situation, where the upper plate is colder then the lower plate, we have heavy fluid above lighter fluid and the fluid is only conditionally stable. If the buoyance force due to Archimedes dominates the drag force due to viscosity, the fluid will be unstable. If the drag force dominates, the fluid will be stable. Thus, if the fluid is inviscid, it will be unstable if T2>T1 no matter how small the temperature difference is.

When the temperature gradient is over a critical value, rolls of convection will appear between the two plates, a phenomenon which is called the Rayleigh-Benard instability. Thus, as the Rayleigh number defined by

where

are the thermal coefficient, temperature difference an kinematic viscosity and k is the thermal diffusivity is increased, there is a transition from heat transfer through conduction through heat transfer through convection. This is a buoyance driven convection and is called thermal convection or Rayleigh-Benard convection. The critical Rayleigh number for which this transition occurs can be estimated from the Nusselt number, defined by
Nu= flow/flow without conduction.
We have that Nu=1 for R < Rc and then increases for R > Rc.

The Model

[Problem description] [Introduction] [The Model] [Numerical Methods] [Conclusion]

Boussinesq approximation

We will use Boussinesq approximation on our study of the R-B instability. This means essentially that the density is assumed to a function only of the temperature and only in the buoyance term, while other fluid parameters such as viscosity and thermal diffusivity are regarded as constants.

Such systems are governed by the Boussinesq equations. It is known that sufficiently close to the onset of convection and in a sufficiently wide cell, the fluid motion is weakly three-dimensional in that only slow spatial and temporal modulations of locally two-dimensional roll pattern can occur. The five three-dimensional Boussinesq equations can then be replaced by a single two-dimensional equation that describes the evolution of a complex amplitude or envelope function whose phase and magnitude slowly modulate the sinusoidal oscillations of the rolls.

2D equation for the 3D convection

A lowest-order nontrivial amplitude equation is obtained by multiple-scale analysis of Boussinesq equations near onset. The perturbation is made in terms of the small parameter

where Rc is the critical Rayleigh number at the Rayleigh-Benard instability for a laterally infinite container. By including the rapid oscillations in a real amplitude field Psi(x,y,tau) where tau =epsilon-1t is the slow time variable, a rotational invariant model equation can be derived. This gives us the following model problem

The model equation can be variationally derived from a Liapunov functional

which decreases monotonically for all initial conditions. F is an "energy" measure for solutions and provides a means of ordering possible roll patterns according to their stability.


Numerical methods

[Finite Difference Scheme]
[Spectral Legendre/ Chebyshev Galerkin method]
[Comparison of Finite Difference and Spectral Method]
[Time dependent equation]

[Problem description] [Introduction] [The Model] [Numerical Methods] [Conclusion]

When solving the model(1) numerically, the biharmonic BVP

with boundary conditions

will have to be solved for every time step. The first part of Exercise 2 will therefore be finding numerical codes for solving this boundary value problem. We are going to compare the 2.nd order Finite Difference Scheme and a Spectral Garlekin Method in solving this problem.


Finite Difference Scheme

Problem 1:

From netlib we down loaded the code bihar for solving our biharmonic boundary value problem.


Problem 2:

To check our program later we chose a particular solution for our problem


Differentiation with respect to x and y gives us

and we get

for the right side of the equation. The boundary conditions becomes

We now have a problem with the known solution u(x,y).


Problem 3:

We made a main program in Fortran77. This program
  • Give values to all arbitrary parameters

  • Set up the boundary conditions (f2 and f3)

  • Call the subroutine bihar

  • Call the subroutine approxit which approximates the solution at an arbitrary point (x,y)

  • Call a subroutine error which calculates the max error, L2 error and relative error

  • Write output data to a file

The plots in this part were done in matlab. The matlab-programs can be found here.

Problem 4:

We solved our problem from 2 numerically with Nx=Ny=49. Figure 2 shows the solution u of the equation. Figure 3 shows the absolute error, e =|uc-ue|. For the our choice for Nx and Ny we got L2=3.9885412672635798E-05. Look in our log-files for the other results. Figure 4 shows how the L2 error of bihar varies for Ndof = (1+2n)(1+2n), n=1,...,1000.

Fig 2: The solution uc



Fig 3: The absolute error



Fig 4: The L2 error vs Ndof




Spectral Legendre/ Chebyshev Galerkin method

Problem 1:

We down loaded the spectral Legendre/Chebyshev Galerkin method Bispec for solving the biharmonic boundary value problem in quadratic domains.

Problem 2:

To check our program later we chose a particular solution for our problem


For the right hand side of the equation we now get

With a = c = -1, b = d = 1 the boundary conditions becomes

We now have a problem with the known solution u(x,y).

Problem 3:

We made a main program in Fortran90. This program
  • Give values to all arbitrary parameters

  • Set up the boundary conditions

  • Call the subroutine leg, which solves the generalized biharmonic equation

  • Call the subroutine funkval which approximates the solution at an arbitrary point (x,y)

  • Call a subroutine error which calculates the max error, L2 error and relative error

  • Write output data to a file

The plots in this part were done in matlab. The matlab-programs can be found here.

Problem 4:

We solved our problem from 2 numerically with Nx=Ny=49 and c=2. Figure 5 shows the solution u of the equation. Figure 6 shows the absolute error, e =|uc-ue|. Look in our log-files for the other results. Figure 7 shows how the L2 error of bispec varies with Ndof.

Fig 5: The solution uc



Fig 6: The absolute error



Fig 7: The L2 error vs Ndof



The plot of the absolute error looked very special in this case. We have a increasing error near the x- and y-axis (but it still not larger than 4*10-14. The reason for this may be that the grid points lay further from each other near the axis. Since we also have that the testfunction is zero on these axis, this could make the error blow up. It would have been interesting to check the absolute error for a finer grid or other test functions to see if the same phenomena occured.


Comparison of Finite Difference and Spectral Method

Look in our log-files for the results from this part of the exercise. The plots were done in matlab. The matlab-programs can be found here.

Problem 1:

Nx and Ny are the number of degrees of freedom in the x- and y- direction, respectively. The total number of degrees of freedom is given by


The computational cost of solving (1) with the Finite Difference Method is O(N2log(N)) and O(N3) with the Spectral Method.

Problem 2:

We chose the exact solution

to problem (1) and let Nx = Ny, a = c = -1, b = d =1. Figure 8 shows the L2 error vs. N using logarithmic scale along the abscissa. The Finite Difference Method use uniformly spaced nodes, and since the sin-function is periodic we get large error for small N. The Spectral Method does not use uniformly spaced nodes, therefore we do not have the same problem using this method. For the Finite Difference method we have a accuracy of order O(dx2) while the Spectral method exponentially reaches the numerical error bound.


Fig 8: The L2 error vs N


Figure 9 shows the CPU-time vs N. Notice the oscillations in the CPU-time for the Finite Difference Method. These are typical for methods using FFT (fast fourier transform). The FFT (O(N2log(N)) can only be used for some values near a number power of 2, otherwise (the slowest case occurs if N is a prime) the method uses a regular solver (O(N3)).


Fig 9: The CPU-time vs N


The two figures shows that even if the Finite Difference Method is faster for large N, the Spectral Method is much more accurate. This suggests that the Spectral Method converges faster even tough it is O(N3).

Problem 3:

Let

be an exact solution to problem (1). We chose Nx = Ny = 200, a = c = -1, b = d = 1. Let c = 1, 5, 10, 15, 20, . . ., 100. Figure 10 shows the L2 error vs. c for the Finite Difference Method (blue) and the Spectral Method (red) using logarithmic scale along the abscissa. Figure 11 shows the CPU-time vs c.


Fig 10: The L2 error vs c


Figure 8 shows the CPU-time vs c.


Fig 11: The CPU-time vs c


Problem 4:

For given c we get (2c)2 'rolls'. Figure 12 and 13 shows the level curves of u(x,y) for c=1 and c=2. We can see that c=1 and c=2 gives us 4 and 16 'rolls' respectively.



Fig 12: The level curves of u(x,y) for c=1. Fig 13: The level curves of u(x,y) for c=2.


By plotting the relative error vs c we can determine how many degrees of freedom we need to resolve a 'roll' with a certain accuracy, using the two methods. If ca is the value for c that gives us the wanted accuracy the degrees of freedom are

Figure 14 shows the relative error vs c for Ndof=2002 .


Fig 14: The relative error vs c


This gives us the following results:

Finite difference Method:


Spectral Method:


NB1:The values for Nf we got by reading ca from the plot are not very accurate, but may indicate the size of the exact Nf.

NB2: If we compute the L2-norm of our exact solution we find that it is constant and equals 3/4. We could therefore have used the result form 3. directly if we had scaled the L2-error with 3/4.


Problem 5:

We are going to compare the CPU time needed to resolve the "roll" patterns corresponding to exact solution br>

to (1) with 1% accuracy and 0.01% accuracy. For given accuracy and c we have to chose Ndof according to the formulae for Nf given in 3. Since the computed N-values are not very accurate, we chose to let the program run for different values of N for the two test function and stop when the wanted accuracy (test on relative error) was reached (noted by Nactual). The results are shown in the table below. (Since it is difficult to get 0.01% accuracy with the Finite Difference Method we use 0.1% accuracy). We did the runs on an Onyx2 4xR10000 (250MHz) CPUs and used the compiler options '-n32 -O3 -mips4' for the MIPSpro compiler.

Finite difference Method:


Spectral Method:


These results shows that the Spectral Method reaches the wanted accuracy much faster than the Finite Difference Method. (As expected form the results in 2). Thus the Spectral Method converges faster even tough it is O(N3).



Time dependent equation

Having discretized the partial differential equation (PDE) in (1) spatially, gives an initial value problem of ordinary differential equations (ODE) with dimension Ndof to be solved. The linear 4.th order term is very stiff and should be integrated implicitly.

Look in our log-files for the results from this part of the exercise. The plots were done in matlab. The matlab-programs can be found here.

Problem 1:

We will discretizise the partial differential equation

Integrating the 4.th order term implicitly using implicit Euler and the nonlinear term using explicit Euler gives us the scheme

or by rearranging the terms

Inserting for

gives us

We now have a problem of the form


to be solved for every time step, where

Problem 2:

We made two programs, biharmonic_movie based the FD code bihar and bispec_movie based on the spectral code bispec, which solves the IBVP (1).

We had to try different values for L and N (degrees of freedom) to get the program to simulate the physical situation. In the physical case you can observe that the "waves" rotates. This happens because of physical "noise" in the system. In the numerical case we do not have this kind of noise, so we are therefore dependent of some numerical noise to get the waves to rotate. We therefore had to try out different values for L and N to get the small instability that made the rotation possible. We found (after some work!) that N=26 and L=16pi gave us the wanted result for the Spectral Method. The system reached the steady state situation after about tau=2700. When using the Finite Difference Method we chose the same region. We set the degrees of freedom to N=101 (since this method is less accurate). We did not chose a larger N since this would make the runs take a long time.

Problem 3:

We varied epsilon from 0.01 to 0.25 and caluclated Psi2 averaged over the domain. This was done after the system had reached the steady state situation. We used L=16pi, d=4pi and initial condition

Figure 15 shows epsilon vs. Psi2 averaged over the domain for the Finite Difference Method (blue) and Spectral Method (red).


Fig 15: Psi2 vs epsilon


For epsilon smaller than 0.03/ 0.04 the process is diffusive and the steady state is at zero-level. For epsilon 0.03/ 0.04 transition from conduction to convection occurs. Beyond that we see only conduction.

Problem 4:

Let epsilon = 0.25, L=16pi, d=4pi. The initial condition is

We calculated Psi2 averaged over the domain for each time step. We also calculated the Liapunov Functional (energy measure)

Figure 16 shows F(tau) vs tau for the Finite Difference Method (blue) and the Spectral Method (red).



Fig 16: F vs tau


The reason why F is not always (but just almost) monotone decreasing is because of the error in computing the integral. If we had chosen a finer grid we would have gotten an even better result. The difference in F(tau) for the two methods after steady state is reached occures because of approximation error (the Finite Diffence method is less accurate than the Spectral Method).

Problem 5:

We calculated Psi2 in 4. at every time step. Figure 17 shows tau vs. Psi2 for the Finite Difference Method (blue) and Spectral Method (red).



Fig 17: Psi2 vs tau


The difference in F(tau) for the two methods after steady state is reached occures because of approximation error (the Finite Diffence method is less accurate than the Spectral Method). Notice that it takes longer time for the waves to rotate when the Spectral Method is used, than with the Finite Difference Method. This is because the Spectral Method is more stable. It therefore takes longer time to create enough numerical noise to get the waves to rotate when this method is used.

The animations were generated by first writing jpeg images (Finite Difference Method and Spectral Method) with matlab. Then we transformed the images to pnm file format and used the Berkeley mpeg_encoder with a parameter file to make a mpeg movie. To run the animation we used mpeg_play.

  • The animation of Psi for the Finite Difference . N = 29. From tau=0 to tau=750 with a step of 50.

  • The animation of Psi for the Spectral Method. N = 26. From tau=0 to tau=3000 with a step of 51.


    		•

    Conclusion

    [Problem description] [Introduction] [The Model] [Numerical Methods] [Conclusion]
    For both methods the L2-error decreased with increasing number of degrees of freedom. When we compared the two methods we observed that the error decreased much faster with increasing N when using the Spectral method. In fact the Spectral Method reaches exponentially the numerical error bound. This suggested that the Spectral Method converged faster than the Finite Difference Method, but since the CPU-time increased faster for the Spectral Method with increasing N we could not be sure which method that converges faster. The next problem was therefore to compare the CPU-time the two methods used for solving a test problem to a given accuracy. The result form this was that the Spectral Method converged much faster (this was as expected since the difference in decrease in error was much greater than the difference in CPU-time for the two methods with increasing N).

    The Finite Difference Method is faster for a given N, but considering the error the Spectral Method computes much faster for a given accuracy. We observed that the Finite Difference Method has problems with N not near numbers power of 2, otherwise it uses FFT to calculate.The computational work for the Finite Difference Method is O(N2log(N)) when FFT is used and O(N3) for the Spectral Method.

    Although the Spectral Method converges faster than the Finite Difference Method, we had trouble making this method simulate the physical situation correct. The method was in fact too accurate and stable! In the physical case you can observe that the "waves" rotates. This happens because of physical "noise" in the system. In the numerical case we do not have this kind of noise, so we are therefore dependent of some numerical noise to get the waves to rotate. We therefore had to try out different values for L and N to get the small instability that made the rotation possible.

    For small epsilon the process is diffusive and the steady state is at zero-level. For epsilon 0.03/ 0.04 (in our runs) transition from conduction to convection occurs. Beyond that we see only convection.

    The Finite Difference Method was quiet easy to use and we had only some minor problems with the indices. For the Spectral Method we spent considerable time to implement our code and we had to do several work-arounds. The non-uninform coordinates and region-transformations also made the implementation more complex.

    We conclude that the Spectral Method gives better results, but the use is not straight forward. For our time dependent problem it was in fact easier to get the Finite Difference Method to simulate the physical problem correctly because instability (physical or numerical) is needed.


    strekOffisiell side

    Last update by matthey@ii.uib.no and ingab@ii.uib.no