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

Chapter 2
Pre-processing

To set up the problem we have three steps. The first is setting up a mesh as discussed in section 2.1. The second one is setting the expansion bases as explained in section 2.2. We also need to configure the problem initial conditions, boundary conditions and parameters which are discussed in 2.3.

2.1 Mesh generation

The first pre-processing step consists in generating the mesh in a Nektar + + compatible format. One option to do this is to use the open-source mesh-generator Gmesh to first create the geometry. The mesh format provided by Gmesh is not consistent with the Nektar + + solvers and, therefore, it needs to be converted. An example of how to do this can be found in the Advection Solver Tutorial.

For two-dimensional simulations, the mesh definition contains 6 tags encapsulated within the GEOMETRY tag. The first tag, VERTEX, contains the spatial coordinates of the vertices of the various elements of the mesh. The second tag, EDGE contains the lines connecting the vertices. The third tag, ELEMENT, defines the elements (note that in this case we have only quadrilateral - e.g. <Q ID="85"> - elements). The fourth tag, CURVED, is used to describe the control points for the curve. Note this tag is only necessary if curved edges or faces are present in the mesh and may otherwise be obmitted. The fifth tag, COMPOSITE, is constituted by the physical regions of the mesh called composite, where the composites formed by elements represent the solution sub-domains - i.e. the mesh sub-domains where we want to solve our set of equations (note that we will use these composites to define expansion bases on each sub-domain in section 2.3) - while the composites formed by edges are the boundaries of the domain where we need to apply suitable boundary conditions (note that we will use these composites to specify the boundary conditions in section 2.3). Finally, the sixth tag, DOMAIN, formally specifies the overall solution domain as the union of the composites forming the solution subdomains (note that the specification of different subdomain - i.e. composites - in this case is not necessary since they are constituted by same element shapes). For additional details on the GEOMETRY tag refer to the User-Guide. 

1<?xml version="1.0" encoding="utf-8" ?> 
2<NEKTAR> 
3   <GEOMETRY DIM="2" SPACE="2"> 
4      <VERTEX> 
5        <V ID="0">-1.00000000e+01 1.00000000e+01 0.00000000e+00 </V> 
6        ... 
7        <V ID="706">-4.93844170e-01 -7.82172325e-02 0.00000000e+00 </V> 
8      </VERTEX> 
9      <EDGE> 
10        <E ID="0"> 0 1 </E> 
11        ... 
12        <E ID="1346"> 706 668 </E> 
13      </EDGE> 
14      <ELEMENT> 
15        <Q ID="0"> 0 1 2 3 </Q> 
16        ... 
17        <Q ID="639"> 1345 1346 1269 615 </Q> 
18     </ELEMENT> 
19     <CURVED> 
20        <E ID="0" EDGEID="1344" NUMPOINTS="4" TYPE="PolyEvenlySpaced"> ... 
21        ... 
22        <E ID="1346" EDGEID="235" NUMPOINTS="4" TYPE="PolyEvenlySpaced"> ... 
23      </CURVED> 
24      <COMPOSITE> 
25        <C ID="100"> E[1268,1271,...,1344,1346] </C> 
26        <C ID="101"> E[3,6,...,1256,1266] </C> 
27        ... 
28        <C ID="0"> Q[0-639] </C> 
29      </COMPOSITE> 
30      <DOMAIN> C[0] </DOMAIN> 
31   </GEOMETRY> 
32</NEKTAR>

Note: In this case the mesh has been defined under the GEOMETRY tag with the EXPANSIONS definition and the CONDITIONS section in the same .xml file. However, the mesh can be a separate input .xml format containing only the geometry definition. Also, note that this mesh is in uncompressed format. In order to reduce the size of a large mesh compressed format should be used.

2.2 Expansion bases

We need to specify the expansion bases we want to use in each of the composites or sub-domains (COMPOSITE="..") introduced in section 2.1: 

1<EXPANSIONS> 
2   <E COMPOSITE="C[0]" NUMMODES="3" FIELDS="rho,rhou,rhov,E" 
3   TYPE="MODIFIED" /> 
4</EXPANSIONS>

For this case there is only one composite, COMPOSITE="C[0]", where NUMMODES is the number of coefficients we want to use for the basis functions (that is commonly equal to P+1 where P is the polynomial order of the basis functions), TYPE allows selecting the basis functions, FIELDS is the solution variable of our problem and COMPOSITE are the mesh regions. For additional details on the EXPANSIONS tag refer to the User-Guide.

