[prev] [prev-tail] [tail] [up]

Chapter 3
Linear Stability Analysis

In this chapter, we will compute the leading eigenvalues and eigenvectors for Rayleigh–Bénard Convection (RBC). We have the analytical solution for the base flow, which is simply the solution of the conduction equation for the temperature field in the absence of fluid motion (the same solution used as the initial condition in the previous chapter). Therefore, U = 0 and T = 1 - y. Our goal is to reproduce the value of the critical Rayleigh number, Rac = 1707.762, first predicated by Chandrasekhar1 .

First, we will perform linear stability analysis for Ra = 1900 and Pr = 0.71. In the $NEKTUTORIAL/LSA/Ra_1900 folder there is a file called rbc-LSA.xml. This file is similar to rbc-DNS.xml, but contains additional instructions to perform the linear stability analysis.

Task: 3.1 Configure the following additional SOLVERINFO options which are related to the stability analysis.

  1. set the EvolutionOperator to Direct in order to activate the linearised governing equations (1.8)–(1.9).

  2. set the Driver to ModifiedArnoldi in order to use the ModifiedArnoldi eigenvalue analysis.

Task: 3.2 Set the parameters for the Arnoldi algorithm.

Next, we need to define the InitialConditions. We will use the random initial condition, which is implemented by using following within the CONDITIONS section. 

1<FUNCTION NAME="InitialConditions"> 
2   <E VAR="u" VALUE="awgn(0.0001)" /> 
3   <E VAR="v" VALUE="awgn(0.0001)" /> 
4   <E VAR="T" VALUE="awgn(0.0001)" /> 
5   <E VAR="p" VALUE="0" /> 
6</FUNCTION>

Task: 3.3 Define a function called BaseFlow. Set the following function for the base flow:

U = 0 (3.1)
V = 0 (3.2)
T = 1 - y (3.3)
P = 0 (3.4)

Task: 3.4 Run the solver to perform the analysis

$NEK/IncNavierStokesSolver rbc-LSA.xml

The simulation should converge in 8 iterations, and the terminal screen should look similar to the one below at the end: 

EV: 0 1.13325     0           1.25088     0 
Writing: "rbc-LSA_eig_0.fld" (0.015323s, XML)

Task: 3.5 Plot the leading eigenvector in Paraview or VisIt. This should look like the solution shown in figures 3.1.

PIC

Figure 3.1 The velocity vector field u′ on top of the density plots of T′ for the leading eigenmode.

3.1 Estimate the value of critical Rayleigh number

Task: 3.6 Now again run the linear stability analysis for Ra = 1600,1700,1800. Value of growth rates for these simulation are given in Table 3.1.

Ra σ
  
1600 -0.731021
1700 -0.0519525
1800 0.608177

Table 3.1 Growth rate of the most unstable mode as a function of Rayleigh number.

Task: 3.7 Use linear interpolation to estimate the Rayleigh critical. You should get Rac = 1707.87, which makes a good agreement with the value obtained by Chandrasekhar (1961).
For linear interpolation you can use following Python script: 

In [1]: Ra = np.array([1600, 1700, 1800, 1900]) 
 
In [2]: sigma = np.array([-0.731021, -0.0519525, 0.608177, 1.25088]) 
 
In [3]: np.lnterp(0, sigma, Ra) 
 
Out[3]: 1707.8700467105318

This completes the tutorial.

[prev] [prev-tail] [front] [up]