Global Stability Analysis: Channel Flow

PIC

Global Stability Analysis: Channel Flow

Department of Aeronautics, Imperial College London, UKScientific Computing and Imaging Institute, University of Utah, USA

December 2, 2016

Nektar++ Tutorial

1 Introduction x1.1 Mesh generation 2 Computation of the base flow 3 Stability Analysis x3.1 Velocity Correction Scheme x3.2 Coupled Linearised Navier-Stokes algorithm x3.3 Three-dimensional Channel flow x3.4 Solutions

1 Introduction

[next] [prev] [prev-tail] [tail] [up] Chapterx1Introduction The aim of this tutorial is to introduce the user to the spectral/hp element framework Nektar++ and its use for global stability computations. Information on how to install the libraries, solvers, and utilities on your own computer is available on the webpage www.nektar.info. Task: 1.1 Prepare for the tutorial. Make sure that you have:

Installed and tested Nektar++ v4.3.5from a binary package, or compiled it from source. By default binary packages will install all executables in /usr/bin. If you compile from source they will be in the sub-directory dist/bin of the build directory you created in the Nektar++ source tree. We will refer to the directory containing the executables as $NEK for the remainder of the tutorial.

Downloaded the tutorial files: flow-stability-channel.tar.gz

     Unpack it using unzip flow-stability-channel.tar.gz to produce a directory
     flow-stability-channel with subdirectories called tutorial and complete We

will refer to the tutorial directory as $NEKTUTORIAL.

Task: 1.2 Additionally, you should also install

a visualization package capable of reading VTK files, such as ParaView (which can be downloaded from here) or VisIt (downloaded from here). Alternatively, you can generate Tecplot formatted .dat files for use with Tecplot.

a plotting program capable of reading data from ASCII text files, such as GNUPlot or MATLAB.

Optionally, you can install the open-source mesh generator Gmshx(which can be downloaded from here) to generate the meshes for the tutorial examples yourself. However, pre-generated meshes are provided.

In this tutorial, we will cover the stability analysis of a two-dimensional channel flow, through both a splitting scheme (the Velocity Correction Scheme) and the direct inversion algorithm (also referred to as the Coupled Linearised Navier-Stokes solver). We will briefly show the application of the stability tools presented to a three-dimensional channel flow test case.

Linear stability analysis is a technique that allows us to determine the asymptotic stability of a flow. By decomposing the velocity and pressure in the Navier-Stokes equations as a summation of a base flow (U, P) and perturbation (u′, p′), such that u = U + ϵu′, p = P + ϵp′, with ϵ ≪ 1, we derive the linearised Navier-Stokes equations,

\frac{∂u′}{∂t} + U · ∇u′+ u′· ∇U  = −∇p′ + \frac{1}{Re}∇^{2}u′+ f′,  (1.1)
∇ · u′  = 0.  (1.2)

We will consider a parallel base flow through a 2-D channel (known as Poiseuille flow) at Reynolds number Re = 7500. The velocity has the following analytic form:

. U = (y + 1)(1 − y)e_{x}

(1.3)

The domain is Ω = [−π, π] × [−1, 1] and it is composed by 48 quadrilateral elements as shown in figure 1.1. The problem has been made non-dimensional using the centreline velocity and the channel half-height.

PIC Figurex1.1: 48 quadrilaterals mesh This mesh was created using the software Gmshxand the first step is to convert it into a suitable input format so that it can be processed by the Nektar++ libraries. The files in the $NEKTUTORIAL directory are as follows:

Folder geometry

Channel.geo - Gmshxfile that contains the geometry of the problem Channel.msh - Gmshxgenerated mesh data listing mesh vertices and elements.

Folder base Channel-Base.xml - Nektar++ session file, generated with the $NEK/NekMesh utility, for computing the base flow.

Folder stability/VCS

Channel-VCS.xml - Nektar++ session file, generated with $NEK/NekMesh, for performing the stability analysis. Channel-VCS.rst - Nektar++ field file that contains a set of initial conditions closer to the solution in order to achieve faster convergence.

Folder stability/Coupled