Tip: One source of instability is aliasing effects which arise from the nonlinearity of the underlying problem. Dealiasing techniques based on the concept of consistent integration can be applied in order to improve the robustness of the solver. For further information about dealisaing techniques, please check Dealiasing techniques for high-order spectral element methods on regular and irregular grids.

An example of dealiasing technique on for quadrilateral elements: 

1<EXPANSIONS> 
2      <E COMPOSITE="C[0]" BASISTYPE="GLL_Lagrange,GLL_Lagrange" 
3        NUMMODES="5,5" POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre" 
4        NUMPOINTS="10,10" FIELDS="rho,rhou,rhov,E" /> 
5</EXPANSIONS>

2.3 Configuring problem definitions

We will now proceed to set up the various problem parameters, the solver settings, initial and boundary conditions.

Parameters

The case will be run at Mach number equal to M∞ = 0.2, Reynolds number ReL=1 = 200 and Pr = 0.72, with the pressure set to p∞ = 101325 Pa and the density equal to ρ = 1.225 Kg∕m3. The cylinder is defined as an isothermal wall with imposed temperature Twall = 300.15 K.

Within PARAMETERS tag, we can can also define the final physical time of the simulation, FinTime, the number of steps NumSteps, the step-interval when an output file is written IO_CheckSteps and the step-interval when information about the simulation is printed to the screen IO_InfoSteps.

Task: 2.1 In the .xml file under the tag PARAMETERS, define all the flow parameters as described above. These are declared as Mach, Re, Pr, pinf, rhoinf and Twall. Define the number of steps NumSteps as the ratio of the FinalTime to the time-step TimeStep.

Warning: Do not define both Prandtl number and the thermal conductivity parameters. They are correlated and defining both will prevent the simulation to start.

1<PARAMETERS> 
2    <P> TimeStep = 0.00001 </P> 
3    <P> FinTime = 0.01 </P> 
4    <P> NumSteps = FinTime/TimeStep </P> 
5    <P> IO_CheckSteps = 100 </P> 
6    <P> IO_InfoSteps = 100 </P> 
7    <P> GasConstant = 287.058 </P> 
8    <P> Gamma = 1.4 </P> 
9    <P> pInf = 101325 </P> 
10    <P> rhoInf = 1.225 </P> 
11    <P> Mach = 0.2 </P> 
12    <P> cInf = sqrt(Gamma * pInf / rhoInf) </P> 
13    <P> uInf = Mach*cInf </P> 
14    <P> vInf = 0.0 </P> 
15    <P> Twall = 300.15 </P> 
16    <P> Re = 200 </P> 
17    <P> L = 1 </P> 
18    <P> mu = rhoInf * L * uInf / Re </P> 
19    <P> Pr = 0.72 </P> 
20   </PARAMETERS>

Solver Settings

We now declare how the flow will be solved. We want to include the effects of fluid viscosity and heat conduction and consequently the equation type we are going to use is the Navier-Stokes equations.

Note: In Nektar + + the spatial discretization of the compressible Navier-Stokes equations is projected in the polynomial space via discontinuous projection. Specifically we make use of either of the discontinuous Galerkin (DG) method or the Flux-Reconstruction (FR) approach. Consequently, set the Projection to DisContinuous, as Continuous Projection is not supported in the Compressible Flow Solver.

We must specify the advection type which will be the classical DG in weak form. Note Nektar + + also presents the FRDG scheme, which recovers the DG scheme with exact mass matrix, the FRHU scheme, which recovers the DG scheme with lumped mass matrix and the FRSD scheme, which recovers a spectral difference scheme. We must also define the diffusion operator we want to use, which will be local Discontinuous Galerkin and the time integration method which will be the Classical Runge Kutta of order 4.

Tip: When selecting the Advection Type scheme, bear in mind that:

Additionally, we need to define the Upwind Type (i.e. Riemann solver) we want to use for the advection operator. For this problem we will use HLLC (Harten, Lax, van Leer+Contact) Riemann solver. Also, we will use the constant viscosity type.

Note: A Riemann problem is solved at each interface of the computational domain for the advection term. Nektar + + provides ten different Riemann solvers, one exact and nine approximated. The exact one solves the problem using a Newton iterative method. The high accuracy of this method may imply a high computational cost. The approximated Riemann solvers do not take into account the full Riemann problem, these simplifications of the exact solver provide lower computational cost but lower accuracy.

Task: 2.2 In the .xml file under the tag SOLVERINFO, define all the solver parameters as described above. These are declared as EQType, Projection, AdvectionType, DiffusionType, TimeIntegrationMethod, UpwindType, ViscosityType.

