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

11.3 Setting up the Simulation

In the following the possible options are shown for setting up a simulation using the incompressible Navier-Stokes.

11.3.1 Expansions and Approximation Spaces

The Expansions are introduced generally in section 3.2. The Expansion section for an incompressible flow simulation can be set as for other solvers regardless of the projection type. Here an example for a 3D simulation with all variables having a same polynomial order (for 2D simulations the specified fields would be just u,v,p). 

1<EXPANSIONS> 
2  <E COMPOSITE="C[0]" NUMMODES="6" FIELDS="u,v,w,p" TYPE="MODIFIED" /> 
3</EXPANSIONS>

In case of a simulation using the Direct Solver we need to set FIELDS=u,v as the pressure expansion order will be automatically set to fulfil the inf-sup condition. Possible choices for the expansion TYPE are:

Basis TYPE
Modal MODIFIED
Nodal GLL_LAGRANGE
Nodal SEM GLL_LAGRANGE_SEM

For well resolved simulations it appears that often using the same polynomial space for the pressure and velocity does give suitable answer but this does not satisfy the so-called LBB or inf-sup condition. Therefore, it is potentially better to specify an equivalent of the Taylor Hood approximation and use one higher polynomial order for velocity than the pressure with a continuous expansion. To specify this type of expansion you can use an expansion section of the form: 

1  <EXPANSIONS> 
2        <E COMPOSITE="C[0]" NUMMODES="8" FIELDS="u,v" TYPE="MODIFIED" /> 
3        <E COMPOSITE="C[0]" NUMMODES="7" FIELDS="p"   TYPE="MODIFIEDQUADPLUS1" /> 
4  </EXPANSIONS>

In the above example the “u,v” fields are specified to have a polynomial order of 7 using a modified expansion. Implicitly this form of the expansion definition uses a quadrature order of 9. The above definition then also uses a modified expansion for pressure but of polynomial order 6. Since currently for this solver to run we need to use a consistent quadrature order for both the velocity and pressure fields we specify the MODIFIEDQUADPLUS1 to tell the solver to use an additional quadrature point and therefore also use 9 quadrature points in each 1D direction for the pressure.

In other cases it is sometimes useful to run with an even higher quadrature order, for example to handle highly deformed elements where the Jacobian is represented by a polynomial expansion. This can be done by using a more detailed definition of the expansion of the form: 

1  <EXPANSIONS> 
2      <E COMPOSITE="C[0]" BASISTYPE="Modified_A,Modified_B" NUMMODES="8,8" POINTSTYPE="GaussLobattoLegendre,GaussRadauMAlpha1Beta0" NUMPOINTS="9,8"  FIELDS="u,v" /> 
3      <E COMPOSITE="C[0]" BASISTYPE="Modified_A,Modified_B" NUMMODES="7,7" POINTSTYPE="GaussLobattoLegendre,GaussRadauMAlpha1Beta0" NUMPOINTS="9,8"  FIELDS="p" /> 
4  </EXPANSIONS>

In this example we have specified an 8th order expansion for “u,v” and a 7th order expansion for “p”. The BasisType is given as “Modified_A, Modified_B” which is for a triangular expansion (note that for a quadrilateral expansion it would have been “Modified_A,Modified_A”) and so the number of quadrature points in this case is 9 in the first direction which uses Gauss-Lobatto-Legendre points but only 8 in the second direction since this uses a Gauss-Radau formula with α = 1,β = 0 weights (see [22] for details).

Further information is also available in Section 3.2.

11.3.2 Governing Equation

The first thing in setting up a simulation is to define what governing eqautions are to be solved. This can be done using the EqType tag in the SOLVERINFO section.

Possible values are:

Equations EQTYPE Dim. Projections Alg.
Steady Stokes (SS) SteadyStokes All CG VCS
Steady Oseen (SO) SteadyOseen All CG DS
Unsteady Stokes (US) UnsteadyStokes All CG VCS
Steady Linearised NS (SLNS) SteadyLinearisedNS All CG DS
Unsteady Linearised NS (ULNS) UnsteadyLinearisedNS All CG VCS,DS
Unsteady NS (UNS) UnsteadyNavierStokes All CG,CG-DG VCS