Channel-Coupled.xml - Nektar++ session file, generated with $NEK/NekMesh, for performing the stability analysis.

Folder stability3D

PPF_R10000_3D.xml - Nektar++ session file for performing coupled 3D stability analysis at Re = 10000. PPF_R15000_3D.xml - Nektar++ session file for performing coupled 3D stability analysis at Re = 15000.

1.1 Mesh generation

The first step is to generate a mesh that is readable by Nektar++. The files necessary in this section can be found in $NEKTUTORIAL/geometry/. To achieve this task we use Gmshxin conjunction with the Nektar++ pre-processing utility called $NEK/NekMesh. Specifically, we first generate the mesh in figure 1.1 using Gmshxand successively we convert it into a suitable Nektar++ format using $NEK/NekMesh.

Task: 1.3 Convert the Gmsh geometry provided into the XML Nektar++ format and with two periodic boundaries

Channel.msh can be generated using Gmsh by running the following command:

      gmsh -2 Channel.geo

Channel.xml can be generated using the $NEK/NekMesh pre-processing tool:

      $NEK/NekMesh Channel.msh Channel.xml

Channel-al.xml can be generated using the module peralign available with the pre-processing tool $NEK/NekMesh:

      $NEK/NekMesh -m peralign:surf1=2:surf2=3:dir=x Channel.xml
      Channel-al.xml

where surf1 and surf2 correspond to the periodic physical surface IDs specified in Gmsh (in our case the inflow has a physical ID=2 while the outflow has a physical ID=3) and dir is the periodicity direction (i.e. the direction normal to the inflow and outflow boundaries - in our case x).

Examine the Channel.xml and Channel-al.xml files you have just created. Only the mesh and default expansions are defined at present and the only difference between the two files is the ordering of the edges in the section composite ID=3 which has been re-ordered in order to apply periodic boundary conditions correctly.

Warning: There is currently an issue when using the coupled solver and periodic edges which is being investigated. For achieving the correct channel flow stability results when using the Coupled Linearised Navier-Stokes algorithm (see section 3.2), please use the files provided in

the folder $NEKTUTORIAL/stability/Coupled.

[next] [prev] [prev-tail] [front] [up] 2 Computation of the base flow

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

Chapterx2Computation of the base flow We must first create an appropriate base flow. Since, in hydrodynamic stability theory, it is assumed that the base flow is incompressible, this can be computed using the incompressible Navier-Stokes solver ($NEK/IncNavierStokesSolver). Tip: Note that the incompressible Navier-Stokes solver ($NEK/IncNavierStokesSolver) executable encapsulates the nonlinear equations as well as the stability tools. Therefore, when you setup either a nonlinear incompressible problem or an incompressible stability problem you should use the same executable - i.e.:

 $NEK/IncNavierStokesSolver file.xml.

The instructions for running one or the other are specified on the XML file. For the problem considered here, the specified boundary conditions will be no-slip on the walls and periodic for the inflow/outflow. In this case, since it is not a constant pressure gradient that drives the flow, it is necessary to use a constant body-force in the streamwise direction. It can be shown that this should be equal to 2ν.

In the folder $NEKTUTORIAL/base you will find the file Channel-Base.xml which contains the geometry along with the necessary parameters to solve the problem. The GEOMETRY section defines the mesh of the problem and it is generated automatically as you have seen in the previous task. The expansion type and order is specified in the EXPANSIONS section. An expansion basis is applied to a geometry composite, where by composite we mean a collection of mesh entities (specifically here, a collection of mesh elements), specified in the GEOMETRY section. A default entry is always included by the $NEK/NekMesh utility. In this case the composite C[0] refers to the set of all elements. The FIELDS attribute specifies the fields for which this expansion should be used. The TYPE attribute specifies the kind of the polynomial basis functions to be used in the expansion. For example,

1 < EXPANSIONS> x 2 xxx< ExCOMPOSITE= "C[0]"xNUMMODES= "11"xFIELDS= "u,v,p"xTYPE= "GLL_LAGRANGE "/> x 3 </ EXPANSIONS> .

