compressible-flow.tex 25.3 KB
Newer Older
1
\chapter{Compressible Flow Solver}
2

3
\section{Synopsis}
4
The CompressibleFlowSolver allows us to solve
5 6 7
the unsteady compressible Euler and Navier-Stokes
equations for 1D/2D/3D problems using a discontinuous
representation of the variables. In the following we describe
8 9
both the compressible Euler and the Navier-Stokes equations.

10
\subsection{Euler equations}
11 12
The Euler equations can be expressed as a hyperbolic
conservation law in the form
13
\begin{equation}\label{eq:euler}
14 15
\frac{\partial \mathbf{q} }{\partial t} + \frac{\partial \mathbf{f}_i}{\partial x}
+ \frac{\partial \mathbf{g}_i}{\partial y} +
16 17
\frac{\partial \mathbf{h}_i}{\partial z} = 0,
\end{equation}
18 19 20 21
where $\mathbf{q} $ is the vector of the conserved variables,
$\mathbf{f}_i =  \mathbf{f}_i (\mathbf{q})$, $\mathbf{g}_i
= \mathbf{g}_i (\mathbf{q})$ and $\mathbf{h}_i =
\mathbf{h}_i (\mathbf{q})$ are  the vectors of the
22 23 24 25 26 27
inviscid fluxes
\begin{equation}
\mathbf{q} =\left\{
\begin{array}{c}
\rho \\
\rho u \\
28
\rho v \\
29 30 31 32 33 34 35
\rho w \\
E
\end{array} \right\}, \hspace{0.5 cm}
\mathbf{f}_i =\left\{
\begin{array}{c}
\rho u \\
p + \rho u^2 \\
36
\rho u v \\
37 38 39 40 41 42 43
\rho u w \\
u (E + p)
\end{array} \right\}, \hspace{0.5 cm}
\mathbf{g}_i =\left\{
\begin{array}{c}
\rho v \\
\rho u v \\
44
p + \rho v^2 \\
45 46 47 48 49 50 51
\rho v w \\
v (E + p)
\end{array} \right\}, \hspace{0.5 cm}
\mathbf{h}_i =\left\{
\begin{array}{c}
\rho w \\
\rho u w \\
52
\rho v w \\
53 54 55 56 57
p + \rho w^2 \\
w (E + p)
\end{array} \right\},
\label{Euler}
\end{equation}
58 59
where $\rho$ is the density, $u$, $v$ and $w$ are the velocity components in $x$, $y$ and
$z$ directions, $p$ is the pressure and $E$ is the total energy. In this work we considered
60 61
a perfect gas law for which the pressure is related to the total energy by the following expression
\begin{equation}
62
E = \frac{p}{\gamma - 1} + \frac{1}{2} \rho(u^2 + v^2 + w^2),
63 64 65
\label{energy}
\end{equation}
where $\gamma$ is the ratio of specific heats.
66
\subsection{Compressible Navier-Stokes equations}
67 68 69 70
The Navier-Stokes equations include the effects of fluid viscosity and heat
conduction and are consequently composed by an inviscid and a viscous
flux. They depend not only on the conserved variables but also, indirectly,
on their gradient. The second order partial differential equations for the
71 72
three-dimensional case can be written as:
\begin{equation}
73 74 75
\frac{\partial \mathbf{q} }{\partial t} +
\frac{\partial \mathbf{f}}{\partial x} +
\frac{\partial \mathbf{g}}{\partial y} +
76 77
\frac{\partial \mathbf{h}}{\partial z} = 0,
\end{equation}
78 79 80 81
where $\mathbf{q} $ is the vector of the conserved variables, $\mathbf{f} =
\mathbf{f} (\mathbf{q}, \nabla (\mathbf{q}))$,
$\mathbf{g}=  \mathbf{g} (\mathbf{q}, \nabla (\mathbf{q}))$ and $\mathbf{h} =
\mathbf{h} (\mathbf{q}, \nabla (\mathbf{q}))$
82 83 84 85 86
are the vectors of the  fluxes which can also be written as:
\begin{equation}
\begin{array}{l}
\mathbf{f} = \mathbf{f}_i - \mathbf{f}_v,
\mathbf{g} = \mathbf{g}_i - \mathbf{g}_v,
87
\mathbf{h} = \mathbf{h}_i - \mathbf{h}_v,
88 89
\end{array}
\end{equation}
90 91
where $\mathbf{f}_i$, $\mathbf{g}_i $ and $\mathbf{h}_i$ are the inviscid fluxes
of Eq. (\ref{Euler}) and $\mathbf{f}_v$, $\mathbf{g}_v $ and $\mathbf{h}_v$ are
92 93 94 95 96 97 98
the viscous fluxes which take the following form:
\begin{equation}
\begin{array}{c} \vspace{0.2 cm}
\mathbf{f}_v =\left\{
\begin{array}{c}
0 \\
\tau_{xx} \\
99
\tau_{yx} \\
100 101 102 103 104 105 106
\tau_{zx} \\
u \tau_{xx} + v \tau_{yx} + w \tau_{zx} + k T_x
\end{array} \right\},  \hspace{0.2 cm}
\mathbf{g}_v =\left\{
\begin{array}{c}
0 \\
\tau_{xy} \\
107
\tau_{yy} \\
108 109 110 111 112 113 114
\tau_{zy} \\
u \tau_{xy} + v \tau_{yy} + w \tau_{zy} + k T_y
\end{array} \right\}, \\ \vspace{0.2 cm}
\mathbf{h}_v =\left\{
\begin{array}{c}
0 \\
\tau_{xz} \\
115
\tau_{yz} \\
116 117 118
\tau_{zz} \\
u \tau_{xz} + v \tau_{yz} + w \tau_{zz} + k T_z
\end{array} \right\},
119
\end{array}
120 121
\label{viscous flux}
\end{equation}
122 123 124 125
where $\tau_{xx}$, $\tau_{xy}$, $\tau_{xz}$, $\tau_{yx}$, $\tau_{yx}$,
$\tau_{yy}$, $\tau_{yz}$, $\tau_{zx}$, $\tau_{zy}$ and $\tau_{zz}$ are
the components of the stress tensor\footnote{Note that we use Stokes
hypothesis $\lambda = -2/3$.}
126 127
\begin{equation}
\begin{array}{cc} \vspace{0.2 cm}
128 129 130 131 132
\tau_{xx} = 2 \mu \left( u_x - \frac{u_x + v_y + w_z}{3}  \right), &  \tau_{yy} =
2 \mu \left( v_y - \frac{u_x + v_y + w_z}{3}  \right), \\
\vspace{0.2 cm}\tau_{zz} = 2 \mu \left( w_z - \frac{u_x + v_y + w_z}{3}  \right),
& \tau_{xy} = \tau_{yx} = \mu(v_x + u_y), \\
\vspace{0.2 cm} \tau_{yz} = \tau_{zy} = \mu(w_y + v_z),  & \tau_{zx} = \tau_{xz} =
133 134 135
\mu(u_z + w_x). \\
\end{array}
\end{equation}
136
where $\mu$ is the dynamic viscosity calculated using
137 138
the Sutherland's law and $k$ is the thermal conductivity.