For example, for solving the unsteady Navier-Stokes equations, the following should be included inside the SOLVERINFO tags 

1<SOLVERINFO> 
2<I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes"/> 
3</SOLVERINFO>

Please note that, we will have only one SOLVERINFO tags in the xml file which inclues all the settings. If more than one SOLVERINFO is defined in the xml session file, the last one will override the previous ones and the settings could be lost or be overwritten with the defaults.

11.3.3 Solution Algorithms

Depend on the type of the equations that will be solved, a consistent sovler scheme should be defined in the SOLVERINFO using SolverType tag.

Available options are as shown next:

Algorithm SolverType Dimensions Projections
Velocity Correction Scheme (VCS) VelocityCorrectionScheme 2D, Quasi-3D, 3D CG, CG-DG
Implicit VCS VCSImplicit 2D and 3D CG
VCS with coordinate transformation VCSMapping 2D, Quasi-3d and 3D CG, CG-DG
VCS with weak pressure VCSWeakPressure 2D, Quasi-3D, 3D CG, CG-DG
Smoothed Profile Method (SPM) SmoothedProfileMethod 2D, Quasi-3D, 3D CG, CG-DG
Direct solver CoupledLinearisedNS 2D, Quasi-3D CG

For exmaple, for the unsteady incompressible Navier-Stokes, using the VCS scheme otherwise known as Velocity Correction Scheme, the follwoing should be added inside the SOLVERINFO tags. 

1<I PROPERTY="SolverType"  VALUE="VelocityCorrectionScheme"/>

or if we want to improve the timestepping stability and hence being able to run with higher CFL condition, this could be done using the implicit VCS as: 

1<I PROPERTY="SolverType"  VALUE="VCSImplicit"/>

Along with setting the solver scheme, we need to set the following tags in SOLVERINFO tags

11.3.4 TimeIntegrationScheme

11.3.5 Stabilisation Techniques

There are various techniques and options available to stablise the simulation as explained in this section. Should any of these options to be used, they are needed to be included in the SOLVERINFO tags.

11.3.5.1 Dealiasing

One source of instablity in the simulation could be related to aliasing. This can be overcome using “Dealiasing” technique. The following options are available and if needed should be included in the SOLVERINFO section of the xml file.

11.3.5.2 Spectral Vanishing Viscousity

Spectral Vanishing Viscousity (SVV) activates a stabilization technique which increases the viscosity on the modes with the highest frequencies. In the Nektar++ there are two kinds of SVV are supported, one acts on the spectral/hp expansions and the other on the Fourier expansions. Additionally, the SVV supports several different kernels as will be shortly explained next.

To activate SVV for the spectral/hp expansions, i.e. in 2D and 3D simulaitons as well as the xy domain in Quasi-3D simulations, the following should be included in the SOLVERINFO tag. 

1<I PROPERTY="SpectralVanishingViscosity" VALUE="True"/>

Using the VALUE="True" will activates the Exponential Kernel by default. There are indeed tree kernels available as summarized next:

SVV Kernel SpectralVanishingViscosity
Exponential Kernel True
Power Kernel PowerKernel
DG Kernel DGKernel

The Exponential kernel is based on the work of Maday et al. [30], its extension to 2D can be found in [23]. A diffusion coefficient can be specified which defines the base magnitude of the viscosity; this parameter is scaled by h∕p. SVV viscosity is activated for expansion modes greater than the product of the cut-off ratio and the expansion order. The Power kernel is a smooth function with no cut-off frequency; it focusses on a narrower band of higher expansion modes as the polynomial order increases. The cut-off ratio parameter for the Power kernel corresponds to the power ratio, see Moura et al. [32].

It is recommended that for the spectral/hp SVV one uses DGKernel as follows: 

1<I PROPERTY="SpectralVanishingViscosity" VALUE="DGKernel"/>

also note that the DGKernel is only supported for the spectral/hp expansions.