Note that all the results obtained in this tutorial refers to the expansion parameters just

defined (i.e. NUMMODES="11" FIELDS="u,v,p" TYPE="GLL_LAGRANGE").

In order to complete the problem definition and generate the base flow we need to specify a section called CONDITIONS in the session file. If we examine Channel-Base.xml, we can see how to define the conditions of the particular problem to solve.

In particular, the CONDITIONS section contains the following entries:

Solver information (SOLVERINFO) such as the equation, the projection type (Continuous or Discontinuous Galerkin), the evolution operator (Nonlinear for non-linear Navier-Stokes,

Direct1 1in this case the term Direct refers to the direct stability analysis (opposed to the adjoint analysis) and it has no relation with the Coupled Linearised Navier-Stokes algorithm that will be explained in the next section , Adjoint or TransientGrowth for linearised forms) and the analysis driver to use (Standard, Arpack or ModifiedArnoldi), along with other properties. The solver properties are specified as quoted attributes and have the form

1 < SOLVERINFO> x 2 xxx< IxPROPERTY= "[STRING]"xVALUE= "[STRING]"x/> x 3 xxx...x 4 </ SOLVERINFO>

 Task: 2.1 In the SOLVERINFO section of Channel-Base.xml:

Note: The bits to be completed are identified by xin this file.

set EQTYPE to UnsteadyNavierStokes to select the unsteady incompressible Navier-Stokes equations,

set the EvolutionOperator to Nonlinear in order to select the non-linear Navier-Stokes,

set the Projection property to Continuous in order to select the continuous Galerkin approach,

set the Driver to Standard in order to perform standard time-integration.

The parameters (PARAMETERS) are specified as name-value pairs:

1 < PARAMETERS> x 2 xxx< P> x[KEY]x= x[VALUE]x</ P> x 3 xxx...x 4 </ PARAMETERS>

Parameters may be used within other expressions, such as function definitions, boundary conditions or the definition of other subsequently defined parameters.

Task: 2.2 Declare the two parameters Re, that represents the Reynolds number, and Kinvis, which represents the kinematic viscosity. Now set the Reynolds number to 7500 and the kinematic viscosity to 1/Re - i.e.

 <P> Re = 7500 </P>
<P> Kinvis = 1/Re </P>

Note that you can put previously defined parameters in the VALUE entry which can be an expression.

The declaration of the variable(s) (VARIABLES).

1 < VARIABLES> x 2 xxx< VxID= "0"> xux</ V> x 3 xxx< VxID= "1"> xvx</ V> x 4 xxx< VxID= "2"> xpx</ V> x 5 </ VARIABLES>

The specification of boundary regions (BOUNDARYREGIONS) in terms of composites defined in the GEOMETRY section and the conditions applied on those boundaries (BOUNDARYCONDITIONS). Boundary regions have the form

1 < BOUNDARYREGIONS> x 2 xxx< BxID= "[INDEX]"> x[COMPOSITE− ID]x</ B> x 3 xxx...x 4 </ BOUNDARYREGIONS>

The boundary conditions enforced on a region take the following format and must define the condition for each variable specified in the VARIABLES section to ensure the problem is well-posed.

1 < BOUNDARYCONDITIONS > x 2 xxx< REGIONxREF= "[B− REGION− INDEX]"> x 3 xxxxxx< [TYPE]xVAR= "[VARIABLE_1]"xVALUE= "[EXPRESSION_1]"/> x 4 xxxxxx< [TYPE]xVAR= "[VARIABLE_2]"xVALUE= "[EXPRESSION_2]"/> x 5 xxxxxx...x 6 xxx</ REGION> x 7 xxx...x 8 </ BOUNDARYCONDITIONS >

The REF attribute for a boundary condition region should correspond to the ID="[INDEX]" of the desired boundary region specified in the BOUNDARYREGIONS section.

The definition of the (time- and) space-dependent functions (FUNCTION), in terms of x, y, z and t, such as initial conditions, forcing functions, and exact solutions. The VARIABLES represent the components of the specific function in a specified direction and they must be the same for every function.