139
\subsection{Numerical discretisation}
140
In Nektar++ the spatial discretisation of the Euler and of the Navier-Stokes
141
equations is projected in the polynomial space via a discontinuous projection.
142 143 144 145 146 147 148 149 150 151 152
Specifically we make use either of the discontinuous Galerkin (DG) method
or the Flux Reconstruction (FR) approach.
In both the approaches the physical domain $\Omega$ is divided into a mesh
of $N$ non-overlapping elements $\Omega_{e}$ and the solution is allowed
to be discontinuous at the boundary between two adjacent elements.
Since the Euler as well as the Navier-Stokes equations are defined locally
(on each element of the computational domain), it is necessary to define
a term to couple the elements of the spatial discretisation in order to allow
information to propagate across the domain. This term, called numerical
interface flux, naturally arises from the discontinuous Galerkin formulation
as well as from the Flux Reconstruction approach.
153

154
For the advection term it is common to solve a Riemann problem at each
155
interface of the computational domain through exact or approximated Riemann
156 157 158 159 160 161 162 163 164 165 166 167
solvers. In Nektar++ there are different Riemann solvers, one exact and nine
approximated. The exact Riemann solver applies an iterative procedure to
satisfy conservation of mass, momentum and energy and the equation of
state. The left and right states are connected either with the unknown variables
through the Rankine-Hugoniot relations, in the case of shock, or the isentropic
characteristic equations, in the case of rarefaction waves. Across the contact
surface, conditions of continuity of pressure and velocity are employed.
Using these equations the system can be reduced to a non-linear algebraic
equation in one unknown (the velocity in the intermediate state) that is solved
iteratively using a Newton method. Since the exact Riemann solver gives
a solution with an order of accuracy that is related to the residual in the
Newton method, the accuracy of the method may come at high computational
168 169
cost. The approximated Riemann solvers are simplifications of the exact solver.