Using any of the SVV Kernels, the cut-off ratio and diffusion coefficient will be set by default. However, as mentioned above these are could be modified by the user. For the DGKernel, the default diffusion coefficient is one and can be changed by setting the SVVDiffCoeff parameter in the <PARAMETERS> tag. It is recommended that the defalut values for the DGKernel won’t be changed, however, should the user wanted to changed it for exmaple, to reduce the amount of the diffusion coefficient to “0.75”, this can be done as shown below: 

1<P> SVVDiffCoeff  = 0.75 </P>

For the Exponential Kernel there are two parameters that can be adjusted. The first parameter, i.e. SVVCutOffRatio, controls the cut-off ratio and the second parameter is the SVVDiffCoeff that controls the diffusion coefficient.

A practical hint especially for high-Reynolds number flows is that if the simulation is not very well stable using the default values of exponential kernel, one could start with the cut-off ratio around "0.5" and the diffusion coefficient of "1" and gradually increases the cut-off ratio and/or reducing the diffusion coefficient to get the stable solution with the least amount of dissipation.

In a Quasi-3D simulation, in addition to the spectral/hp modes, one can activate the SVV for the Fourier expansions too. In this case, these two kinds of SVV can be independently activated using the SpectralVanishingViscositySpectralHP for spectral/hp and SpectralVanishingViscosityHomo1D for Forier expansions. Additionally, the cut-off ratio and diffusion coefficients are also can be set independently. For the spectral/hp this would be the same as explained before. For the Fourier expansions the cut-off ratio and diffusion coefficient can be set using SVVCutoffRatioHomo1D and SVVDiffCoeffHomo1D parameters respectively. These paramters need to be set in the PARAMETERS tag.

1<P> SVVDiffCoeff  = 0.1 </P> 
2<P> SVVCutOffRatio= 0.5 </P> 
3<P> SVVDiffCoeffHomo1D  = 0.1 </P> 
4<P> SVVCutOffRatioHomo1D  = 0.5 </P>

11.3.5.3 GJPStabilisation

GJPStabilisation activates the gradient jump penalty[33] stabilization technique. The GJP is less dissipative than the spectral vanishing viscousity (SVV) stabilisation and gives a better turbulent structures specially when the simulaiton is under-resolved or in other words when the polynomial order is low. It is expected that with increasing the polynomial order and approaching towards a resolved solution both SVV and GJP stabilisation solutions approaches each other. Since the GJP is less dissipative than the SVV, it means that it probably requires a smaller Δt for the stable solution. Further, GJP only affects the spectral/hp discretisations and doesn’t have any stabilisation effect on the Fourier discritisation. This means that in the Quasi-3D simulations, if GJP is activated we still needs to use the SpectralVanishingViscousityHomo1D to stablise the solution in the spanwise direction which is discretised using the Fourier method. For the 2D and full 3D problems, GJP stabilisation should suffices though. 

1<I PROPERTY="GJPStabilisation" VALUE="SemiImplicit"/>

Another option for the GJP stabilisation is to use Explicit for its value instead of SemiImplicit. If the Explicit value is set, there is an additional option that can be set for the explicit GJP stabilisation in the solver info as follwos. 

1<I PROPERTY="GJPNormalVelocity" VALUE="True"/>

Using the GJPNormalVelocity, the GJP formulation will uses the average velocity normal to the edge/face as the velocity scaling of jump term.

11.3.6 Simulation Parameters

In Setting up the simulation, various parmeters can be used to set required parameters such as TimeStep or determine the controls such as number of I/O. The following parameters can be specified in the PARAMETERS section of the session file. Note that some of these have already been discussed in previous sections and only repeated here for the sake of completeness.

11.3.7 Boundary conditions

The process of setting up the incompressible solver requires setting apropriate boundary conditions. This includes consistent pressure boundary conditions, outflow boundary conditions for truncated domains. Additionally, for the biomechanical simulaitons such as pulsatile flows, e.g. flow in arteries a specific type of boundary condition such as Wormersley conditions may be required. These boundary conditions are explained in this section.

11.3.7.1 Pressure boundary conditions

In order to specify the pressure boundary conditions given by equation (11.4) or for the equivalent conditions in the VCSWeakPressure scheme the USERDEFINEDTYPE condition “H” can be used. Therefore a zero velocity wall boundary condition on boundary region 0 in two-dimensions can be specified as 