1<SOLVERINFO> 
2      <I PROPERTY="EQType" VALUE="NavierStokesCFE" /> 
3      <I PROPERTY="Projection" VALUE="DisContinuous" /> 
4      <I PROPERTY="AdvectionType" VALUE="WeakDG" /> 
5      <I PROPERTY="DiffusionType" VALUE="LDGNS" /> 
6      <I PROPERTY="TimeIntegrationMethod" VALUE="ClassicalRungeKutta4"/> 
7      <I PROPERTY="UpwindType" VALUE="HLLC" /> 
8      <I PROPERTY="ProblemType" VALUE="General" /> 
9      <I PROPERTY="ViscosityType" VALUE="Constant" /> 
10</SOLVERINFO>

Variables

In the VARIABLES tag we set the solution variable. For the 2D case we have: 

1<VARIABLES> 
2    <V ID="0"> rho </V> 
3    <V ID="1"> rhou </V> 
4    <V ID="2"> rhov </V> 
5    <V ID="3"> E </V> 
6</VARIABLES>

Note again the weak enforcement of the boundary conditions. The BCs are applied to the fluxes rather than to the non conservative variables of the problem. For further understanding, please check A guide to the Implementation of the Boundary Conditions.

Boundary Conditions

The BOUNDARYREGIONS tag specifies the regions where to apply the boundary conditions. 

1<BOUNDARYREGIONS> 
2    <B ID="0"> C[100] </B> 
3    <B ID="1"> C[101] </B> 
4    <B ID="2"> C[102] </B> 
5    <B ID="3"> C[103] </B> 
6</BOUNDARYREGIONS>

The next tag is BOUNDARYCONDITIONS by which the boundary conditions are actually specified for each boundary ID specified in the BOUNDARYREGIONS tag. The boundary conditions have been set as explained in section 1.3 

1<!-- Wall --> 
2  <REGION REF="0"> 
3      <D VAR="rho" USERDEFINEDTYPE="WallViscous" VALUE="0" /> 
4      <D VAR="rhou" USERDEFINEDTYPE="WallViscous" VALUE="0" /> 
5      <D VAR="rhov" USERDEFINEDTYPE="WallViscous" VALUE="0" /> 
6      <D VAR="E" USERDEFINEDTYPE="WallViscous" VALUE="0" /> 
7  </REGION> 
8<!-- Farfield --> 
9    <REGION REF="1"> 
10      <D VAR="rho" VALUE="rhoInf" /> 
11      <D VAR="rhou" VALUE="rhoInf*uInf" /> 
12      <D VAR="rhov" VALUE="rhoInf*vInf" /> 
13      <D VAR="E" VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInf*uInf+vInf*vInf)" /> 
14    </REGION> 
15<!-- Inflow --> 
16    <REGION REF="2"> 
17      <D VAR="rho" VALUE="rhoInf" /> 
18      <D VAR="rhou" VALUE="rhoInf*uInf" /> 
19      <D VAR="rhov" VALUE="rhoInf*vInf" /> 
20      <D VAR="E" VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInf*uInf+vInf*vInf)" /> 
21    </REGION> 
22 
23 
24 
25<!-- Outflow --> 
26    <REGION REF="3"> 
27      <D VAR="rho" VALUE="rhoInf" /> 
28      <D VAR="rhou" VALUE="rhoInf*uInf" /> 
29      <D VAR="rhov" VALUE="rhoInf*vInf" /> 
30      <D VAR="E" VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInf*uInf+vInf*vInf)" /> 
31    </REGION>

Note: As explained in section 2.3 Continuous Projection is not supported in the Compressible Flow Solver. Therefore, boundary conditions are specified through Dirichlet BCs and Neumann BCs are not supported.

The initial conditions have been set as explained in section 1.3. 

1<FUNCTION NAME="InitialConditions"> 
2      <E VAR="rho" VALUE="rhoInf"/> 
3      <E VAR="rhou" VALUE="rhoInf*uInf" /> 
4      <E VAR="rhov" VALUE="rhoInf*vInf" /> 
5      <E VAR="E" VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInf*uInf+vInf*vInf)"/> 
6</FUNCTION>

2.4 Artificial Viscosity

In order to stabilise the flow in the presence of flow discontinuities we utilise a shock capturing technique which makes use of artificial viscosity to damp oscillations in the solution, in conjunction with a discontinuity sensor to decide where the addition of artificial viscosity is needed.

Tip: In order to turn the NonSmooth artificial viscosity on:

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