170 171
Concerning the diffusion term, the coupling between the elements is achieved
by using a local discontinuous Galerkin (LDG) approach as well as five different
172 173 174
FR diffusion terms.

The boundary conditions are also implemented by exploiting the numerical
175 176
interface fluxes just mentioned.
For a more detailed description of the above the interested reader can refer
177 178
to \cite{DeGMen14} and \cite{MenDeG14}.

179
\section{Usage}
180 181 182
\begin{lstlisting}[style=BashInputStyle]
CompressibleFlowSolver session.xml
\end{lstlisting}
183 184


185
\section{Session file configuration}
186
In the following we describe the session file configuration. Specifically we consider the
187 188
sections under the tag \inltt{<CONDITIONS>} in the session (.xml) file.
\subsection*{Parameters}
189 190 191 192 193 194 195 196 197 198 199 200
Under this section it is possible to set the parameters of the simulation.
\begin{lstlisting}[style=XmlStyle]
<PARAMETERS>
  <P> TimeStep            = 0.0000001                  </P>
  <P> FinTime             = 1.0                        </P>
  <P> NumSteps            = FinTime/TimeStep           </P>
  <P> IO_CheckSteps       = 5000                       </P>
  <P> IO_InfoSteps        = 1                          </P>
  <P> Gamma               = 1.4                        </P>
  <P> pInf                = 101325                     </P>
  <P> rhoInf              = 1.225                      </P>
  <P> TInf                = pInf/(287.058*rhoInf)      </P>
Douglas Serson's avatar
Douglas Serson committed
201
  <P> Twall               = pInf/(287.058*rhoInf)+15.0 </P>
202 203 204 205 206 207 208 209 210 211 212
  <P> uInf                = 147.4                      </P>
  <P> vInf                = 0.0                        </P>
  <P> wInf                = 0.0                        </P>
  <P> mu                  = 1e-5                       </P>
  <P> Pr                  = 0.72                       </P>
  <P> thermalConductivity = 0.02                       </P>
