Examples

9.4 Examples

9.4.1 Shock capturing

Compressible flows can be characterised by abrupt changes in flow variables within the flow domain often referred to as shocks. These discontinuities can lead to numerical instabilities (Gibbs phenomena). This problem is prevented by locally adding a diffusion term to the equations to damp the numerical oscillations.

9.4.1.1 Non-smooth artificial viscosity model

For the non-smooth artificial viscosity model the added artificial viscosity is constant in each element and discontinuous between the elements. The Euler system is augmented by an added Laplacian term on right hand side of equation 9.1 [39]. The diffusivity of the system is controlled by a variable viscosity coefficient ε. For consistency ε is proportional to the element size and inversely proportional to the polynomial order. Finally, from physical considerations ε needs to be proportional to the maximum characteristic speed of the problem. The final form of the artificial viscosity is

ε = ε hλ S, 0p max

(9.8)

where S is a sensor.

To enable the non-smooth viscosity model, the following line has to be added to the SOLVERINFO section:

1<SOLVERINFO> 
2  <I PROPERTY="ShockCaptureType"      VALUE="NonSmooth" /> 
3<SOLVERINFO>

9.4.1.2 Modal sensor

As shock sensor, a modal resolution-based indicator is used

( ) s = log ⟨q --q˜,q---˜q⟩ , e 10 ⟨q,q⟩

(9.9)

where ⟨⋅,⋅⟩ represents a L² inner product, q and are the full and truncated expansions of a state variable (in our case density)

N(P) N (P -1) ∑ ∑ q(x) = ^qiϕi, ˜q(x) = ^qiϕi, i=1 i=1

(9.10)

then the constant element-wise sensor is computed as follows

