compressible-flow.tex 25.3 KB
 1 \chapter{Compressible Flow Solver}  Chris Cantwell committed Aug 08, 2014 2   3 \section{Synopsis}  Gianmarco Mengaldo committed Aug 10, 2014 4 The CompressibleFlowSolver allows us to solve  Giacomo Castiglioni committed Jan 10, 2019 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  Gianmarco Mengaldo committed Aug 10, 2014 8 9 both the compressible Euler and the Navier-Stokes equations.  10 \subsection{Euler equations}  Giacomo Castiglioni committed Jan 10, 2019 11 12 The Euler equations can be expressed as a hyperbolic conservation law in the form  Dirk Ekelschot committed Aug 14, 2014 13 \label{eq:euler}  Giacomo Castiglioni committed Jan 10, 2019 14 15 \frac{\partial \mathbf{q} }{\partial t} + \frac{\partial \mathbf{f}_i}{\partial x} + \frac{\partial \mathbf{g}_i}{\partial y} +  Gianmarco Mengaldo committed Aug 10, 2014 16 17 \frac{\partial \mathbf{h}_i}{\partial z} = 0,  Giacomo Castiglioni committed Jan 10, 2019 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  Gianmarco Mengaldo committed Aug 10, 2014 22 23 24 25 26 27 inviscid fluxes \mathbf{q} =\left\{ \begin{array}{c} \rho \\ \rho u \\  Giacomo Castiglioni committed Jan 10, 2019 28 \rho v \\  Gianmarco Mengaldo committed Aug 10, 2014 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 \\  Giacomo Castiglioni committed Jan 10, 2019 36 \rho u v \\  Gianmarco Mengaldo committed Aug 10, 2014 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 \\  Giacomo Castiglioni committed Jan 10, 2019 44 p + \rho v^2 \\  Gianmarco Mengaldo committed Aug 10, 2014 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 \\  Giacomo Castiglioni committed Jan 10, 2019 52 \rho v w \\  Gianmarco Mengaldo committed Aug 10, 2014 53 54 55 56 57 p + \rho w^2 \\ w (E + p) \end{array} \right\}, \label{Euler}  Giacomo Castiglioni committed Jan 10, 2019 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  Gianmarco Mengaldo committed Aug 10, 2014 60 61 a perfect gas law for which the pressure is related to the total energy by the following expression  Giacomo Castiglioni committed Jan 10, 2019 62 E = \frac{p}{\gamma - 1} + \frac{1}{2} \rho(u^2 + v^2 + w^2),  Gianmarco Mengaldo committed Aug 10, 2014 63 64 65 \label{energy} where $\gamma$ is the ratio of specific heats.  66 \subsection{Compressible Navier-Stokes equations}  Giacomo Castiglioni committed Jan 10, 2019 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  Gianmarco Mengaldo committed Aug 10, 2014 71 72 three-dimensional case can be written as:  Giacomo Castiglioni committed Jan 10, 2019 73 74 75 \frac{\partial \mathbf{q} }{\partial t} + \frac{\partial \mathbf{f}}{\partial x} + \frac{\partial \mathbf{g}}{\partial y} +  Gianmarco Mengaldo committed Aug 10, 2014 76 77 \frac{\partial \mathbf{h}}{\partial z} = 0,  Giacomo Castiglioni committed Jan 10, 2019 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}))$  Gianmarco Mengaldo committed Aug 10, 2014 82 83 84 85 86 are the vectors of the fluxes which can also be written as: \begin{array}{l} \mathbf{f} = \mathbf{f}_i - \mathbf{f}_v, \mathbf{g} = \mathbf{g}_i - \mathbf{g}_v,  Giacomo Castiglioni committed Jan 10, 2019 87 \mathbf{h} = \mathbf{h}_i - \mathbf{h}_v,  Gianmarco Mengaldo committed Aug 10, 2014 88 89 \end{array}  Giacomo Castiglioni committed Jan 10, 2019 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  Gianmarco Mengaldo committed Aug 10, 2014 92 93 94 95 96 97 98 the viscous fluxes which take the following form: \begin{array}{c} \vspace{0.2 cm} \mathbf{f}_v =\left\{ \begin{array}{c} 0 \\ \tau_{xx} \\  Giacomo Castiglioni committed Jan 10, 2019 99 \tau_{yx} \\  Gianmarco Mengaldo committed Aug 10, 2014 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} \\  Giacomo Castiglioni committed Jan 10, 2019 107 \tau_{yy} \\  Gianmarco Mengaldo committed Aug 10, 2014 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} \\  Giacomo Castiglioni committed Jan 10, 2019 115 \tau_{yz} \\  Gianmarco Mengaldo committed Aug 10, 2014 116 117 118 \tau_{zz} \\ u \tau_{xz} + v \tau_{yz} + w \tau_{zz} + k T_z \end{array} \right\},  Giacomo Castiglioni committed Jan 10, 2019 119 \end{array}  Gianmarco Mengaldo committed Aug 10, 2014 120 121 \label{viscous flux}  Giacomo Castiglioni committed Jan 10, 2019 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$.}  Gianmarco Mengaldo committed Aug 10, 2014 126 127  \begin{array}{cc} \vspace{0.2 cm}  Giacomo Castiglioni committed Jan 10, 2019 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} =  Gianmarco Mengaldo committed Aug 10, 2014 133 134 135 \mu(u_z + w_x). \\ \end{array}  Giacomo Castiglioni committed Jan 10, 2019 136 where $\mu$ is the dynamic viscosity calculated using  Gianmarco Mengaldo committed Aug 10, 2014 137 138 the Sutherland's law and $k$ is the thermal conductivity.  139 \subsection{Numerical discretisation}  Giacomo Castiglioni committed Jan 10, 2019 140 In Nektar++ the spatial discretisation of the Euler and of the Navier-Stokes  Gianmarco Mengaldo committed Aug 10, 2014 141 equations is projected in the polynomial space via a discontinuous projection.  Giacomo Castiglioni committed Jan 10, 2019 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.  Gianmarco Mengaldo committed Aug 10, 2014 153   Giacomo Castiglioni committed Jan 10, 2019 154 For the advection term it is common to solve a Riemann problem at each  Gianmarco Mengaldo committed Aug 10, 2014 155 interface of the computational domain through exact or approximated Riemann  Giacomo Castiglioni committed Jan 10, 2019 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  Gianmarco Mengaldo committed Aug 10, 2014 168 169 cost. The approximated Riemann solvers are simplifications of the exact solver.  Giacomo Castiglioni committed Jan 10, 2019 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  Gianmarco Mengaldo committed Aug 10, 2014 172 173 174 FR diffusion terms. The boundary conditions are also implemented by exploiting the numerical  Giacomo Castiglioni committed Jan 10, 2019 175 176 interface fluxes just mentioned. For a more detailed description of the above the interested reader can refer  Gianmarco Mengaldo committed Aug 10, 2014 177 178 to \cite{DeGMen14} and \cite{MenDeG14}.  179 \section{Usage}  Chris Cantwell committed Aug 19, 2014 180 181 182 \begin{lstlisting}[style=BashInputStyle] CompressibleFlowSolver session.xml \end{lstlisting}  Gianmarco Mengaldo committed Aug 10, 2014 183 184   185 \section{Session file configuration}  Gianmarco Mengaldo committed Aug 10, 2014 186 In the following we describe the session file configuration. Specifically we consider the  187 188 sections under the tag \inltt{} in the session (.xml) file. \subsection*{Parameters}  Gianmarco Mengaldo committed Aug 10, 2014 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]