</PARAMETERS>
\end{lstlisting}
\begin{itemize}
\item \inltt{TimeStep} is the time-step we want to use;
\item \inltt{FinTime} is the final physical time at which we want our simulation to stop;
213
\item \inltt{NumSteps} is the equivalent of \inltt{FinTime} but instead of specifying the
214 215
physical final time we specify the number of time-steps;
\item \inltt{IO\_CheckSteps} sets the number of steps between successive checkpoint files;
216
\item \inltt{IO\_InfoSteps} sets the number of steps between successive info stats are printed
217 218 219 220 221
to screen;
\item \inltt{Gamma} ratio of the specific heats. Default value = 1.4;
\item \inltt{pInf} farfield pressure (i.e. $p_{\infty}$). Default value = 101325 $Pa$;
\item \inltt{rhoInf} farfield density (i.e. $\rho_{\infty}$). Default value = 1.225 $Kg/m^{3}$;
\item \inltt{TInf} farfield temperature (i.e. $T_{\infty}$). Default value = 288.15 $K$;
222
\item \inltt{Twall} temperature at the wall when isothermal boundary
Douglas Serson's avatar
Douglas Serson committed
223
conditions are employed (i.e. $T_{w}$). Default value = 300.15$K$;
224
\item \inltt{uint} farfield $X$-component of the velocity (i.e. $u_{\infty}$). Default value = 0.1 $m/s$;
225 226 227 228
\item \inltt{vInf} farfield $Y$-component of the velocity (i.e. $v_{\infty}$). Default value = 0.0 $m/s$;
\item \inltt{wInf} farfield $Z$-component of the velocity (i.e. $w_{\infty}$). Default value = 0.0 $m/s$;
\item \inltt{mu} dynamic viscosity (i.e. $\mu_{\infty}$). Default value = 1.78e-05 $Pa s$;
\item \inltt{Pr} Prandtl number. Default value = 0.72;
Douglas Serson's avatar
Douglas Serson committed
229
\item \inltt{thermalConductivity} thermal conductivity (i.e. $\kappa_{\infty}$). This can be set as an
230
 alternative to \inltt{Pr}, in which case the Prandtl number is calculated from $\kappa_{\infty}$
Douglas Serson's avatar
Douglas Serson committed
231
 (it is only possible to set one of them). By default, this is obtained from the Prandtl number;
232 233
\end{itemize}

234
\subsection*{Solver info}
235
Under this section it is possible to set the solver information.
236
\begin{lstlisting}[style=XmlStyle]
237
<SOLVERINFO>
238 239 240 241 242 243 244 245 246
  <I PROPERTY="EQType"                VALUE="NavierStokesCFE"      />
  <I PROPERTY="Projection"            VALUE="DisContinuous"        />
  <I PROPERTY="AdvectionType"         VALUE="WeakDG"               />
  <I PROPERTY="DiffusionType"         VALUE="LDGNS"                />
  <I PROPERTY="TimeIntegrationMethod" VALUE="ClassicalRungeKutta4" />
  <I PROPERTY="UpwindType"            VALUE="ExactToro"            />
  <I PROPERTY="ProblemType"           VALUE="General"              />
  <I PROPERTY="ViscosityType"         VALUE="Constant"             />
  <I PROPERTY="EquationOfState"       VALUE="IdealGas"             />