1 < FUNCTIONxNAME= "[NAME]"> x 2 xxx< ExVAR= "[VARIABLE_1]"xVALUE= "[EXPRESSION]"/> x 3 xxx< ExVAR= "[VARIABLE_2]"xVALUE= "[EXPRESSION]"/> x 4 xxx...x 5 </ FUNCTION>

Alternatively, one can specify the function using an external Nektar++ field file. For example, this will be used to specify the InitialConditions or ExactConditions.

1 < FUNCTIONxNAME= "[NAME]"> x 2 xxx< FxFILE= "[FILENAME]"/> x 3 </ FUNCTION>

Task: 2.3 Define a function called ExactSolution. For the Poiseuille flow with a streamwise forcing term the exact solution is:

\begin{aligned}U & = (y + 1)(1 − y) & (2.1) \\ V & = 0 & (2.2) \\ P & = 0 & (2.3) \end{aligned}

Note: You need to use the first definition of FUNCTION where you can set an

     EXPRESSION.

Tip: If you are interested in the meaning of the other parameters and options present in the XML file, they should be available in the User-Guide. If not - just ask and we should be able to answer! Task: 2.4 Define a body forcing function in the streamwise direction (called BodyForce):

f_{x} = 2ν = 2*Kinvis.

Note that for using the body force you need the following additional tag outside the section

CONDITIONS:

1 < FORCING> x 2 xxx< FORCExTYPE= "Body"> x 3 xxxxxx< BODYFORCE> xBodyForcex</ BODYFORCE> x 4 xxx</ FORCE> x 5 </ FORCING>

It is possible to specify an arbitrary initial condition. In this case, it was decided to start from the exact solution of the problem in order to have a steady state in just few iterations. If the initial condition is not specified, it will be set to zero by default.

This completes the specification of the problem.

Task: 2.5 Compute the base flow using the Channel-Base.xml session file by typing:

$NEK/IncNavierStokesSolver Channel-Base.xml

At the end of the simulation, the fields will be written to a binary file Channel-Base.fld and the L_{2} error (using the given exact solution) and the L_{∞} error will be printed on the terminal for each of the variables.

In particular, the terminal screen should look like this:

=======================================================================xxxxxxxxxxxxxxxxxEquationType:xUnsteadyNavierStokesxxxxxxxxxxxxxxxxxSessionxName:xChannel-BasexxxxxxxxxxxxxxxxxSpatialxDim.:x2xxxxxxxxxxxMaxxSEMxExp.xOrder:x11xxxxxxxxxxxxxxxExpansionxDim.:x2xxxxxxxxxxxxxxProjectionxType:xContinuousxGalerkinxxxxxxxxxxxxxxxxxxxxAdvection:xexplicitxxxxxxxxxxxxxxxxxxxxDiffusion:xexplicitxxxxxxxxxxxxxxxxxxxxTimexStep:x0.001xxxxxxxxxxxxxxxxxNo.xofxSteps:x1000xxxxxxxxxxCheckpointsx(steps):x500xxxxxxxxxxxxxIntegrationxType:xIMEXOrder3x=======================================================================xInitialxConditions:xxx-xFieldxu:x(y+1)*(1-y)xxx-xFieldxv:x0xxx-xFieldxp:x0xWriting:x"Channel-Base_0.chk"xSteps:x100xxxxxxTime:x0.1xxxxxxxxxxCPUxTime:x1.296sxSteps:x200xxxxxxTime:x0.2xxxxxxxxxxCPUxTime:x0.440151sxSteps:x300xxxxxxTime:x0.3xxxxxxxxxxCPUxTime:x0.440857sxSteps:x400xxxxxxTime:x0.4xxxxxxxxxxCPUxTime:x0.438776sxSteps:x500xxxxxxTime:x0.5xxxxxxxxxxCPUxTime:x0.441416sxWriting:x"Channel-Base_1.chk"xSteps:x600xxxxxxTime:x0.6xxxxxxxxxxCPUxTime:x0.439318sxSteps:x700xxxxxxTime:x0.7xxxxxxxxxxCPUxTime:x0.438448sxSteps:x800xxxxxxTime:x0.8xxxxxxxxxxCPUxTime:x0.443955sxSteps:x900xxxxxxTime:x0.9xxxxxxxxxxCPUxTime:x0.443197sxSteps:x1000xxxxxTime:x1xxxxxxxxxxxxCPUxTime:x0.440219sxWriting:x"Channel-Base_2.chk"xTime-integrationxx:x5.26234sxWriting:x"Channel-Base.fld"x-------------------------------------------xTotalxComputationxTimex=x6sx-------------------------------------------xLx2xerrorx(variablexu)x:x1.7664e-12xLxinfxerrorx(variablexu)x:x3.59475e-12xLx2xerrorx(variablexv)x:x4.79197e-13xLxinfxerrorx(variablexv)x:x1.12599e-11xLx2xerrorx(variablexp)x:x1.68712e-11xLxinfxerrorx(variablexp)x:x5.2737e-12