1  <BOUNDARYCONDITIONS> 
2    <REGION REF="0"> 
3      <D VAR="u" VALUE="0" /> 
4      <D VAR="v" VALUE="0" /> 
5      <N VAR="p" USERDEFINEDTYPE="H" VALUE="0"  /> 
6    </REGION> 
7  </BOUNDARYCONDITIONS>

11.3.7.2 Outflow boundary conditions

The most straightforward outflow condition is to specify fully developed conditions of ∇un+1 ⋅n = 0 and p = 0 which can be specified as 

1  <BOUNDARYCONDITIONS> 
2    <REGION REF="0"> 
3      <N VAR="u" VALUE="0" /> 
4      <N VAR="v" VALUE="0" /> 
5      <D VAR="p" VALUE="0"  /> 
6    </REGION> 
7  </BOUNDARYCONDITIONS>

However when energetic vortices pass through an outflow region one can experience instabilities as identified by the work of Dong, Karnidakis and Chryssostomidis [10]. In this paper they suggest to impose a pressure Dirichlet outflow condition of the form

pn+1 = νn ⋅∇u *,n+1 ⋅n - 1-| u*,n+1 |2 S (n ⋅u *,n+1)+ fn+1 ⋅n 2 o b
(11.27)

with a step function defined by So(n ⋅u) = 12(1 - tanh un0⋅uδ), where u0 is the characteristic velocity scale and δ is a non-dimensional positive constant chosen to be sufficiently small. fb is the forcing term in this case the analytical conditions can be given but if these are not known explicitly, it is set to zero, i.e. fb = 0. (see the test KovaFlow_m8_short_HOBC.xml for a non-zero example). Note that in the paper [10] they define this term as the negative of what is shown here so that it could be use used to impose a default pressure values. This does however mean that the forcing term is imposed through the velocity components u,v by specifying the entry VALUE (An example can be found in ChanFlow_m3_VCSWeakPress_ConOBC.xml). For the velocity component one can specify

1[ 1 ] ∇un+1 ⋅ n = --pn+1n + --| u *,n+1 |2 So(n ⋅u*,n+1) - ν(∇ ⋅u*,n+1)n- fnb+1 ν 2
(11.28)

This condition can be enforced using the USERDEFINEDTYPE “HOutflow”, i.e. 

1  <BOUNDARYCONDITIONS> 
2    <REGION REF="0"> 
3      <N VAR="u" USERDEFINEDTYPE="HOutflow" VALUE="0" /> 
4      <N VAR="v" USERDEFINEDTYPE="HOutflow" VALUE="0" /> 
5      <D VAR="p" USERDEFINEDTYPE="HOutflow" VALUE="0" /> 
6    </REGION> 
7  </BOUNDARYCONDITIONS>

Note that in the moving body work of Bao et al. [4] some care must be made to identify when the flow over the boundary is incoming or outgoing and so a modification of the term

1-| u*,n+1 |2 So(n ⋅u*,n+1) 2

is replaced with

1 ( ) -- (θ + α2) | u *,n+1 |2 +(1 - θ + α1 )(u*,n+1 ⋅n )u*,n+1 So(n ⋅u*,n+1) 2

where the default values are given by θ = 1,α1 = 0,α2 = 0 and these values can be set through the parameters OutflowBC_theta, OutflowBC_alpha1 and OutflowBC_alpha2.

Dong has also suggested convective like outflow conditions in [9] which can be enforced through a Robin type specification of the form

n+1 [ ∂u---- + γ0D0-un+1 = 1- fn+1 + E (n,u*,n+1)+ pn+1n - ν(∇ ⋅u *,n+1)n] + D0-^u ∂n Δt ν Δt
(11.29)