247 248 249 250 251 252 253
</SOLVERINFO>
\end{lstlisting}
\begin{itemize}
\item \inltt{EQType} is the tag which specify the equations we want solve:
\begin{itemize}
\item \inltt{NavierStokesCFE} (Compressible Navier-Stokes equations);
\item \inltt{EulerCFE} (Compressible Euler equations).
254 255
\item \inltt{IsentropicVortex} (Isentropic vortex test case).
\item \inltt{RinglebFlow} (Ringleb flow test case).
256 257 258 259
\end{itemize}
\item \inltt{Projection} is the type of projection we want to use:
\begin{itemize}
\item \inltt{DisContinuous}.\\
260
Note that the Continuous projection is not supported in the Compressible Flow Solver.
261
\end{itemize}
262
\item \inltt{AdvectionType} is the advection operator we want to use.
263 264 265 266 267 268 269 270
\begin{itemize}
\item \inltt{WeakDG} (classical DG in weak form);
\item \inltt{FRDG} (Flux-Reconstruction recovering nodal DG scheme);
\item \inltt{FRSD} (Flux-Reconstruction recovering a spectral difference (SD) scheme);
\item \inltt{FRHU} (Flux-Reconstruction recovering Huynh (G2) scheme);
\item \inltt{FRcmin} (Flux-Reconstruction with $c = c_{min}$);
\item \inltt{FRcinf} (Flux-Reconstruction with $c = \infty$).
\end{itemize}
271
Note that only \inltt{WeakDG} is fully supported, the other operators work only with quadrilateral elements ($2D$ or $2.5D$).
Douglas Serson's avatar
Douglas Serson committed
272 273
\item \inltt{DiffusionType} is the diffusion operator we want to use
for the Navier-Stokes equations:
274
\begin{itemize}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
275
\item \inltt{LDGNS} (LDG with primitive variables. The penalty term is inversely proportional to the element size, proportional to the local viscosity for the momentum equations and to the thermal conductivity for the energy equation, and proportional to an optional parameter \inltt{LDGNSc11} which is by default set to one; proportionality to polynomial order can be manually imposed by setting the parameter \inltt{LDGNSc11} equal to $p^2$);
Douglas Serson's avatar
Douglas Serson committed
276 277 278 279 280
\item \inltt{LFRDGNS} (Flux-Reconstruction recovering nodal DG scheme);
\item \inltt{LFRSDNS} (Flux-Reconstruction recovering a spectral difference (SD) scheme);
\item \inltt{LFRHUNS} (Flux-Reconstruction recovering Huynh (G2) scheme);
\item \inltt{LFRcminNS} (Flux-Reconstruction with $c = c_{min}$);
\item \inltt{LFRcinfNS} (Flux-Reconstruction with $c = \infty$).
281
\end{itemize}
282
Note that only \inltt{LDGNS} is fully supported, the other operators work only with quadrilateral elements ($2D$ or $2.5D$).
283
\item \inltt{TimeIntegrationMethod} is the time-integration scheme we want to use.
284 285 286
Note that only an explicit discretisation is supported:
\begin{itemize}
\item \inltt{ForwardEuler};
287 288
\item \inltt{RungeKutta2\_SSP};
\item \inltt{RungeKutta3\_SSP};
289 290
\item \inltt{ClassicalRungeKutta4}.
\end{itemize}
291
\item \inltt{UpwindType} is the numerical interface flux (i.e. Riemann solver)
292 293 294
we want to use for the advection operator:
\begin{itemize}
\item \inltt{AUSM0};
295 296 297 298
\item \inltt{AUSM1};
\item \inltt{AUSM2};
\item \inltt{AUSM3};
\item \inltt{Average};
299 300 301 302 303 304 305 306 307 308 309
\item \inltt{ExactToro};
\item \inltt{HLL};
\item \inltt{HLLC};
\item \inltt{LaxFriedrichs};
\item \inltt{Roe}.
\end{itemize}
\item \inltt{ViscosityType} is the viscosity type we want to use:
\begin{itemize}
\item \inltt{Constant} (Constant viscosity);
\item \inltt{Variable} (Variable viscosity through the Sutherland's law.);
\end{itemize}
310 311 312 313 314 315 316 317 318
\item \inltt{EquationOfState} allows selecting an equation of state for
accounting for non-ideal gas behaviour:
\begin{itemize}
\item \inltt{IdealGas} (default option);
\item \inltt{VanDerWaals} (requires additional parameters \inltt{Tcrit} and \inltt{Pcrit});
\item \inltt{RedlichKwong} (requires additional parameters \inltt{Tcrit} and \inltt{Pcrit});
\item \inltt{PengRobinson} (requires additional parameters \inltt{Tcrit}, \inltt{Pcrit} and
\inltt{AcentricFactor});
\end{itemize}
319 320
\end{itemize}

321
\subsection*{Boundary conditions}
322 323 324
In this section we can specify the boundary conditions for our problem.
First we need to define the variables under the section \inltt{VARIABLES}.
For a 1D problem we have:
325
\begin{lstlisting}[style=XmlStyle]
326 327 328 329 330 331 332 333
<VARIABLES>
  <V ID="0"> rho  </V>
  <V ID="1"> rhou </V>
  <V ID="4"> E    </V>
</VARIABLES>
\end{lstlisting}

For a 2D problem we have
334
\begin{lstlisting}[style=XmlStyle]
335 336 337 338 339 340 341 342 343
<VARIABLES>
  <V ID="0"> rho  </V>
  <V ID="1"> rhou </V>
  <V ID="2"> rhov </V>
  <V ID="4"> E    </V>
</VARIABLES>
\end{lstlisting}

For a 3D problem we have:
344
\begin{lstlisting}[style=XmlStyle]
345 346 347 348 349 350 351 352 353 354 355
<VARIABLES>
  <V ID="0"> rho  </V>
  <V ID="1"> rhou </V>
  <V ID="2"> rhov </V>
  <V ID="3"> rhow </V>
  <V ID="4"> E    </V>
</VARIABLES>
\end{lstlisting}

After having defined the variables depending on the dimensions of the problem we want to solve
it is necessary to specify the boundary regions on which we want to define the boundary conditions:
356
\begin{lstlisting}[style=XmlStyle]
357 358 359 360 361 362 363 364 365
<BOUNDARYREGIONS>
  <B ID="0"> C[100] </B>
</BOUNDARYREGIONS>
\end{lstlisting}

Finally we can specify the boundary conditions on the regions specified under \inltt{BOUNDARYREGIONS}.
In the following some examples for a 2D problem:
\begin{itemize}
\item Slip wall boundary conditions:
366
\begin{lstlisting}[style=XmlStyle]
367 368 369 370 371 372 373 374 375 376 377
<BOUNDARYCONDITIONS>
  <REGION REF="0">
    <D VAR="rho"  USERDEFINEDTYPE="Wall" VALUE="0" />
    <D VAR="rhou" USERDEFINEDTYPE="Wall" VALUE="0" />
    <D VAR="rhov" USERDEFINEDTYPE="Wall" VALUE="0" />
    <D VAR="E"    USERDEFINEDTYPE="Wall" VALUE="0" />
  </REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}