The final step regarding the base flow is to visualise the flow fields. Specifically, we need to convert convert the .fld file into a format readable by a visualisation post-processing tool. In this tutorial we decided to convert the .fld file into a VTK format and to use the open-source visualisation package called Paraviewx.

Task: 2.6 Convert the file:

$NEK/FieldConvert Channel-Base.xml Channel-Base.fld Channel-Base.vtu

Now open Paraviewxand use File ->Open, to select the VTK file, click the ’Apply’ button to render the geometry, and select each field in turn from the left-most drop-down menu on the toolbar to visualise the output.

Note: You can also open this type of file in VisIt.

In figure 2.1 we show how the base flow just computed should look like.

PIC Figurex2.1: u-component of the velocity Tip: Note that Nektar++ supports also Tecplot. To obtain a Tecplot-readable file you can run the following command:

 $NEK/FieldConvert Channel-Base.xml Channel-Base.fld Channel-Base.dat

[prev] [prev-tail] [front] [up] 3 Stability Analysis

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

Chapterx3Stability Analysis After having computed the base flow it is now possible to calculate the eigenvalues and the eigenmodes of the linearised Navier-Stokes equations. Two different algorithms can be used to solve the equations:

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 ARP ACK and the modified Arnoldi

algorithm1 ^{1International Journal for Numerical Methods in Fluids, 2008; 57:1435-1458} 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.

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

set the EvolutionOperator to Direct in order to activate the forward linearised Navier-Stokes system. set the Driver to Arpack in order to use the ARP ACK eigenvalue analysis.

Instruct ARPACK to converge onto specific eigenvalues through the solver property

ArpackProblemType. In particular, set ArpackProblemType to LargestMag to get

the eigenvalues with the largest magnitude (that determines the stability of the flow).

Note: It is also possible to select the eigenvalue with the largest real part by setting ArpackProblemType to (LargestReal) or with the largest imaginary part

by setting ArpackProblemType to (LargestImag).

Task: 3.2 Set the parameters for the IRAM algorithm.

kdim=16: dimension of Krylov-space,

nvec=2: number of requested eigenvalues,

nits=500: number of maximum allowed iterations,

evtol=1e-6: accepted tolerance on the eigenvalues and it determines the stopping criterion of the method.

Task: 3.3 Configure the two FUNCTION called InitialConditions and BaseFlow.

A  restart  file  is  provided  to

accelerate communications. Set the InitialConditions function to be read from Channel-VCS.rst. The solution will then converge after 16 iterations after it has populated the Krylov subspace.

Note: The restart file is a field file (same format as .fld files) that contains the eigenmode of the system.

Note: Since the simulations often take hundreds of iterations to converge, we will not initialise the IRAM method with a random vector during this tutorial. Normally, a random vector would be used by setting the SolverInfo option

InitialVector to Random.

The base flow file (Channel-Base.fld), computed in the previous section, should

     be copied into the Channel/Stability folder and renamed Channel-VCS.bse.

Now specify a function called BaseFlow which reads this file.