TimeStep = 0.0000001

FinTime = 1.0

NumSteps = FinTime/TimeStep

IO_CheckSteps = 5000

IO_InfoSteps = 1

Gamma = 1.4

pInf = 101325

rhoInf = 1.225

TInf = pInf/(287.058*rhoInf)

 Douglas Serson committed Jul 21, 2016 201 

Twall = pInf/(287.058*rhoInf)+15.0

 Gianmarco Mengaldo committed Aug 10, 2014 202 203 204 205 206 207 208 209 210 211 212 

uInf = 147.4

vInf = 0.0

wInf = 0.0

mu = 1e-5

Pr = 0.72

thermalConductivity = 0.02

\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;  Giacomo Castiglioni committed Jan 10, 2019 213 \item \inltt{NumSteps} is the equivalent of \inltt{FinTime} but instead of specifying the  Gianmarco Mengaldo committed Aug 10, 2014 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;  Giacomo Castiglioni committed Jan 10, 2019 216 \item \inltt{IO\_InfoSteps} sets the number of steps between successive info stats are printed  Gianmarco Mengaldo committed Aug 10, 2014 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$;  Giacomo Castiglioni committed Jan 10, 2019 222 \item \inltt{Twall} temperature at the wall when isothermal boundary  Douglas Serson committed Jul 21, 2016 223 conditions are employed (i.e. $T_{w}$). Default value = 300.15$K$;  Chris Cantwell committed Sep 24, 2017 224 \item \inltt{uint} farfield $X$-component of the velocity (i.e. $u_{\infty}$). Default value = 0.1 $m/s$;  Gianmarco Mengaldo committed Aug 10, 2014 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 committed Sep 06, 2016 229 \item \inltt{thermalConductivity} thermal conductivity (i.e. $\kappa_{\infty}$). This can be set as an  Giacomo Castiglioni committed Jan 10, 2019 230  alternative to \inltt{Pr}, in which case the Prandtl number is calculated from $\kappa_{\infty}$  Douglas Serson committed Sep 06, 2016 231  (it is only possible to set one of them). By default, this is obtained from the Prandtl number;  Gianmarco Mengaldo committed Aug 10, 2014 232 233 \end{itemize}  234 \subsection*{Solver info}  Gianmarco Mengaldo committed Aug 10, 2014 235 Under this section it is possible to set the solver information.  Giacomo Castiglioni committed Jan 10, 2019 236 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 237   Andrea Cassinelli committed Sep 20, 2019 238 239 240 241 242 243 244 245 246   Gianmarco Mengaldo committed Aug 10, 2014 247 248 249 250 251 252 253  \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).  Douglas Serson committed Jul 20, 2016 254 255 \item \inltt{IsentropicVortex} (Isentropic vortex test case). \item \inltt{RinglebFlow} (Ringleb flow test case).  Gianmarco Mengaldo committed Aug 10, 2014 256 257 258 259 \end{itemize} \item \inltt{Projection} is the type of projection we want to use: \begin{itemize} \item \inltt{DisContinuous}.\\  Giacomo Castiglioni committed Jan 10, 2019 260 Note that the Continuous projection is not supported in the Compressible Flow Solver.  Gianmarco Mengaldo committed Aug 10, 2014 261 \end{itemize}  Giacomo Castiglioni committed Sep 19, 2019 262 \item \inltt{AdvectionType} is the advection operator we want to use.  Gianmarco Mengaldo committed Aug 10, 2014 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}  Giacomo Castiglioni committed Sep 19, 2019 271 Note that only \inltt{WeakDG} is fully supported, the other operators work only with quadrilateral elements ($2D$ or $2.5D$).  Douglas Serson committed Sep 13, 2016 272 273 \item \inltt{DiffusionType} is the diffusion operator we want to use for the Navier-Stokes equations:  Gianmarco Mengaldo committed Aug 10, 2014 274 \begin{itemize}  Giacomo Castiglioni committed Nov 07, 2019 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 committed Sep 13, 2016 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$).  Gianmarco Mengaldo committed Aug 10, 2014 281 \end{itemize}  Giacomo Castiglioni committed Sep 19, 2019 282 Note that only \inltt{LDGNS} is fully supported, the other operators work only with quadrilateral elements ($2D$ or $2.5D$).  Giacomo Castiglioni committed Jan 10, 2019 283 \item \inltt{TimeIntegrationMethod} is the time-integration scheme we want to use.  Gianmarco Mengaldo committed Aug 10, 2014 284 285 286 Note that only an explicit discretisation is supported: \begin{itemize} \item \inltt{ForwardEuler};  Giacomo Castiglioni committed Sep 19, 2019 287 288 \item \inltt{RungeKutta2\_SSP}; \item \inltt{RungeKutta3\_SSP};  Gianmarco Mengaldo committed Aug 10, 2014 289 290 \item \inltt{ClassicalRungeKutta4}. \end{itemize}  Giacomo Castiglioni committed Jan 10, 2019 291 \item \inltt{UpwindType} is the numerical interface flux (i.e. Riemann solver)  Gianmarco Mengaldo committed Aug 10, 2014 292 293 294 we want to use for the advection operator: \begin{itemize} \item \inltt{AUSM0};  Giacomo Castiglioni committed Jan 10, 2019 295 296 297 298 \item \inltt{AUSM1}; \item \inltt{AUSM2}; \item \inltt{AUSM3}; \item \inltt{Average};  Gianmarco Mengaldo committed Aug 10, 2014 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}  Douglas Serson committed Jan 15, 2018 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}  Gianmarco Mengaldo committed Aug 10, 2014 319 320 \end{itemize}  321 \subsection*{Boundary conditions}  Gianmarco Mengaldo committed Aug 10, 2014 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:  Giacomo Castiglioni committed Jan 10, 2019 325 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 326 327 328 329 330 331 332 333  rho rhou E \end{lstlisting} For a 2D problem we have  Giacomo Castiglioni committed Jan 10, 2019 334 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 335 336 337 338 339 340 341 342 343  rho rhou rhov E \end{lstlisting} For a 3D problem we have:  Giacomo Castiglioni committed Jan 10, 2019 344 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 345 346 347 348 349 350 351 352 353 354 355  rho rhou rhov rhow E \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:  Giacomo Castiglioni committed Jan 10, 2019 356 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 357 358 359 360 361 362 363 364 365  C[100] \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:  Dave Moxey committed Jul 23, 2015 366 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 367 368 369 370 371 372 373 374 375 376 377  \end{lstlisting} \item No-slip wall boundary conditions:  Dave Moxey committed Jul 23, 2015 378 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 379 380 381 382 383 384 385 386 387 388  \end{lstlisting}  Dave Moxey committed Jul 23, 2015 389 390 391 392 393 394 395 396 397 398 399 400 \item Adiabatic wall boundary conditions: \begin{lstlisting}[style=XmlStyle] \end{lstlisting}  Gianmarco Mengaldo committed Aug 10, 2014 401 \item Farfield boundary conditions (including inviscid characteristic boundary conditions):  Giacomo Castiglioni committed Jan 10, 2019 402 403 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 404 405 406 407   Giacomo Castiglioni committed Jan 10, 2019 408  \end{lstlisting} \item Pressure outflow boundary conditions:  Giacomo Castiglioni committed Jan 10, 2019 415 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 416 417 418 419 420   Giacomo Castiglioni committed Jan 10, 2019 421   Gianmarco Mengaldo committed Aug 10, 2014 422 423 424  \end{lstlisting}  Giacomo Castiglioni committed Jan 10, 2019 425 where pOut is the target static pressure at the boundary.  Gianmarco Mengaldo committed Aug 10, 2014 426 427 \end{itemize}  428 \subsection*{Initial conditions and exact solution}  Gianmarco Mengaldo committed Aug 10, 2014 429 Under the two following sections it is possible to define the initial conditions and the exact solution (if existent).  Giacomo Castiglioni committed Jan 10, 2019 430 \begin{lstlisting}[style=XmlStyle]  Gianmarco Mengaldo committed Aug 10, 2014 431 432 433 434   Giacomo Castiglioni committed Jan 10, 2019 435   Giacomo Castiglioni committed Jan 10, 2019 438   Gianmarco Mengaldo committed Aug 10, 2014 439 440 441 442   Giacomo Castiglioni committed Jan 10, 2019 443  \end{lstlisting}  Chris Cantwell committed Aug 19, 2014 447   448 449 \section{Examples} \subsection{Shock capturing}  Giacomo Castiglioni committed Sep 19, 2019 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 committed May 02, 2019 451   452 \subsubsection{Non-smooth artificial viscosity model}  Giacomo Castiglioni committed May 02, 2019 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  Dirk Ekelschot committed Aug 14, 2014 458   Giacomo Castiglioni committed May 02, 2019 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, where $S$ is a sensor. As shock sensor, a modal resolution-based indicator is used \label{eq:sensor} s_e = log_{10}\left( \frac{\langle q - \tilde{q}, q - \tilde{q} \rangle}{\langle q, q \rangle} \right) , 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) 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 , then the constant element-wise sensor is computed as follows S_\varepsilon = \left \{ \begin{array}{lll} 0 & \mbox{if} & s_e s_0+\kappa  Giacomo Castiglioni committed Jan 10, 2019 478  \end{array}  Giacomo Castiglioni committed May 02, 2019 479  \right.,  Dirk Ekelschot committed Aug 14, 2014 480   Giacomo Castiglioni committed May 02, 2019 481 where $s_0 = s_\kappa - 4.25\;log_{10}(p)$.  Chris Cantwell committed Sep 24, 2017 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 committed Sep 19, 2019 485   Andrea Cassinelli committed Sep 20, 2019 486   Giacomo Castiglioni committed Sep 19, 2019 487   Chris Cantwell committed Sep 24, 2017 488 \end{lstlisting}  Giacomo Castiglioni committed May 02, 2019 489 The diffusivity and the sensor can be controlled by the following parameters:  Chris Cantwell committed Sep 24, 2017 490 491 492 \begin{lstlisting}[style=XmlStyle]

Skappa = -1.3

 Giacomo Castiglioni committed Jan 10, 2019 493 

Kappa = 0.2

 Giacomo Castiglioni committed May 02, 2019 494 

mu0 = 1.0

 Chris Cantwell committed Sep 24, 2017 495 496 
\end{lstlisting}  Giacomo Castiglioni committed May 02, 2019 497   Dirk Ekelschot committed Aug 14, 2014 498 499 \begin{figure}[!htbp] \begin{center}  Chris Cantwell committed Jul 25, 2015 500 501 \includegraphics[width = 0.47 \textwidth]{img/Mach_P4.pdf} \includegraphics[width = 0.47 \textwidth]{img/ArtVisc_P4.pdf}  Giacomo Castiglioni committed May 02, 2019 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}  Dirk Ekelschot committed Aug 14, 2014 503 504 505 \label{fig:} \end{center} \end{figure}  Chris Cantwell committed Sep 24, 2017 506   507 \subsection{Variable polynomial order}  Giacomo Castiglioni committed Jan 10, 2019 508 A sensor based $p$-adaptive algorithm is implemented to optimise the computational cost and accuracy.  Dirk Ekelschot committed Aug 14, 2014 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.\\ \label{eq:pswitch} p_e =\left \{ \begin{array}{l}  Giacomo Castiglioni committed Jan 10, 2019 514  p_e-1\ \ \ \mbox{if}\ \ s_e>s_{ds}\\  Dirk Ekelschot committed Aug 14, 2014 515 516 517  p_e+1\ \ \ \mbox{if}\ \ s_{sm }  Douglas Serson committed Apr 13, 2017 566 \end{lstlisting}  Chris Cantwell committed Aug 08, 2014 567   Chris Cantwell committed Sep 24, 2017 568 where \inltt{NUMMODES} corresponds to $P$+1, where $P$ is the order of the polynomial  Giacomo Castiglioni committed Jan 10, 2019 569 used to approximate the solution. \inltt{NUMPOINTS} specifies the number of quadrature  Chris Cantwell committed Sep 24, 2017 570 points.