\item No-slip wall boundary conditions:
378
\begin{lstlisting}[style=XmlStyle]
379 380 381 382 383 384 385 386 387 388
<BOUNDARYCONDITIONS>
  <REGION REF="0">
    <D VAR="rho"  USERDEFINEDTYPE="WallViscous" VALUE="0" />
    <D VAR="rhou" USERDEFINEDTYPE="WallViscous" VALUE="0" />
    <D VAR="rhov" USERDEFINEDTYPE="WallViscous" VALUE="0" />
    <D VAR="E"    USERDEFINEDTYPE="WallViscous" VALUE="0" />
  </REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}

389 390 391 392 393 394 395 396 397 398 399 400
\item Adiabatic wall boundary conditions:
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
  <REGION REF="0">
    <D VAR="rho"  USERDEFINEDTYPE="WallAdiabatic" VALUE="0" />
    <D VAR="rhou" USERDEFINEDTYPE="WallAdiabatic" VALUE="0" />
    <D VAR="rhov" USERDEFINEDTYPE="WallAdiabatic" VALUE="0" />
    <D VAR="E"    USERDEFINEDTYPE="WallAdiabatic" VALUE="0" />
  </REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}

401
\item Farfield boundary conditions (including inviscid characteristic boundary conditions):
402 403
\begin{lstlisting}[style=XmlStyle]
<BOUNDARYCONDITIONS>
404 405 406 407
  <REGION REF="0">
    <D VAR="rho"  VALUE="rhoInf" />
    <D VAR="rhou" VALUE="rhoInf*uInf" />
    <D VAR="rhov" VALUE="rhoInf*vInf" />
408
    <D VAR="E"
409 410 411 412 413 414
    VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInf*uInf+vInf*vInf+wInf*wInf)"/>
  </REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}

\item Pressure outflow boundary conditions:
415
\begin{lstlisting}[style=XmlStyle]
416 417 418 419 420
<BOUNDARYCONDITIONS>
  <REGION REF="0">
    <D VAR="rho"  USERDEFINEDTYPE="PressureOutflow" VALUE="0" />
    <D VAR="rhou" USERDEFINEDTYPE="PressureOutflow" VALUE="0" />
    <D VAR="rhov" USERDEFINEDTYPE="PressureOutflow" VALUE="0" />
421
    <D VAR="E"    USERDEFINEDTYPE="PressureOutflow" VALUE="pOut" />
422 423 424
  </REGION>
