Chapter 3
Stability Analysis

• the Velocity Correction Scheme (VelocityCorrectionScheme) and
• the Coupled Linearised Navier-Stokes algorithm (CoupledLinearisedNS).

We will consider both cases, highlighting the similarities and differences of these two methods. In this tutorial we will use the Implicitly Restarted Arnoldi Method (IRAM), which is implemented in the open-source library ARPACK and the modified Arnoldi algorithm¹ that is also available in Nektar++.

3.1 Velocity Correction Scheme

First, we will compute the leading eigenvalues and eigenvectors using the velocity correction scheme method. In the $NEKTUTORIAL/stability folder there is a file called Channel-VCS.xml. This file is similar to Channel-Base.xml, but contains additional instructions to perform the direct stability analysis.

Note: The entire GEOMETRY section, and EXPANSIONS section must be identical to that used to compute the base flow.

At the end of the simulation, the terminal screen should look like this:

The eigenvalues are computed in the exponential form Me^iθ where M = |λ| is the magnitude, while θ = arctan(λ_i∕λ_r) is the phase:

(3.1)

It is interesting to consider more general quantities that do not depend on the time length chosen for each iteration T. For this purpose we consider the growth rate σ = ln(M)∕T and the frequency ω = θ∕T.

Figures ?? and ?? show the profile of the computed eigenmode. The eigenmodes associated with the computed eigenvalues are stored in the files Channel_VCS_eig_0.fld and Channel_VCS_eig_1.fld. It is possible to convert this file into VTK format in the same way as previously done for the base flow.

3.2 Coupled Linearised Navier-Stokes algorithm

Note: Remember to use the files provided in the folder Stability/Coupled for this case.

It is possible to perform the same stability analysis using a different method based on the Coupled Linearised Navier-Stokes algorithm. This method requires the solution of the full velocity-pressure system, meaning that the velocity matrix system and the pressure system are coupled, in contrast to the velocity correction scheme/splitting schemes.

Inside the folder $NEKTUTORIAL/stability there is a file called Channel-Coupled.xml that contains all the necessary parameters that should be defined. In this case we will specify the base flow through an analytical expression. Even in this case, the geometry, the type and number of modes are the the same of the previous simulations.

It is important to note that the use of the coupled solver requires that only the velocity component variables are specified, while the pressure is implicitly evaluated.

This has already been set up in the XML file. This is necessary to tell Nektar++ to use the previous solution as the right hand side vector for each Arnoldi iteration.

The terminal screen should look like this:

Using the Stokes algorithm, we are computing the leading eigenvalue of the inverse of the operator L⁻¹. Therefore the eigenvalues of L are the inverse of the computed values² . However, it is interesting to note that these values are different from those calculated with the Velocity Correction Scheme, producing an apparent inconsistency. However, this can be explained considering that the largest eigenvalues associated to the operator L correspond to the ones that are clustered near the origin of the complex plane if we consider the spectrum of L⁻¹. Therefore, eigenvalues with a smaller magnitude may be present but are not associated with the largest-magnitude eigenvalue of operator L. One solution is to consider a large Krylov dimension specified by kdim and the number of eigenvalues to test using nvec. This will however take more iterations. Another alternative is to use shifting but in this case it will make a real problem into a complex one (we shall show an example later). Finally, another alternative is to search for the eigenvalue with a different criterion, for example, the largest imaginary part.

In this case, it is easy to to see that the eigenvalues of the evolution operator L are the same ones computed in the previous section with the time-stepping approach (apart from round-off errors). It is interesting to note that this method converges much quicker that the time-stepping algorithm. However, building the coupled matrix that allows us to solve the problem can take a non-negligible computational time for more complex cases.

3.3 Three-dimensional Channel flow

Now that we have presented the various stability-analysis tools present in Nektar++, we conclude showing the capabilities of the code in three spatial dimensions. In the folder
$NEKTUTORIAL/stability3D there are the files that are required for the stability analysis - note that we do not show the geometry and the base flow generation (we will be using the exact solution) since we have already presented these features in the previous tutorials.

The case considered is similar to the channel flow presented in section ??. However, in this case the Reynolds number is set to 10000. In order to run a three-dimensional simulation, we can either run the full 3D solver by creating a 3D geometry or use a 2D geometry and specify the use of a Fourier expansion in the third direction. The last method is also known as 3D homogenous 1D approach. Here we will present this approach.

Specifically, we use a 2D geometry and we add the various parameters necessary to use the Fourier expansion. Note that in the 2D plane we will use a MODIFIED expansion basis with NUMMODES=11.

Now run the solver - the terminal screen should look like this:

Now convert the two files containing the eigenvectors and visualise them in Paraview or VisIt - the solution should look like the one below:

3.4 Solutions

Completed solutions to the tutorials are available in the completed directory.