Compressible NavierStokes miniapp¶
This example is located in the subdirectory examples/fluids
.
It solves the timedependent NavierStokes equations of compressible gas dynamics in a static Eulerian threedimensional frame using unstructured highorder finite/spectral element spatial discretizations and explicit or implicit highorder timestepping (available in PETSc).
Moreover, the NavierStokes example has been developed using PETSc, so that the pointwise physics (defined at quadrature points) is separated from the parallelization and meshing concerns.
Running the miniapp¶
The NavierStokes miniapp is controlled via commandline options. The following options are common among all problem types:
Option 
Description 
Default value 


CEED resource specifier 


Run in test mode 


Test absolute tolerance 


Test filename 


Problem to solve ( 


Use implicit time integartor formulation 


Polynomial degree of tensor product basis (must be >= 1) 


Number of extra quadrature points 


Use regular refinement for visualization 


Frequency of output, in number of steps. 


Continue from previous solution 


Output directory 


Use wall boundary conditions on this list of faces 


An array of constrained component numbers for wall BCs 


Use slip boundary conditions, for the x component, on this list of faces 


Use slip boundary conditions, for the y component, on this list of faces 


Use slip boundary conditions, for the z component, on this list of faces 


Use inflow boundary conditions on this list of faces 


Use outflow boundary conditions on this list of faces 


View PETSc 


View PETSc performance log 


View comprehensive information about runtime options 
For the case of a square/cubic mesh, the list of face indices to be used with bc_wall
, bc_inflow
, bc_outflow
and/or bc_slip_x
, bc_slip_y
, and bc_slip_z
are:
PETSc Face Name 
Cartesian direction 
Face ID 

faceMarkerBottom 
z 
1 
faceMarkerRight 
+x 
2 
faceMarkerTop 
+z 
3 
faceMarkerLeft 
x 
4 
PETSc Face Name 
Cartesian direction 
Face ID 

faceMarkerBottom 
z 
1 
faceMarkerTop 
+z 
2 
faceMarkerFront 
y 
3 
faceMarkerBack 
+y 
4 
faceMarkerRight 
+x 
5 
faceMarkerLeft 
x 
6 
Advection¶
For testing purposes, there is a reduced mode for pure advection, which holds density \(\rho\) and momentum density \(\rho \bm u\) constant while advecting “total energy density” \(E\). These are available in 2D and 3D.
2D advection¶
For the 2D advection problem, the following additional commandline options are available:
Option 
Description 
Default value 
Unit 


Characteristic radius of thermal bubble 



1 meter in scaled length units 



1 second in scaled time units 



1 kilogram in scaled mass units 



Strong (1) or weak/integrated by parts (0) residual 



Stabilization method ( 



Scale coefficient for stabilization tau (nondimensional) 



Wind type in Advection ( 



Constant wind vector when 



Total energy of inflow wind when 


An example of the rotation
mode can be run with:
./navierstokes problem advection2d dm_plex_box_faces 20,20 dm_plex_box_lower 0,0 dm_plex_box_upper 1000,1000 bc_wall 1,2,3,4 wall_comps 4 wind_type rotation implicit stab supg
and the translation
mode with:
./navierstokes problem advection2d dm_plex_box_faces 20,20 dm_plex_box_lower 0,0 dm_plex_box_upper 1000,1000 units_meter 1e4 wind_type translation wind_translation 1,.5 bc_inflow 1,2,3,4
Note the lengths in dm_plex_box_upper
are given in meters, and will be nondimensionalized according to units_meter
.
3D advection¶
For the 3D advection problem, the following additional commandline options are available:
Option 
Description 
Default value 
Unit 


Characteristic radius of thermal bubble 



1 meter in scaled length units 



1 second in scaled time units 



1 kilogram in scaled mass units 



Strong (1) or weak/integrated by parts (0) residual 



Stabilization method ( 



Scale coefficient for stabilization tau (nondimensional) 



Wind type in Advection ( 



Constant wind vector when 



Total energy of inflow wind when 









An example of the rotation
mode can be run with:
./navierstokes problem advection dm_plex_box_faces 10,10,10 dm_plex_dim 3 dm_plex_box_lower 0,0,0 dm_plex_box_upper 8000,8000,8000 bc_wall 1,2,3,4,5,6 wall_comps 4 wind_type rotation implicit stab su
and the translation
mode with:
./navierstokes problem advection dm_plex_box_faces 10,10,10 dm_plex_dim 3 dm_plex_box_lower 0,0,0 dm_plex_box_upper 8000,8000,8000 wind_type translation wind_translation .5,1,0 bc_inflow 1,2,3,4,5,6
Inviscid Ideal Gas¶
Isentropic Euler vortex¶
For the Isentropic Vortex problem, the following additional commandline options are available:
Option 
Description 
Default value 
Unit 


Location of vortex center 



1 meter in scaled length units 



1 second in scaled time units 



Background velocity vector 



Strength of vortex < 10 



Stabilization constant 

This problem can be run with:
./navierstokes problem euler_vortex dm_plex_box_faces 20,20,1 dm_plex_box_lower 0,0,0 dm_plex_box_upper 1000,1000,50 dm_plex_dim 3 bc_inflow 4,6 bc_outflow 3,5 bc_slip_z 1,2 mean_velocity .5,.8,0.
Sod shock tube¶
For the Shock Tube problem, the following additional commandline options are available:
Option 
Description 
Default value 
Unit 


1 meter in scaled length units 



1 second in scaled time units 



Use YZB discontinuity capturing 



Stabilization method ( 

This problem can be run with:
./navierstokes problem shocktube yzb stab su bc_slip_z 3,4 bc_slip_y 1,2 bc_wall 5,6 dm_plex_dim 3 dm_plex_box_lower 0,0,0 dm_plex_box_upper 1000,100,100 dm_plex_box_faces 200,1,1 units_second 0.1
Newtonian viscosity, Ideal Gas¶
For the Density Current, Channel, and Blasius problems, the following common commandline options are available:
Option 
Description 
Default value 
Unit 


1 meter in scaled length units 



1 second in scaled time units 



1 kilogram in scaled mass units 



1 Kelvin in scaled temperature units 



Stabilization method ( 



Stabilization constant, \(c_\tau\) 



Stabilization time constant, \(C_t\) 



Stabilization viscous constant, \(C_v\) 



Stabilization continuity constant, \(C_c\) 



Stabilization momentum constant, \(C_m\) 



Stabilization energy constant, \(C_E\) 



Heat capacity at constant volume 



Heat capacity at constant pressure 



Gravitational acceleration 



Stokes hypothesis second viscosity coefficient 



Shear dynamic viscosity coefficient 



Thermal conductivity 



Developer option to test properties 

boolean 

State variables to solve solution with. 

string 
Density current¶
The Density Current problem the following commandline options are available in addition to the Newtonian Ideal Gas options:
Option 
Description 
Default value 
Unit 


Location of bubble center 



Axis of density current cylindrical anomaly, or 



Characteristic radius of thermal bubble 



Reference potential temperature 



Perturbation of potential temperature 



Atmospheric pressure 



BruntVaisala frequency 


This problem can be run with:
./navierstokes problem density_current dm_plex_box_faces 16,1,8 degree 1 dm_plex_box_lower 0,0,0 dm_plex_box_upper 2000,125,1000 dm_plex_dim 3 rc 400. bc_wall 1,2,5,6 wall_comps 1,2,3 bc_slip_y 3,4 mu 75
Channel flow¶
The Channel problem the following commandline options are available in addition to the Newtonian Ideal Gas options:
Option 
Description 
Default value 
Unit 


Maximum/centerline velocity of the flow 



Reference potential temperature 



Atmospheric pressure 



Multiplier for body force ( 
1 
This problem can be run with the channel.yaml
file via:
./navierstokes options_file channel.yaml
problem: 'channel'
mu: .01
umax: 40
implicit: true
ts:
type: 'beuler'
adapt_type: 'none'
dt: 5e6
dm_plex_box_lower: 0,0,0
dm_plex_box_upper: .01,.01,.001
dm_plex_dim: 3
degree: 1
dm_plex_box_faces: 10,10,1
bc_slip_z: 1,2
bc_wall: 3,4
wall_comps: 1,2,3
dm_plex_box_bd: 'periodic,none,none'
Blasius boundary layer¶
The Blasius problem the following commandline options are available in addition to the Newtonian Ideal Gas options:
Option 
Description 
Default value 
Unit 


Freestream velocity 



Freestream temperature 



Wall temperature 



Boundary layer height at the inflow 



Atmospheric pressure 



Height at which 



Number of elements to keep below 



Growth rate of the elements in the refinement region 



Downward angle of the top face of the domain. This face serves as an outlet. 



Whether to use stg for the inflow conditions 



Path to file with y node locations. If empty, will use mesh warping instead. 



Number of Chebyshev terms 



Prefix for Chebyshev snes solve 
This problem can be run with the blasius.yaml
file via:
./navierstokes options_file blasius.yaml
problem: 'blasius'
implicit: true
ts:
adapt_type: 'none'
type: 'beuler'
dt: 0.2e5
max_time: 1.0e3
output_freq: 10
## Linear Settings:
degree: 1
dm_plex_box_faces: 40,60,1
platemesh_nDelta: 45
# # Quadratic Settings:
# degree: 2
# dm_plex_box_faces: 20,30,1
# platemesh:
# nDelta: 22
# growth: 1.1664 # 1.08^2
stab: 'supg'
Ctau_t: 1
#Ctau_v: 36,60,128 is what PHASTA has for p=1,2,3
# Linear Settings:
Ctau_v: 36
Ctau_C: 0.25
Ctau_M: 0.25
Ctau_E: 0.125
# # Quadratic Settings:
# Ctau_v: 60
# Ctau_C: 0.125
# Ctau_M: 0.125
# Ctau_E: 0.125
q_extra: 0
dm_plex_box_lower: 0,0,0
dm_plex_box_upper: 4.2e3,4.2e3,5.e5
dm_plex_dim: 3
# Faces labeled 1=z 2=z+ 3=y 4=y+ 5=x+ 6=x
bc_slip_z: 1,2
bc_wall: 3
wall_comps: 1,2,3
bc_inflow: 6
bc_outflow: 5,4
g: 0,0,0
stg:
use: false
inflow_path: "./STGInflow_blasius.dat"
mean_only: true
STG Inflow for Flat Plate¶
Using the STG Inflow for the blasius problem adds the following commandline options:
Option 
Description 
Default value 
Unit 


Path to the STGInflow file 



Path to the STGRand file 



Growth rate of the wavemodes 



Convective velocity, \(U_0\) 



Only impose the mean velocity (no fluctutations) 



Strongly enforce the STG inflow boundary condition 

This problem can be run with the blasius.yaml
file via:
./navierstokes options_file blasius.yaml stg_use true
Note the added stg_use true
flag
This overrides the stg: use: false
setting in the blasius.yaml
file, enabling the use of the STG inflow.
The NavierStokes equations¶
The mathematical formulation (from [GRL10], cf. SE3) is given in what follows. The compressible NavierStokes equations in conservative form are
where \(\bm{\sigma} = \mu(\nabla \bm{u} + (\nabla \bm{u})^T + \lambda (\nabla \cdot \bm{u})\bm{I}_3)\) is the Cauchy (symmetric) stress tensor, with \(\mu\) the dynamic viscosity coefficient, and \(\lambda =  2/3\) the Stokes hypothesis constant. In equations (15), \(\rho\) represents the volume mass density, \(U\) the momentum density (defined as \(\bm{U}=\rho \bm{u}\), where \(\bm{u}\) is the vector velocity field), \(E\) the total energy density (defined as \(E = \rho e\), where \(e\) is the total energy), \(\bm{I}_3\) represents the \(3 \times 3\) identity matrix, \(g\) the gravitational acceleration constant, \(\bm{\hat{k}}\) the unit vector in the \(z\) direction, \(k\) the thermal conductivity constant, \(T\) represents the temperature, and \(P\) the pressure, given by the following equation of state
where \(c_p\) is the specific heat at constant pressure and \(c_v\) is the specific heat at constant volume (that define \(\gamma = c_p / c_v\), the specific heat ratio).
The system (15) can be rewritten in vector form
for the state variables 5dimensional vector
where the flux and the source terms, respectively, are given by
Let the discrete solution be
with \(P=p+1\) the number of nodes in the element \(e\). We use tensorproduct bases \(\psi_{kji} = h_i(X_0)h_j(X_1)h_k(X_2)\).
For the time discretization, we use two types of time stepping schemes.
Explicit timestepping method
The following explicit formulation is solved with the adaptive RungeKuttaFehlberg (RKF45) method by default (any explicit timestepping scheme available in PETSc can be chosen at runtime)
\[ \bm{q}_N^{n+1} = \bm{q}_N^n + \Delta t \sum_{i=1}^{s} b_i k_i \, , \]where
\[ \begin{aligned} k_1 &= f(t^n, \bm{q}_N^n)\\ k_2 &= f(t^n + c_2 \Delta t, \bm{q}_N^n + \Delta t (a_{21} k_1))\\ k_3 &= f(t^n + c_3 \Delta t, \bm{q}_N^n + \Delta t (a_{31} k_1 + a_{32} k_2))\\ \vdots&\\ k_i &= f\left(t^n + c_i \Delta t, \bm{q}_N^n + \Delta t \sum_{j=1}^s a_{ij} k_j \right)\\ \end{aligned} \]and with
\[ f(t^n, \bm{q}_N^n) =  [\nabla \cdot \bm{F}(\bm{q}_N)]^n + [S(\bm{q}_N)]^n \, . \]Implicit timestepping method
This time stepping method which can be selected using the option
implicit
is solved with Backward Differentiation Formula (BDF) method by default (similarly, any implicit timestepping scheme available in PETSc can be chosen at runtime). The implicit formulation solves nonlinear systems for \(\bm q_N\):(19)¶\[ \bm f(\bm q_N) \equiv \bm g(t^{n+1}, \bm{q}_N, \bm{\dot{q}}_N) = 0 \, , \]where the time derivative \(\bm{\dot q}_N\) is defined by
\[ \bm{\dot{q}}_N(\bm q_N) = \alpha \bm q_N + \bm z_N \]in terms of \(\bm z_N\) from prior state and \(\alpha > 0\), both of which depend on the specific time integration scheme (backward difference formulas, generalized alpha, implicit RungeKutta, etc.). Each nonlinear system (19) will correspond to a weak form, as explained below. In determining how difficult a given problem is to solve, we consider the Jacobian of (19),
\[ \frac{\partial \bm f}{\partial \bm q_N} = \frac{\partial \bm g}{\partial \bm q_N} + \alpha \frac{\partial \bm g}{\partial \bm{\dot q}_N}. \]The scalar “shift” \(\alpha\) scales inversely with the time step \(\Delta t\), so small time steps result in the Jacobian being dominated by the second term, which is a sort of “mass matrix”, and typically wellconditioned independent of grid resolution with a simple preconditioner (such as Jacobi). In contrast, the first term dominates for large time steps, with a condition number that grows with the diameter of the domain and polynomial degree of the approximation space. Both terms are significant for timeaccurate simulation and the setup costs of strong preconditioners must be balanced with the convergence rate of Krylov methods using weak preconditioners.
To obtain a finite element discretization, we first multiply the strong form (17) by a test function \(\bm v \in H^1(\Omega)\) and integrate,
with \(\mathcal{V}_p = \{ \bm v(\bm x) \in H^{1}(\Omega_e) \,\, \bm v(\bm x_e(\bm X)) \in P_p(\bm{I}), e=1,\ldots,N_e \}\) a mapped space of polynomials containing at least polynomials of degree \(p\) (with or without the higher mixed terms that appear in tensor product spaces).
Integrating by parts on the divergence term, we arrive at the weak form,
where \(\bm{F}(\bm q_N) \cdot \widehat{\bm{n}}\) is typically replaced with a boundary condition.
Note
The notation \(\nabla \bm v \!:\! \bm F\) represents contraction over both fields and spatial dimensions while a single dot represents contraction in just one, which should be clear from context, e.g., \(\bm v \cdot \bm S\) contracts over fields while \(\bm F \cdot \widehat{\bm n}\) contracts over spatial dimensions.
We solve (20) using a Galerkin discretization (default) or a stabilized method, as is necessary for most realworld flows.
Galerkin methods produce oscillations for transportdominated problems (any time the cell Péclet number is larger than 1), and those tend to blow up for nonlinear problems such as the Euler equations and (lowviscosity/poorly resolved) NavierStokes, in which case stabilization is necessary. Our formulation follows [HST10], which offers a comprehensive review of stabilization and shockcapturing methods for continuous finite element discretization of compressible flows.
SUPG (streamlineupwind/PetrovGalerkin)
In this method, the weighted residual of the strong form (17) is added to the Galerkin formulation (20). The weak form for this method is given as
(21)¶\[ \begin{aligned} \int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t}  \bm{S}(\bm{q}_N) \right) \,dV  \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\ + \int_{\partial \Omega} \bm v \cdot \bm{F}(\bm{q}_N) \cdot \widehat{\bm{n}} \,dS & \\ + \int_{\Omega} \mathcal{P}(\bm v)^T \, \left( \frac{\partial \bm{q}_N}{\partial t} \, + \, \nabla \cdot \bm{F} \, (\bm{q}_N)  \bm{S}(\bm{q}_N) \right) \,dV &= 0 \, , \; \forall \bm v \in \mathcal{V}_p \end{aligned} \]This stabilization technique can be selected using the option
stab supg
.SU (streamlineupwind)
This method is a simplified version of SUPG (21) which is developed for debugging/comparison purposes. The weak form for this method is
(22)¶\[ \begin{aligned} \int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t}  \bm{S}(\bm{q}_N) \right) \,dV  \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\ + \int_{\partial \Omega} \bm v \cdot \bm{F}(\bm{q}_N) \cdot \widehat{\bm{n}} \,dS & \\ + \int_{\Omega} \mathcal{P}(\bm v)^T \, \nabla \cdot \bm{F} \, (\bm{q}_N) \,dV & = 0 \, , \; \forall \bm v \in \mathcal{V}_p \end{aligned} \]This stabilization technique can be selected using the option
stab su
.
In both (22) and (21), \(\mathcal P\) is called the perturbation to the testfunction space, since it modifies the original Galerkin method into SUPG or SU schemes. It is defined as
where parameter \(\bm{\tau} \in \mathbb R^{3}\) (spatial index) or \(\bm \tau \in \mathbb R^{5\times 5}\) (field indices) is an intrinsic time scale matrix. Most generally, we consider \(\bm\tau \in \mathbb R^{3,5,5}\). This expression contains the advective flux Jacobian, which may be thought of as mapping from a 5vector (state) to a \((5,3)\) tensor (flux) or from a \((5,3)\) tensor (gradient of state) to a 5vector (time derivative of state); the latter is used in (23) because it’s applied to \(\nabla\bm v\). The forward variational form can be readily expressed by differentiating \(\bm F_{\text{adv}}\) of (18)
where \(\diff P\) is defined by differentiating (16). This action is also readily computed by forwardmode AD, but since \(\bm v\) is a test function, we actually need the action of the adjoint to use (23) in finite element computation; that can be computed by reversemode AD. We may equivalently write the stabilization term as
where \(\bm r\) is the strong form residual and \(\bm\tau\) is a \(5\times 5\) matrix.
Stabilization scale \(\bm\tau\)
A velocity vector \(\bm u\) can be pulled back to the reference element as \(\bm u_{\bm X} = \nabla_{\bm x}\bm X \cdot \bm u\), with units of reference length (nondimensional) per second. To build intuition, consider a boundary layer element of dimension \((1, \epsilon)\), for which \(\nabla_{\bm x} \bm X = \bigl(\begin{smallmatrix} 2 & \\ & 2/\epsilon \end{smallmatrix}\bigr)\). So a small normal component of velocity will be amplified (by a factor of the aspect ratio \(1/\epsilon\)) in this transformation. The ratio \(\lVert \bm u \rVert / \lVert \bm u_{\bm X} \rVert\) is a covariant measure of (half) the element length in the direction of the velocity. A contravariant measure of element length in the direction of a unit vector \(\hat{\bm n}\) is given by \(\lVert \bigl(\nabla_{\bm X} \bm x\bigr)^T \hat{\bm n} \rVert\). While \(\nabla_{\bm X} \bm x\) is readily computable, its inverse \(\nabla_{\bm x} \bm X\) is needed directly in finite element methods and thus more convenient for our use. If we consider a parallelogram, the covariant measure is larger than the contravariant measure for vectors pointing between acute corners and the opposite holds for vectors between oblique corners.
The cell Péclet number is classically defined by \(\mathrm{Pe}_h = \lVert \bm u \rVert h / (2 \kappa)\) where \(\kappa\) is the diffusivity (units of \(m^2/s\)). This can be generalized to arbitrary grids by defining the local Péclet number
For scalar advectiondiffusion, the stabilization is a scalar
where \(\xi(\mathrm{Pe}) = \coth \mathrm{Pe}  1/\mathrm{Pe}\) approaches 1 at large local Péclet number. Note that \(\tau\) has units of time and, in the transportdominated limit, is proportional to element transit time in the direction of the propagating wave. For advectiondiffusion, \(\bm F(q) = \bm u q\), and thus the perturbed test function is
See [HST10] equations 1517 and 3436 for further discussion of this formulation.
For the NavierStokes and Euler equations, [WJD03] defines a \(5\times 5\) diagonal stabilization \(\mathrm{diag}(\tau_c, \tau_m, \tau_m, \tau_m, \tau_E)\) consisting of
continuity stabilization \(\tau_c\)
momentum stabilization \(\tau_m\)
energy stabilization \(\tau_E\)
The NavierStokes code in this example uses the following formulation for \(\tau_c\), \(\tau_m\), \(\tau_E\):
where \(\bm g = \nabla_{\bm x} \bm{X} \cdot \nabla_{\bm x} \bm{X}\) is the metric tensor and \(\Vert \cdot \Vert_F\) is the Frobenius norm. This formulation is currently not available in the Euler code.
In the Euler code, we follow [HST10] in defining a \(3\times 3\) diagonal stabilization according to spatial criterion 2 (equation 27) as follows.
where \(c_{\tau}\) is a multiplicative constant reported to be optimal at 0.5 for linear elements, \(\hat{\bm n}_i\) is a unit vector in direction \(i\), and \(\nabla_{x_i} = \hat{\bm n}_i \cdot \nabla_{\bm x}\) is the derivative in direction \(i\). The flux Jacobian \(\frac{\partial \bm F_{\text{adv}}}{\partial \bm q} \cdot \hat{\bm n}_i\) in each direction \(i\) is a \(5\times 5\) matrix with spectral radius \((\lambda_{\max \text{abs}})_i\) equal to the fastest wave speed. The complete set of eigenvalues of the Euler flux Jacobian in direction \(i\) are (e.g., [Tor09])
where \(u_i = \bm u \cdot \hat{\bm n}_i\) is the velocity component in direction \(i\) and \(a = \sqrt{\gamma P/\rho}\) is the sound speed for ideal gasses. Note that the first and last eigenvalues represent nonlinear acoustic waves while the middle three are linearly degenerate, carrying a contact wave (temperature) and transverse components of momentum. The fastest wave speed in direction \(i\) is thus
Note that this wave speed is specific to ideal gases as \(\gamma\) is an ideal gas parameter; other equations of state will yield a different acoustic wave speed.
Currently, this demo provides three types of problems/physical models that can be selected at run time via the option problem
.
Advection, the problem of the transport of energy in a uniform vector velocity field, Isentropic Vortex, the exact solution to the Euler equations, and the so called Density Current problem.
Advection¶
A simplified version of system (15), only accounting for the transport of total energy, is given by
with \(\bm{u}\) the vector velocity field. In this particular test case, a blob of total energy (defined by a characteristic radius \(r_c\)) is transported by two different wind types.
Rotation
In this case, a uniform circular velocity field transports the blob of total energy. We have solved (30) applying zero energy density \(E\), and noflux for \(\bm{u}\) on the boundaries.
Translation
In this case, a background wind with a constant rectilinear velocity field, enters the domain and transports the blob of total energy out of the domain.
For the inflow boundary conditions, a prescribed \(E_{wind}\) is applied weakly on the inflow boundaries such that the weak form boundary integral in (20) is defined as
\[ \int_{\partial \Omega_{inflow}} \bm v \cdot \bm{F}(\bm q_N) \cdot \widehat{\bm{n}} \,dS = \int_{\partial \Omega_{inflow}} \bm v \, E_{wind} \, \bm u \cdot \widehat{\bm{n}} \,dS \, , \]For the outflow boundary conditions, we have used the current values of \(E\), following [PMK92] which extends the validity of the weak form of the governing equations to the outflow instead of replacing them with unknown essential or natural boundary conditions. The weak form boundary integral in (20) for outflow boundary conditions is defined as
\[ \int_{\partial \Omega_{outflow}} \bm v \cdot \bm{F}(\bm q_N) \cdot \widehat{\bm{n}} \,dS = \int_{\partial \Omega_{outflow}} \bm v \, E \, \bm u \cdot \widehat{\bm{n}} \,dS \, , \]
Isentropic Vortex¶
Threedimensional Euler equations, which are simplified and nondimensionalized version of system (15) and account only for the convective fluxes, are given by
Following the setup given in [ZZS11], the mean flow for this problem is \(\rho=1\), \(P=1\), \(T=P/\rho= 1\) (Specific Gas Constant, \(R\), is 1), and \(\bm{u}=(u_1,u_2,0)\) while the perturbation \(\delta \bm{u}\), and \(\delta T\) are defined as
where \((\bar{x}, \, \bar{y}) = (xx_c, \, yy_c)\), \((x_c, \, y_c)\) represents the center of the domain, \(r^2=\bar{x}^2 + \bar{y}^2\), and \(\epsilon\) is the vortex strength (\(\epsilon\) < 10). There is no perturbation in the entropy \(S=P/\rho^\gamma\) (\(\delta S=0)\).
Shock Tube¶
This test problem is based on Sod’s Shock Tube (from[sod]), a canonical test case for discontinuity capturing in one dimension. For this problem, the threedimensional Euler equations are formulated exactly as in the Isentropic Vortex problem. The default initial conditions are \(P=1\), \(\rho=1\) for the driver section and \(P=0.1\), \(\rho=0.125\) for the driven section. The initial velocity is zero in both sections. Slip boundary conditions are applied to the side walls and wall boundary conditions are applied at the end walls.
SU upwinding and discontinuity capturing have been implemented into the explicit timestepping operator for this problem. Discontinuity capturing is accomplished using a modified version of the \(YZ\beta\) operator described in [TS07]. This discontinuity capturing scheme involves the introduction of a dissipation term of the form
The shock capturing viscosity is implemented following the first formulation described in {cite} tezduyar2007yzb
. The characteristic velocity \(u_{cha}\) is taken to be the acoustic speed while the reference density \(\rho_{ref}\) is just the local density. Shock capturing viscosity is defined by the following
where,
\(\beta\) is a tuning parameter set between 1 (smoother shocks) and 2 (sharper shocks. The parameter \(h_{SHOCK}\) is a length scale that is proportional to the element length in the direction of the density gradient unit vector. This density gradient unit vector is defined as \(\hat{\bm j} = \frac{\nabla \rho}{\nabla \rho}\). The original formulation of Tezduyar and Senga relies on the shape function gradient to define the element length scale, but this gradient is not available to qFunctions in libCEED. To avoid this problem, \(h_{SHOCK}\) is defined in the current implementation as
where
The constant \(C_{YZB}\) is set to 0.1 for piecewise linear elements in the current implementation. Larger values approaching unity are expected with more robust stabilization and implicit timestepping.
Density Current¶
For this test problem (from [SWW+93]), we solve the full NavierStokes equations (15), for which a cold air bubble (of radius \(r_c\)) drops by convection in a neutrally stratified atmosphere. Its initial condition is defined in terms of the Exner pressure, \(\pi(\bm{x},t)\), and potential temperature, \(\theta(\bm{x},t)\), that relate to the state variables via
where \(P_0\) is the atmospheric pressure. For this problem, we have used noslip and nonpenetration boundary conditions for \(\bm{u}\), and noflux for mass and energy densities.
Channel¶
A compressible channel flow. Analytical solution given in [Whi99]:
where \(H\) is the channel halfheight, \(u_{\max}\) is the center velocity, \(T_w\) is the temperature at the wall, \(Pr=\frac{\mu}{c_p \kappa}\) is the Prandlt number, \(\hat E_c = \frac{u_{\max}^2}{c_p T_w}\) is the modified Eckert number, and \(Re_h = \frac{u_{\max}H}{\nu}\) is the Reynolds number.
Boundary conditions are periodic in the streamwise direction, and noslip and nonpenetration boundary conditions at the walls. The flow is driven by a body force determined analytically from the fluid properties and setup parameters \(H\) and \(u_{\max}\).
Flat Plate Boundary Layer¶
Laminar Boundary Layer  Blasius¶
Simulation of a laminar boundary layer flow, with the inflow being prescribed
by a Blasius similarity
solution. At the inflow,
the velocity is prescribed by the Blasius soution profile, density is set
constant, and temperature is allowed to float. Using weakT: true
, density is
allowed to float and temperature is set constant. At the outlet, a userset
pressure is used for pressure in the inviscid flux terms (all other inviscid
flux terms use interior solution values). The wall is a noslip,
nopenetration, noheat flux condition. The top of the domain is treated as an
outflow and is tilted at a downward angle to ensure that flow is always exiting
it.
Turbulent Boundary Layer¶
Simulating a turbulent boundary layer without modeling the turbulence requires resolving the turbulent flow structures. These structures may be introduced into the simulations either by allowing a laminar boundary layer naturally transition to turbulence, or imposing turbulent structures at the inflow. The latter approach has been taken here, specifically using a synthetic turbulence generation (STG) method.
Synthetic Turbulence Generation (STG) Boundary Condition¶
We use the STG method described in [SSST14]. Below follows a redescription of the formulation to match the present notation, and then a description of the implementation and usage.
Equation Formulation¶
Here, we define the number of wavemodes \(N\), set of random numbers \( \{\bm{\sigma}^n, \bm{d}^n, \phi^n\}_{n=1}^N\), the Cholesky decomposition of the Reynolds stress tensor \(\bm{C}\) (such that \(\bm{R} = \bm{CC}^T\) ), bulk velocity \(U_0\), wavemode amplitude \(q^n\), wavemode frequency \(\kappa^n\), and \(\kappa_{\min} = 0.5 \min_{\bm{x}} (\kappa_e)\).
where \(l_t\) is the turbulence length scale, and \(d_w\) is the distance to the nearest wall.
The set of wavemode frequencies is defined by a geometric distribution:
The wavemode amplitudes \(q^n\) are defined by a model energy spectrum \(E(\kappa)\):
\(\kappa_\eta\) represents turbulent dissipation frequency, and is given as \(2\pi (\nu^3/\varepsilon)^{1/4}\) with \(\nu\) the kinematic viscosity and \(\varepsilon\) the turbulent dissipation. \(\kappa_\mathrm{cut}\) approximates the effective cutoff frequency of the mesh (viewing the mesh as a filter on solution over \(\Omega\)) and is given by:
The enforcement of the boundary condition is identical to the blasius inflow;
it weakly enforces velocity, with the option of weakly enforcing either density
or temperature using the the weakT
flag.
Initialization Data Flow¶
Data flow for initializing function (which creates the context data struct) is given below:
This is done once at runtime. The spatiallyvarying terms are then evaluated at each quadrature point onthefly, either by interpolation (for \(l_t\), \(\varepsilon\), \(C_{ij}\), and \(\overline{\bm u}\)) or by calculation (for \(q^n\)).
The STGInflow.dat
file is a table of values at given distances from the wall.
These values are then interpolated to a physical location (node or quadrature
point). It has the following format:
[Total number of locations] 14
[d_w] [u_1] [u_2] [u_3] [R_11] [R_22] [R_33] [R_12] [R_13] [R_23] [sclr_1] [sclr_2] [l_t] [eps]
where each [ ]
item is a number in scientific notation (ie. 3.1415E0
), and sclr_1
and
sclr_2
are reserved for turbulence modeling variables. They are not used in
this example.
The STGRand.dat
file is the table of the random number set, \(\{\bm{\sigma}^n,
\bm{d}^n, \phi^n\}_{n=1}^N\). It has the format:
[Number of wavemodes] 7
[d_1] [d_2] [d_3] [phi] [sigma_1] [sigma_2] [sigma_3]
The following table is presented to help clarify the dimensionality of the numerous terms in the STG formulation.
Math 
Label 
\(f(\bm{x})\)? 
\(f(n)\)? 

\( \{\bm{\sigma}^n, \bm{d}^n, \phi^n\}_{n=1}^N\) 
RN Set 
No 
Yes 
\(\bm{\overline{u}}\) 
ubar 
Yes 
No 
\(U_0\) 
U0 
No 
No 
\(l_t\) 
l_t 
Yes 
No 
\(\varepsilon\) 
eps 
Yes 
No 
\(\bm{R}\) 
R_ij 
Yes 
No 
\(\bm{C}\) 
C_ij 
Yes 
No 
\(q^n\) 
q^n 
Yes 
Yes 
\(\{\kappa^n\}_{n=1}^N\) 
k^n 
No 
Yes 
\(h_i\) 
h_i 
Yes 
No 
\(d_w\) 
d_w 
Yes 
No 
Meshing¶
The flat plate boundary layer example has custom meshing features to better
resolve the flow. One of those is tilting the top of the domain, allowing for
it to be a outflow boundary condition. The angle of this tilt is controled by
platemesh_top_angle
The primary meshing feature is the ability to grade the mesh, providing better
resolution near the wall. There are two methods to do this; algorithmically, or
specifying the node locations via a file. Algorithmically, a base node
distribution is defined at the inlet (assumed to be \(\min(x)\)) and then
linearly stretched/squeezed to match the slanted top boundary condition. Nodes
are placed such that platemesh_Ndelta
elements are within
platemesh_refine_height
of the wall. They are placed such that the element
height matches a geometric growth ratio defined by platemesh_growth
. The
remaining elements are then distributed from platemesh_refine_height
to the
top of the domain linearly in logarithmic space.
Alternatively, a file may be specified containing the locations of each node.
The file should be newline delimited, with the first line specifying the number
of points and the rest being the locations of the nodes. The node locations
used exactly at the inlet (assumed to be \(\min(x)\)) and linearly
stretched/squeezed to match the slanted top boundary condition. The file is
specified via platemesh_y_node_locs_path
. If this flag is given an empty
string, then the algorithmic approach will be performed.