</BOUNDARYCONDITIONS>
\end{lstlisting}
425
where pOut is the target static pressure at the boundary.
426 427
\end{itemize}

428
\subsection*{Initial conditions and exact solution}
429
Under the two following sections it is possible to define the initial conditions and the exact solution (if existent).
430
\begin{lstlisting}[style=XmlStyle]
431 432 433 434
<FUNCTION NAME="InitialConditions">
  <E VAR="rho"    VALUE="rhoInf"/>
  <E VAR="rhou"   VALUE="rhoInf*uInf"   />
  <E VAR="rhov"   VALUE="rhoInf*vInf"   />
435
  <E VAR="E"
436 437
  VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInf*uInf+vInf*vInf+wInf*wInf)"/>
</FUNCTION>
438

439 440 441 442
<FUNCTION NAME="ExactSolution">
  <E VAR="rho"    VALUE="rhoInf"          />
  <E VAR="rhou"   VALUE="rhoInf*uInf"   />
  <E VAR="rhov"   VALUE="rhoInf*vInf"   />
443
  <E VAR="E"
444 445 446
  VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInf*uInf+vInf*vInf+wInf*wInf)"/>
</FUNCTION>
\end{lstlisting}
447

448 449
\section{Examples}
\subsection{Shock capturing}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
450
Compressible flows can be characterised by abrupt changes in density 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.
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
451

452
\subsubsection{Non-smooth artificial viscosity model}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
453 454 455 456 457
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 \ref{eq:euler} \cite{persson2006sub}.
The diffusivity of the system is controlled by a variable viscosity coefficient $\varepsilon$.
For consistency $\varepsilon$ is proportional to the element size and inversely proportional to the polynomial order.
Finally, from physical considerations $\varepsilon$ needs to be proportional to the maximum characteristic speed of the problem.
The final form of the artificial viscosity is
458
\begin{equation}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477
  \varepsilon = \varepsilon_0 \frac{h}{p} \lambda_{max} S,
\end{equation}
where $S$ is a sensor.

As shock sensor, a modal resolution-based indicator is used
\begin{equation}\label{eq:sensor}
  s_e = log_{10}\left( \frac{\langle q - \tilde{q}, q - \tilde{q}  \rangle}{\langle q, q \rangle} \right) ,
\end{equation}
where $\langle \cdot, \cdot \rangle$ represents a $L^2$ inner product, $q$ and $\tilde{q}$ are the full and truncated expansions of a state variable (in our case density)
\begin{equation}
  q(x) = \sum_{i=1}^{N(P)} \hat{q}_i \phi_i , \quad \tilde{q}(x) = \sum_{i=1}^{N(P-1)} \hat{q}_i \phi_i ,
\end{equation}
then the constant element-wise sensor is computed as follows
\begin{equation}
   S_\varepsilon