n+1 ∂p---- ∂n + --1- νD0pn+1 = -(-ν(∇×∇×u)*,n+1 + N*,n+1)⋅n
--1-- νD0[fn+1 + E(n,u*,n+1 + pn+1n - ν(∇⋅u*,n+1)n] (11.30) 
1<BOUNDARYCONDITIONS> 
2  <REGION REF="0"> 
3   <R VAR="u" USERDEFINEDTYPE="HOutflow" VALUE="0" PRIMCOEFF="D0/TimeStep"/> 
4   <R VAR="v" USERDEFINEDTYPE="HOutflow" VALUE="0" PRIMCOEFF="D0/TimeStep"/> 
5   <R VAR="p" USERDEFINEDTYPE="HOutflow" VALUE="0" PRIMCOEFF="1.0/(D0*Kinvis)"/> 
6  </REGION> 
7</BOUNDARYCONDITIONS>

11.3.7.3 Womersley Boundary Condition

It is possible to define the time-dependent Womersley velocity profile for pulsatile flow in a pipe. The modulation of the velocity profile is based on the desired peak or centerline velocity which can be represented by a Fourier expansion Umax = A(ωn)eiωnt where A are the Fourier modes and ω the frequency. The womersely solution is then defined as:

N 2 ∑ ˜ J0(i3∕2αnr-∕R-) iωnt w (r,t) = A0 (1- (r∕R ) )+ An [1- J0(i3∕2α) ]e n=1

where the womersley number α is defined:

∘ 2πn-- αn = R ---- Tν

and A˜n (n = 1 : N)are the Fourier coefficients scaled in the following way:

˜ ----1---- An = 2An ∕[1- J (i3∕2α )] 0

The Womersley velocity profile is implemented in the following way: 

1<REGION REF="0"> 
2    <D VAR="u" USERDEFINEDTYPE="Womersley:WomParams.xml" VALUE="0" /> 
3    <D VAR="v" USERDEFINEDTYPE="Womersley:WomParams.xml" VALUE="0" /> 
4    <D VAR="w" USERDEFINEDTYPE="Womersley:WomParams.xml" VALUE="0" /> 
5    <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" /> 
6</REGION>

A file containing the Fourier coefficients, Ã, must be in the directory where the solver is called from. The name of the file is defined by the string given in the attribute USERDEFINEDTYPE after the “:” and contains the real and imaginary coefficients. This file has the format 

1<NEKTAR> 
2  <WOMERSLEYBC> 
3    <WOMPARAMS> 
4      <W PROPERTY="Radius" VALUE="0.5" /> 
5      <W PROPERTY="Period" VALUE="1.0" /> 
6      <W PROPERTY="axisnormal" VALUE="0.0,0.0,1.0" /> 
7      <W PROPERTY="axispoint" VALUE="0.0,0.0,0.0" /> 
8    </WOMPARAMS> 
9 
10    <FOURIERCOEFFS> 
11      <F ID="0"> 0.600393641193,    0.0               </F> 
12      <F ID="1"> -0.277707172935,   0.0767582715413   </F> 
13      <F ID="2"> -0.0229953131146,  0.0760936232478   </F> 
14      <F ID="3"> 0.00858135174058,  0.017089888642    </F> 
15      <F ID="4"> 0.0140332527651,   0.0171575122496   </F> 
16      <F ID="5"> 0.0156970122129,   -0.00547357750345 </F> 
17      <F ID="6"> 0.00473626554238,  -0.00498786519876 </F> 
18      <F ID="7"> 0.00204434981523,  -0.00614566561937 </F> 
19      <F ID="8"> -0.000274697215201, 0.000153571881197 </F> 
20      <F ID="9"> -0.000148037910774, 2.68919619581e-05 </F> 
21    </FOURIERCOEFFS> 
22  </WOMERSLEYBC> 
23</NEKTAR>

Each value of à is provided in the FOURIERCOEFFS section and provided as separate entries containing the real and imaginary components, i.e. the mean component provided above is 0.600393641193,0.0.

Similarly in the WOMPARAMS section the key parameters of the boundary condition are also provided as:

11.3.8 Forcing

11.3.8.1 MovingBody

Note: This force type is only supported for the Quasi-3D incompressible Navier-Stokes solver.

This force type allows the user to solve the interaction system of an incompressible fluid flowing past a flexible moving bodies  [37]. By this forcing function, one can eliminate the difficulty of moving mesh by using body-fitted coordinates, so that an additional acceleration term(i.e., forcing term) is introduced to the momentum equations by the non-inertial transform from the deformed and moving coordinate system to non-deformed and stationary one. 

1<FORCE TYPE="MovingBody"> 
2</FORCE>

Available options of the motion type for the moving body include free, constrained and forced vibrations, which can be specified in the TIMEINTEGRATIONSCHEME and SOLVERINFO section. The free type of motion allows the body to move in both streamwise and crossflow directions, while the constrained type limits the motion only in the crossflow direction. For the forced type, the vibration profiles of the body should be specified as a given function or read from input file in MovingBody section. For example: 

1<TIMEINTEGRATIONSCHEME> 
2    <METHOD> IMEX </METHOD> 
3    <ORDER> 2 </ORDER> 
4</TIMEINTEGRATIONSCHEME> 
5 
6<SOLVERINFO> 
7    <I PROPERTY="EQTYPE"                VALUE="UnsteadyNavierStokes"     /> 
8    <I PROPERTY="SolverType"            VALUE="VelocityCorrectionScheme" /> 
9    <I PROPERTY="EvolutionOperator"     VALUE="SkewSymmetric"            /> 
10    <I PROPERTY="Projection"            VALUE="Galerkin"                 /> 
11    <I PROPERTY="GlobalSysSoln"         VALUE="DirectStaticCond"         /> 
12    <I PROPERTY="HOMOGENEOUS"           VALUE="1D"                       /> 
13    <I PROPERTY="USEFFT"                VALUE="FFTW"                     /> 
14    <I PROPERTY="VibrationType"         VALUE="FREE"                     /> 
15</SOLVERINFO>

A moving body type boundary condition should be specified in BOUNDARYCONDITIONS for the velocities on the moving body, 

1<BOUNDARYCONDITIONS> 
2    <REGION REF="0"> 
3        <D VAR="u" USERDEFINEDTYPE="MovingBody" VALUE="0" /> 
4        <D VAR="v" USERDEFINEDTYPE="MovingBody" VALUE="0" /> 
5        <D VAR="w" VALUE="0" /> 
6        <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" /> 
7     </REGION> 
8</BOUNDARYCONDITIONS>

For the simulation of low mass ratio, there is an option to activate fictitious mass method for stabilizing explicit coupling between the fluid solver and structural dynamic solver. Here we need to specify the values of fictitious mass and damping in PARAMETERS, for example, 

1<SOLVERINFO> 
2    <I PROPERTY="FictitiousMassMethod"    VALUE="True" /> 
3</SOLVERINFO> 
4<PARAMETERS> 
5    <P> FictDamp      = 1000    </P> 
6    <P> FictMass       = 1.5    </P> 
7</PARAMETERS>

A filter called MovingBody is encapsulated in this module to evaluate the aerodynamic forces along the moving body surface. The forces for each computational plane are projected along the Cartesian axes and the pressure and viscous contributions are computed in each direction.

The following parameters are supported:

Option name Required Default

Description

OutputFile ✗ session

Prefix of the output filename to which the forces are written.

Frequency ✗ 1

Number of timesteps after which output is written.

Boundary ✓ -

Boundary surfaces on which the forces are to be evaluated.

To enable the filter, add the following to the FORCE tag:: 

1  <FORCE TYPE="MovingBody"> 
2      <PARAM NAME="OutputFile">DragLift</PARAM> 
3      <PARAM NAME="OutputFrequency">10</PARAM> 
4      <PARAM NAME="Boundary"> B[0] </PARAM> 
5  </FORCE>

During the execution a file named DragLift.fce will be created and the value of the aerodynamic forces on boundary 0, defined in the GEOMETRY section, will be output every 10 time steps.evaluates the aerodynamic forces along the moving body surface. The forces for each computational plane are projected along the Cartesian axes and the pressure and viscous contributions are computed in each direction.

Also, to use this module a MAPPING needs to be specified, as described in section 11.7. In the case of free and constrained motion presented here, the functions defined by the mapping act as initial conditions. Also, when using the MovingBody forcing, it is not necessary to set the TIMEDEPENDENT property of the mapping.

11.3.9 Filters

The following filters are supported exclusively for the incompressible Navier-Stokes solver. Further filters from section 3.5 are also available for this solver.

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