Task: 3.4 Run the solver to perform the analysis

 $NEK/IncNavierStokesSolver Channel-VCS.xml

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

    Iterationx16,xoutput:x0,xido=99xConvergedxinx16xiterationsxConvergedxEigenvalues:x2xxxxxxxxxxMagnitudexxxAnglexxxxxxxGrowthxxxxxxFrequencyxEV:x0x1.00112xxxxx0.124946xxxx0.0022353xxx0.249892xWriting:x"Channel-al_eig_0.fld"xEV:x1x1.00112xxxxx-0.124946xxx0.0022353xxx-0.249892xWriting:x"Channel-al_eig_1.fld"xLx2xerrorx(variablexu)x:x0.0367941xLxinfxerrorx(variablexu)x:x0.0678149xLx2xerrorx(variablexv)x:x0.0276887xLxinfxerrorx(variablexv)x:x0.0649249xLx2xerrorx(variablexp)x:x0.00512347xLxinfxerrorx(variablexp)x:x0.00135455

The eigenvalues are computed in the exponential form Me^{iθ} where M = |λ| is the magnitude, while θ = arctan(λ_{i}/λ_{r}) is the phase:

. λ_{1,2} = 1.00112e^{±0.249892i} .

(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.

PIC PIC Figurex3.1: u′- and v′-component of the eigenmode.

Task: 3.5 Verify that for the channel flow case :

\begin{aligned}σ & = 2.2353 × 10^{−3} \\ ω & = ±2.49892 × 10^{−1} \end{aligned}

and that the eigenmodes match those given in figures 3.1. This values are in accordance with the literature, in fact in Canuto et al., 1988 suggests 2.23497 × 10^{−3} and 2.4989154 × 10^{−1} for growth and frequency, respectively.

Tip: Note that Nektar++ implements also the modified Arnoldi algorithm. You can try to use it for this test case by setting Driver to ModifiedArnoldi. You can now try to re-run the simulation and verify that the modified Arnoldi algorithm provides a results that is consistent with the previous computation obtained with Arpack.

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.

 Task: 3.6 Edit the file Channel-Coupled.xml:

Note: As before the bits to be completed are identified by xin this file.

Set the SolverType property to CoupledLinearisedNS in order to solve the linearised Navier-Stokes equations using Nektar + +’s coupled solver.

     the EQTYPE must be set to SteadyLinearisedNS and the Driver to Arpack.

Set the InitialVector property to Random to initialise the IRAM with a random initial vector. In this case the function InitialConditions will be ignored.

To compute the eigenvalues with the largest magnitude, specify LargestMag in the

     property ArpackProblemType.

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.

Task: 3.7 Continue modifying Channel-Coupled.xml:

It is necessary to set up the base flow. For the SteadyLinearisedNS coupled solver, this is defined through a function called AdvectionVelocity. The u component must be set up to 1 − y^{2}, while the v-component to zero.

For the coupled solver, it is also necessary to define the following additional tag outside of the

CONDITIONS tag:

1 < FORCING> x 2 xxx< FORCExTYPE= "StabilityCoupledLNS"> x 3 xxx</ FORCE> x 4 </ FORCING>

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.

Task: 3.8 Now run the solver to compute the eigenvalues

$NEK/IncNavierStokesSolver Channel-Coupled.xml

The terminal screen should look like this:

=======================================================================xxxxxxxxxxxxxxxxxxSolverxType:xCoupledxLinearisedxNSx=======================================================================xxxxxxxxxArnoldixsolverxtypexxxx:xArpackxxxxxxxxxArpackxproblemxtypexxxx:xLMxxxxxxxxxSinglexFourierxmodexxxx:xfalsexxxxxxxxxBetaxsetxtoxZeroxxxxxxx:xfalsexxxxxxxxxShiftx(Real,Imag)xxxxxx:x0,0xxxxxxxxxKrylov-spacexdimensionx:x64xxxxxxxxxNumberxofxvectorsxxxxxx:x4xxxxxxxxxMaxxiterationsxxxxxxxxx:x500xxxxxxxxxEigenvaluextolerancexxx:x1e-06x======================================================xInitialxConditions:xxx-xFieldxu:x0x(default)xxx-xFieldxv:x0x(default)xMatrixxSetupxCosts:x0.565916xMultilevelxcondensation:x0.098134xxxxxxxxxInitalxvectorxxxxxxx:xrandomxIterationx0,xoutput:x0,xido=-1xWriting:x"Channel-Coupled.fld"xIterationx20,xoutput:x0,xido=1xWriting:x"Channel-Coupled.fld"xIterationx40,xoutput:x0,xido=1xWriting:x"Channel-Coupled.fld"xIterationx60,xoutput:x0,xido=1xWriting:x"Channel-Coupled.fld"xIterationx65,xoutput:x0,xido=99xxConvergedxinx65xiterationsxConvergedxEigenvalues:x4xxxxxxxxxxRealxxxxxxxxImaginaryxEV:xx0xx-0.000328987xxxxxxxxxxxx-0xWriting:x"Channel-Coupled_eig_0.fld"xEV:xx1xxx-0.00131595xxxxxxxxxxxx-0xWriting:x"Channel-Coupled_eig_1.fld"xEV:xx2xxx-0.00296088xxxxxxxxxxxx-0xWriting:x"Channel-Coupled_eig_2.fld"xEV:xx3xxx-0.00526379xxxxxxxxxxxx-0xWriting:x"Channel-Coupled_eig_3.fld"xLx2xerrorx(variablexu)x:x2.58891xLxinfxerrorx(variablexu)x:x1.00401xLx2xerrorx(variablexv)x:x0.00276107xLxinfxerrorx(variablexv)x:x0.0033678

Using the Stokes algorithm, we are computing the leading eigenvalue of the inverse of the operator L^{−1} . Therefore the eigenvalues of L are the inverse of the computed

values2 ^{2L is the evolution operator du/dt = Lu} . 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^{−1} . 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.

Task: 3.9 Set up the Solver Info tag ArpackProblemType to LargestImag and run the simulation again.

=======================================================================xxxxxxxxxxxxxxxxxxSolverxType:xCoupledxLinearisedxNSx=======================================================================xxxxxxxxxArnoldixsolverxtypexxxx:xArpackxxxxxxxxxArpackxproblemxtypexxxx:xLIxxxxxxxxxSinglexFourierxmodexxxx:xfalsexxxxxxxxxBetaxsetxtoxZeroxxxxxxx:xfalsexxxxxxxxxShiftx(Real,Imag)xxxxxx:x0,0xxxxxxxxxKrylov-spacexdimensionx:x64xxxxxxxxxNumberxofxvectorsxxxxxx:x4xxxxxxxxxMaxxiterationsxxxxxxxxx:x500xxxxxxxxxEigenvaluextolerancexxx:x1e-06x======================================================xInitialxConditions:xxx-xFieldxu:x0x(default)xxx-xFieldxv:x0x(default)xMatrixxSetupxCosts:x0.557085xMultilevelxcondensation:x0.101482xxxxxxxxxInitalxvectorxxxxxxx:xrandomxIterationx0,xoutput:x0,xido=-1xWriting:x"Channel-Coupled.fld"xIterationx20,xoutput:x0,xido=1xWriting:x"Channel-Coupled.fld"xIterationx40,xoutput:x0,xido=1xWriting:x"Channel-Coupled.fld"xIterationx60,xoutput:x0,xido=1xWriting:x"Channel-Coupled.fld"xIterationx65,xoutput:x0,xido=99xxConvergedxinx65xiterationsxConvergedxEigenvalues:x4xxxxxxxxxxRealxxxxxxxxImaginaryxEV:xx0xxxx0.00223509xxxxxx0.249891xWriting:x"Channel-Coupled_eig_0.fld"xEV:xx1xxxx0.00223509xxxxx-0.249891xWriting:x"Channel-Coupled_eig_1.fld"xEV:xx2xxxx-0.0542748xxxxxx0.300562xWriting:x"Channel-Coupled_eig_2.fld"xEV:xx3xxxx-0.0542748xxxxx-0.300562xWriting:x"Channel-Coupled_eig_3.fld"xLx2xerrorx(variablexu)x:x2.58891xLxinfxerrorx(variablexu)x:x1.00401xLx2xerrorx(variablexv)x:x0.00276107xLxinfxerrorx(variablexv)x:x0.0033678

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.
 Task: 3.10 In the file $NEKTUTORIAL/stability3D/PPF_R10000_3D.xml, make the following

changes:

Add a SOLVERINFO tag called HOMOGENEOUS and set it to 1D.

Add two additional SOLVERINFO tags called ModeType and BetaZero and set them to SingleMode and True, respectively.

Add two PARAMETERS called HomModesZ and LZ and set them to 2 and 1, respectively.

Add two other PARAMETERS called realShift and imagShift and set them to 0.003 and 0.2, respectively.

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

=======================================================================xxxxxxxxxxxxxxxxxxSolverxType:xCoupledxLinearisedxNSx=======================================================================xxxxxxxxxArnoldixsolverxtypexxxx:xModifiedxArnoldixxxxxxxxxSinglexFourierxmodexxxx:xtruexxxxxxxxxBetaxsetxtoxZeroxxxxxxx:xtruex(overridesxLHom)xxxxxxxxxShiftx(Real,Imag)xxxxxx:x0.003,0.2xxxxxxxxxKrylov-spacexdimensionx:x64xxxxxxxxxNumberxofxvectorsxxxxxx:x2xxxxxxxxxMaxxiterationsxxxxxxxxx:x500xxxxxxxxxEigenvaluextolerancexxx:x1e-06x======================================================xInitialxConditions:xxx-xFieldxu:x0x(default)xxx-xFieldxv:x0x(default)xxx-xFieldxw:x0x(default)xWriting:x"PPF_R10000_3D_0.chk"xMatrixxSetupxCosts:x1.97987xMultilevelxcondensation:x0.427631xxxxxxxxxInitalxvectorxxxxxxx:xrandomxIteration:x0xIteration:x1x(residualx:x4.89954)xIteration:x2x(residualx:x3.64295)xIteration:x3x(residualx:x2.54314)x....xIteration:x20x(residualx:x1.35156e-05)xIteration:x21x(residualx:x1.64786e-06)xIteration:x22x(residualx:x1.92473e-07)xWriting:x"PPF_R10000_3D.fld"xLx2xerrorx(variablexu)x:x3.01846xLxinfxerrorx(variablexu)x:x2.25716xLx2xerrorx(variablexv)x:x1.8469xLxinfxerrorx(variablexv)x:x0.985775xLx2xerrorx(variablexw)x:x5.97653e-06xLxinfxerrorx(variablexw)x:x1.2139e-05xEV:xx0xxxxxx0.518448xxxxxx-26.6405xxxx0.00373022xxxxxx0.162477xWriting:x"PPF_R10000_3D_eig_0.fld"xEV:xx1xxxxxx0.518448xxxxxxx26.6405xxxx0.00373022xxxxxx0.237523xWriting:x"PPF_R10000_3D_eig_1.fld"xWarning:xLevelx0xassertionxviolationxComplexxShiftxapplied.xNeedxtoximplementxRitzxre-evaluationxofxeigenvalue.xOnlyxonexhalfxofxcomplexxvaluexwillxbexcorrect

Now convert the two files containing the eigenvectors and visualise them in Paraviewxor VisIt

- the solution should look like the one below:

PIC PIC Figurex3.2: u′- and v′-component of the eigenmode.

 Task: 3.11 The complete input file $NEKTUTORIAL/stability3D/PPF_R15000_3D.xml has

been provided to show a full 3D unstable eigenmode where β is not zero. Run this file and see that you obtain the eigenvalue 0.00248682 ± −0.158348i

Task: 3.12 You can now see what the difference when not using an imaginary shifting. Set the parameters imagShift=0, kdim=384 and nvec=196.

This should take 195 iterations to complete and hidden in the list of eigenvalues should be the unstable values 0.00248662 ± 0.158347i. They were eigevalues 152 and 153 in my run.

3.4 Solutions

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

This completes the tutorial.

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