=  \left \{ \begin{array}{lll}
    0 &	 \mbox{if} & s_e<s_0-\kappa \\
    \frac{1}{2}\left(1+\sin{\frac{\pi\left(s_e-s_0\right)}{2\kappa}}\right) &	\mbox{if} &	| s_e - s_0| \le \kappa\\
    1 &	 \mbox{if} & s_e > s_0+\kappa
478
      \end{array}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
479
    \right.,
480
\end{equation}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
481
where $s_0 = s_\kappa - 4.25\;log_{10}(p)$.
482 483 484

To enable the non-smooth viscosity model, the following line has to be added to the \inltt{SOLVERINFO} section:
\begin{lstlisting}[style=XmlStyle]
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
485
<SOLVERINFO>
486
  <I PROPERTY="ShockCaptureType"      VALUE="NonSmooth" />
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
487
<SOLVERINFO>
488
\end{lstlisting}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
489
The diffusivity and the sensor can be controlled by the following parameters:
490 491 492
\begin{lstlisting}[style=XmlStyle]
<PARAMETERS>
<P> Skappa 	 	= -1.3                     </P>
493
<P> Kappa 	 	= 0.2                      </P>
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
494
<P> mu0 	  	= 1.0                      </P>
495 496
</PARAMETERS>
\end{lstlisting}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
497

498 499
\begin{figure}[!htbp]
\begin{center}
500 501
\includegraphics[width = 0.47 \textwidth]{img/Mach_P4.pdf}
\includegraphics[width = 0.47 \textwidth]{img/ArtVisc_P4.pdf}
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
502
\caption{(a) Steady state solution for $M=0.8$ flow at $\alpha = 1.25^\circ$ past a NACA 0012 profile, (b) Artificial viscosity ($\varepsilon$) distribution}
503 504 505
\label{fig:}
\end{center}
\end{figure}
506

507
\subsection{Variable polynomial order}
508
A sensor based $p$-adaptive algorithm is implemented to optimise the computational cost and accuracy.
509 510 511 512 513
The DG scheme allows one to use different polynomial orders since the fluxes over the elements are determined using a Riemann solver and there is now 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.\\
\begin{equation}\label{eq:pswitch}
  p_e =\left \{ \begin{array}{l}
514
    p_e-1\	\	\ \mbox{if}\		\	 s_e>s_{ds}\\
515 516 517
    p_e+1\	\	\ \mbox{if}\		\	 s_{sm }<s_e<s_{ds}\\
    p_e\	\	\	\	\	\	\	\	\	 \mbox{if}\		\ s_{fl}<s_e<s_{sm}\\
    p_e-1\	\	\ \mbox{if}\		\	 s_e<s_{fl}
518
    \end{array}
519 520 521
    \right.
\end{equation}
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.
522
\subsection{De-Aliasing Techniques}
523 524 525 526
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:
527 528 529 530 531
\begin{itemize}
\item PDE nonlinearities, related to the nonlinear and quasi-linear fluxes.
\item Geometrical nonlinearities, related to the deformed/curves meshes.
\end{itemize}
We consider two de-aliasing strategies based on the concept of consistent integration:
532

533
\begin{itemize}
534 535
\item Local dealiasing: It only targets the PDE-aliasing sources, applying a consistent integration of them locally.
\item 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.
536
\end{itemize}
Douglas Serson's avatar
Douglas Serson committed
537

538 539
Since Nektar++ tackles separately the PDE and geometric aliasing during the
projection and solution of the equations, to consistently
540
integrate all the nonlinearities in the compressible
541 542
NavierStokes equations, the quadrature points should
be selected based on the maximum order of the nonlinearities:
543 544 545
\begin{equation}
Q_{min}= P_{exp}+\frac{max(2P_{exp},P_{geom})}{2} + \frac{3}{2}
\end{equation}
Douglas Serson's avatar
Douglas Serson committed
546

547
where $Q_{min}$ is the minimum required number of quadrature
548 549
points to exactly integrate the highest-degree of nonlinearity,
$P_{exp}$ being the order of the polynomial expansion and $P_{geom}$
Giacomo Castiglioni's avatar
Giacomo Castiglioni committed
550
being the geometric order of the mesh. Bear in mind that we are
551
using a discontinuous discretisation, meaning that aliasing
552
effect are not fully controlled, since the boundary terms
553 554 555 556
introduce non-polynomial functions into the problem.

To enable the global de-aliasing technique, modify the number of quadrature
points by:
Douglas Serson's avatar
Douglas Serson committed
557 558

\begin{lstlisting}[style=XmlStyle]
559 560 561 562 563
<E COMPOSITE="[101]"
   BASISTYPE="Modified_A,Modified_A"
   NUMMODES="7,7"
   POINTSTYPE="GaussLobattoLegendre,GaussLobattoLegendre"
   NUMPOINTS="14,14"
564 565
   FIELDS="rho,rhou,rhov,E"
/>
Douglas Serson's avatar
Douglas Serson committed
566
\end{lstlisting}
567

568
where \inltt{NUMMODES} corresponds to $P$+1, where $P$ is the order of the polynomial
569
used to approximate the solution. \inltt{NUMPOINTS} specifies the number of quadrature
570
points.