(| 0 if s < s - κ { 1( π(se-s0)) e 0 Sε = |( 2 1 + sin 2κ if |se - s0| ≤ κ , 1 if se > s0 + κ

(9.11)

where s₀ = s_κ - 4.25 log₁₀(p).

The modal sensor is enabled by default and it can be explicilty set adding the following line to the SOLVERINFO section:

1<SOLVERINFO> 
2  <I PROPERTY="ShockSensorType"       VALUE="Modal" /> 
3<SOLVERINFO>

The diffusivity and the sensor can be controlled by the following parameters:

1<PARAMETERS> 
2<P> Skappa              = -1.3                     </P> 
3<P> Kappa               = 0.2                      </P> 
4<P> mu0                 = 1.0                      </P> 
5</PARAMETERS>

pict pict

Figure 9.1: (a) Steady state solution for M = 0.8 flow at α = 1.25^∘ past a NACA 0012 profile, (b) Artificial viscosity (ε) distribution

9.4.2 Variable polynomial order

A sensor based p-adaptive algorithm is implemented to optimise the computational cost and accuracy. The DG scheme allows one to use different polynomial orders since the coupling between different elements are determined by common numerical fluxes and there is no further coupling between the elements. Furthermore, the initial p-adaptive algorithm uses the same sensor as the shock capturing algorithm to identify the smoothness of the local solution so it rather straightforward to implement both algorithms at the same time.

The polynomial order in each element can be adjusted based on the sensor value that is obtained. Initially, a converged solution is obtained after which the sensor in each element is calculated. Based on the determined sensor value and the pre-defined sensor thresholds, it is decided to increase, decrease or maintain the degree of the polynomial approximation in each element and a new converged solution is obtained.

( || pe - 1 if se > sds |{ pe + 1 if ssm < se < sds pe = | p if s < s < s ||( e fl e sm pe - 1 if se < sfl

(9.12)

For now, the threshold values s_e, s_ds, s_sm and s_fl are determined empirically by looking at the sensor distribution in the domain. Once these values are set, two .txt files are outputted, one that has the composites called VariablePComposites.txt and one with the expansions called VariablePExpansions.txt. These values have to copied into a new .xml file to create the adapted mesh.

9.4.3 De-Aliasing Techniques

Aliasing effects, arising as a consequence of the nonlinearity of the underlying problem, need to be address to stabilise the simulations. Aliasing appears when nonlinear quantities are calculated at an insufficient number of quadrature points. We can identify two types of nonlinearities:

PDE nonlinearities, related to the nonlinear and quasi-linear fluxes.
Geometrical nonlinearities, related to the deformed/curves meshes.

We consider two de-aliasing strategies based on the concept of consistent integration:

Local dealiasing: It only targets the PDE-aliasing sources, applying a consistent integration of them locally.
Global dealiasing: It targets both the PDE and the geometrical-aliasing sources. It requires a richer quadrature order to consistently integrate the nonlinear fluxes, the geometric factors, the mass matrix and the boundary term.

Since Nektar++ tackles separately the PDE and geometric aliasing during the projection and solution of the equations, to consistently integrate all the nonlinearities in the compressible NavierStokes equations, the quadrature points should be selected based on the maximum order of the nonlinearities:

Qmin = Pexp + max-(2Pexp,Pgeom) + 3- 2 2

(9.13)

where Q_min is the minimum required number of quadrature points to exactly integrate the highest-degree of nonlinearity, P_exp being the order of the polynomial expansion and P_geom being the geometric order of the mesh. Bear in mind that we are using a discontinuous discretisation, meaning that aliasing effect are not fully controlled, since the boundary terms introduce non-polynomial functions into the problem.

To enable the global de-aliasing technique, modify the number of quadrature points by:

1<E COMPOSITE="[101]" 
2   BASISTYPE="Modified_A,Modified_A" 
3   NUMMODES="7,7" 
4   POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre" 
5   NUMPOINTS="14,14" 
6   FIELDS="rho,rhou,rhov,E" 
7/>

where NUMMODES corresponds to P+1, where P is the order of the polynomial used to approximate the solution. NUMPOINTS specifies the number of quadrature points.

9.4.4 Implicit solver

In this example, we solve a compressible flow past a circular cylinder using an implicit discontinuous Galerkin compressible flow solver as shown in figure 9.2. For the implicit time-integration schemes, TimeStep or CFL can be adopted to control the time step. For the case of using CFL, the CFL number can grow from CFL to CFLEnd by a ratio of CFLGrowth to adjust the time step at different stages of the simulation.

The CFL number, growing CFL number and the maximum value of the CFL number are controlled by the following parameters as also described in section 9.3

1<PARAMETERS> 
2  <P> CFL            = 0.1                  </P> 
3  <P> CFLGrowth      = 1.1                  </P> 
4  <P> CFLEnd         = 2.0                  </P> 
5</PARAMETERS>

In this case, the numerical simulation starts from a CFL number of 0.1 and grows by a ratio of 1.1 up to a maximum value of CFL number of 2.0. Note that the CFLEnd parameter may assume higher values, depending on the strategy adopted. In addition, there is no need to define the TimeStep parameter, since the time step size is calculated based on the CFL number in each time step.

Since we are solving an implicit time-integration scheme, we must specify to the solver information the EQType which corresponds to an implicit solver. It should be noted that currently the CFS solver only supports NavierStokesImplicitCFE and EulerImplicitCFE.

1<PARAMETERS> 
2  <I PROPERTY="EQTYPE"                        VALUE="NavierStokesImplicitCFE" /> 
3  <I PROPERTY="TimeIntegrationMethod"         VALUE="DIRKOrder2"              /> 
4</PARAMETERS>

There are some other parameters controlling the performance of the implicit solver. NonlinIterTolRelativeL2 determines the convergence tolerance of the nonlinear system solver relative to the initial nonlinear system residual. NekNonlinSysMaxIterations is the maximum iteration number of the nonlinear system solver. LinSysRelativeTolInNonlin determines the convergence tolerance of linear system solver in each nonlinear iteration. NekLinSysMaxIterations is the maximum iteration number of the linear system solver in each nonlinear iteration. LinSysMaxStorage determines the maximum number of variable vector allowed to store in the linear system solver. Specifically for GMRES solver, the GMRES solver will be restarted if NekLinSysMaxIterations is larger than LinSysMaxStorage and thus LinSysMaxStorage determines the storage consumption of the GMRES solver. Regarding the parameters to control the preconditioners in GMRES, PreconMatFreezNumb specifies the number of time steps to freeze the preconditioning matrices, in other words, the preconditioning matrices will be updated based on the number of time steps provided by the user. PreconItsStep determines the number of preconditioning iterations to calculate the preconditioned vector. The default parameters are listed below.

1<PARAMETERS> 
2  <P> NonlinIterTolRelativeL2   = 1E-3    </P> 
3  <P> NekNonlinSysMaxIterations = 10      </P> 
4  <P> LinSysRelativeTolInNonlin = 5.0E-2  </P> 
5  <P> NekLinSysMaxIterations    = 30      </P> 
6  <P> LinSysMaxStorage          = 30      </P> 
7  <P> PreconMatFreezNumb        = 200     </P> 
8  <P> PreconItsStep             = 7       </P> 
9</PARAMETERS>

Here, we choose to solve the compressible Navier-Stokes equations and use the 2nd order Singly Diagonally Implicit Runge–Kutta (SDIRK) method.

Figure 9.2: Laminar flow past a circular cylinder at Re = 200 and M = 0.2.

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