Lines Matching +full:- +full:a
1 (example-petsc-navier-stokes)=
3 # Compressible Navier-Stokes mini-app
6 …-dependent Navier-Stokes equations of compressible gas dynamics in a static Eulerian three-dimensi…
7 Moreover, the Navier-Stokes example has been developed using PETSc, so that the pointwise physics (…
9 ## Running the mini-app
12 :start-after: <!-- fluids-inclusion -->
14 ## The Navier-Stokes equations
17 The compressible Navier-Stokes equations in conservative form are
22 …abla \cdot \left( \frac{\bm{U} \otimes \bm{U}}{\rho} + P \bm{I}_3 -\bm\sigma \right) - \rho \bm{b}…
23 …t} + \nabla \cdot \left( \frac{(E + P)\bm{U}}{\rho} -\bm{u} \cdot \bm{\sigma} - k \nabla T \right)…
25 $$ (eq-ns)
27 …stress tensor, with $\mu$ the dynamic viscosity coefficient, and $\lambda = - 2/3$ the Stokes hypo…
28 …-ns`, $\rho$ represents the volume mass density, $U$ the momentum density (defined as $\bm{U}=\rho…
31 P = \left( {c_p}/{c_v} -1\right) \left( E - {\bm{U}\cdot\bm{U}}/{(2 \rho)} \right) \, ,
32 $$ (eq-state)
36 The system {eq}`eq-ns` can be rewritten in vector form
39 \frac{\partial \bm{q}}{\partial t} + \nabla \cdot \bm{F}(\bm{q}) -S(\bm{q}) = 0 \, ,
40 $$ (eq-vector-ns)
42 for the state variables 5-dimensional vector
60 - \bm{\sigma} \\
61 - \bm{u} \cdot \bm{\sigma} - k \nabla T
70 $$ (eq-ns-flux)
81 We use tensor-product bases $\psi_{kji} = h_i(X_0)h_j(X_1)h_k(X_2)$.
83 To obtain a finite element discretization, we first multiply the strong form {eq}`eq-vector-ns` by …
86 … \left(\frac{\partial \bm{q}_N}{\partial t} + \nabla \cdot \bm{F}(\bm{q}_N) - \bm{S}(\bm{q}_N) \ri…
89 …}(\Omega_e) \,|\, \bm v(\bm x_e(\bm X)) \in P_p(\bm{I}), e=1,\ldots,N_e \}$ a mapped space of poly…
95 \int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t} - \bm{S}(\bm{q}_N) \right) \…
96 - \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\
100 $$ (eq-weak-vector-ns)
102 where $\bm{F}(\bm q_N) \cdot \widehat{\bm{n}}$ is typically replaced with a boundary condition.
105 …\bm F$ represents contraction over both fields and spatial dimensions while a single dot represent…
111 #### Explicit time-stepping method
113 …formulation is solved with the adaptive Runge-Kutta-Fehlberg (RKF4-5) method by default (any expli…
134 f(t^n, \bm{q}_N^n) = - [\nabla \cdot \bm{F}(\bm{q}_N)]^n + [S(\bm{q}_N)]^n \, .
137 #### Implicit time-stepping method
139 …d using the option `-implicit` is solved with Backward Differentiation Formula (BDF) method by def…
144 $$ (eq-ts-implicit-ns)
152 …e integration scheme (backward difference formulas, generalized alpha, implicit Runge-Kutta, etc.).
153 Each nonlinear system {eq}`eq-ts-implicit-ns` will correspond to a weak form, as explained below.
154 …In determining how difficult a given problem is to solve, we consider the Jacobian of {eq}`eq-ts-i…
160 … by the second term, which is a sort of "mass matrix", and typically well-conditioned independent …
161 …In contrast, the first term dominates for large time steps, with a condition number that grows wit…
162 …Both terms are significant for time-accurate simulation and the setup costs of strong precondition…
167 We solve {eq}`eq-weak-vector-ns` using a Galerkin discretization (default) or a stabilized method, …
169 …-dominated problems (any time the cell Péclet number is larger than 1), and those tend to blow up …
170 …ation follows {cite}`hughesetal2010`, which offers a comprehensive review of stabilization and sho…
172 - **SUPG** (streamline-upwind/Petrov-Galerkin)
174 …ighted residual of the strong form {eq}`eq-vector-ns` is added to the Galerkin formulation {eq}`eq…
179 …\int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t} - \bm{S}(\bm{q}_N) \right) …
180 - \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\
183 \nabla \cdot \bm{F} \, (\bm{q}_N) - \bm{S}(\bm{q}_N) \right) \,dV &= 0
186 $$ (eq-weak-vector-ns-supg)
188 This stabilization technique can be selected using the option `-stab supg`.
190 - **SU** (streamline-upwind)
192 …This method is a simplified version of *SUPG* {eq}`eq-weak-vector-ns-supg` which is developed for …
196 …\int_{\Omega} \bm v \cdot \left( \frac{\partial \bm{q}_N}{\partial t} - \bm{S}(\bm{q}_N) \right) …
197 - \int_{\Omega} \nabla \bm v \!:\! \bm{F}(\bm{q}_N)\,dV & \\
202 $$ (eq-weak-vector-ns-su)
204 This stabilization technique can be selected using the option `-stab su`.
206 In both {eq}`eq-weak-vector-ns-su` and {eq}`eq-weak-vector-ns-supg`, $\bm\tau \in \mathbb R^{5\time…
207 … an ansatz for subgrid state fluctuations $\tilde{\bm q} = -\bm\tau \bm r$ where $\bm r$ is a stro…
208 …riational form can be readily expressed by differentiating $\bm F_{\text{adv}}$ of {eq}`eq-ns-flux`
215 (\diff\bm U \otimes \bm U + \bm U \otimes \diff\bm U)/\rho - (\bm U \otimes \bm U)/\rho^2 \diff\rho…
216 (E + P)\diff\bm U/\rho + (\diff E + \diff P)\bm U/\rho - (E + P) \bm U/\rho^2 \diff\rho
221 where $\diff P$ is defined by differentiating {eq}`eq-state`.
224 A velocity vector $\bm u$ can be pulled back to the reference element as $\bm u_{\bm X} = \nabla_{\…
225 To build intuition, consider a boundary layer element of dimension $(1, \epsilon)$, for which $\nab…
226 So a small normal component of velocity will be amplified (by a factor of the aspect ratio $1/\epsi…
227 The ratio $\lVert \bm u \rVert / \lVert \bm u_{\bm X} \rVert$ is a covariant measure of (half) the …
228 A contravariant measure of element length in the direction of a unit vector $\hat{\bm n}$ is given …
230 If we consider a parallelogram, the covariant measure is larger than the contravariant measure for …
237 $$ (eq-peclet)
239 For scalar advection-diffusion, the stabilization is a scalar
243 $$ (eq-tau-advdiff)
245 where $\xi(\mathrm{Pe}) = \coth \mathrm{Pe} - 1/\mathrm{Pe}$ approaches 1 at large local Péclet num…
246 Note that $\tau$ has units of time and, in the transport-dominated limit, is proportional to elemen…
247 For advection-diffusion, $\bm F(q) = \bm u q$, and thus the SU stabilization term is
251 $$ (eq-su-stabilize-advdiff)
253 where the term in parentheses is a rank-1 diffusivity tensor that has been pulled back to the refer…
254 See {cite}`hughesetal2010` equations 15-17 and 34-36 for further discussion of this formulation.
256 For the Navier-Stokes and Euler equations, {cite}`whiting2003hierarchical` defines a $5\times 5$ di…
261 The Navier-Stokes code in this example uses the following formulation for $\tau_c$, $\tau_m$, $\tau…
281 For Advection-Diffusion, we use a modified version of the formulation for Navier-Stokes:
286 + \frac{\kappa^2 \Vert \bm g \Vert_F ^2}{C_d} \right]^{-1/2}
289 Otherwise, $C_a$ is set via `-Ctau_a` and $C_t$ via `-Ctau_t`.
291 In the Euler code, we follow {cite}`hughesetal2010` in defining a $3\times 3$ diagonal stabilizatio…
295 $$ (eq-tau-conservative)
297 where $c_{\tau}$ is a multiplicative constant reported to be optimal at 0.5 for linear elements, $\…
298 …{\text{adv}}}{\partial \bm q} \cdot \hat{\bm n}_i$ in each direction $i$ is a $5\times 5$ matrix w…
302 \Lambda_i = [u_i - a, u_i, u_i, u_i, u_i+a],
303 $$ (eq-eigval-advdiff)
305 where $u_i = \bm u \cdot \hat{\bm n}_i$ is the velocity component in direction $i$ and $a = \sqrt{\…
306 …ear acoustic waves while the middle three are linearly degenerate, carrying a contact wave (temper…
310 …}} \Bigl( \frac{\partial \bm F_{\text{adv}}}{\partial \bm q} \cdot \hat{\bm n}_i \Bigr) = |u_i| + a
311 $$ (eq-wavespeed)
313 … as $\gamma$ is an ideal gas parameter; other equations of state will yield a different acoustic w…
317 …three types of problems/physical models that can be selected at run time via the option `-problem`.
318 …-advection`, the problem of the transport of energy in a uniform vector velocity field, {ref}`prob…
321 For scale-resolving simulations (such as LES and DNS), statistics for a simulation are more often u…
325 Denote $\langle \phi \rangle$ as the Reynolds average of $\phi$, which in this case would be a aver…
328 \langle \phi \rangle(x,y) = \frac{1}{L_z + (T_f - T_0)}\int_0^{L_z} \int_{T_0}^{T_f} \phi(x, y, z, …
337 The function $\langle \phi \rangle (x,y)$ is represented on a 2-D finite element grid, taken from t…
341 Define a function space on the parent grid as $\mathcal{V}_p^\mathrm{parent} = \{ \bm v(\bm x) \in …
345 To represent these higher-order functions on the parent FEM space, we perform an $L^2$ projection.
353 The projection of a function $u$ onto the parent FEM space would look like:
372 To do this efficiently, **we assume and exploit the full domain grid to be a tensor product in the …
377 To calculate the temporal integral, we do a running average using left-rectangle rule.
378 At the beginning of each simulation, the time integral of a statistic is set to 0, $\overline{\phi}…
379 Periodically, the integral is updated using left-rectangle rule:
383 When stats are written out to file, this running sum is then divided by $T_f - T_0$ to get the time…
388 \bm M [\langle \phi \rangle]_N = \frac{1}{L_z + (T_f - T_0)} \int_\Omega \int_{T_0}^{T_f} \phi(x,y,…
390 where the integral $\int_{T_0}^{T_f} \phi(x,y,z,t) \mathrm{d}t$ is calculated on a running basis.
394 As the simulation runs, it takes a running time average of the statistics at the full domain quadra…
395 This running average is only updated at the interval specified by `-ts_monitor_turbulence_spanstats…
396 …is only solved when statistics are written to file, which is controlled by `-ts_monitor_turbulence…
406 | ----------------- | -------- |
420 To get second-order statistics from these terms, simply use the identity:
423 \langle \phi' \theta' \rangle = \langle \phi \theta \rangle - \langle \phi \rangle \langle \theta \…
426 (differential-filtering)=
434 \overline{\phi} - \nabla \cdot (\beta (\bm{D}\bm{\Delta})^2 \nabla \overline{\phi} ) = \phi
437 …a symmetric positive-definite rank 2 tensor defining the width of the filter, $\bm{D}$ is the filt…
442 - \cancel{\int_{\partial \Omega} \beta v \nabla \overline \phi \cdot (\bm{D}\bm{\Delta})^2 \bm{\hat…
446 The boundary integral resulting from integration-by-parts is crossed out, as we assume that $(\bm{D…
452 It is common to denote a filter width dimensioned relative to the radial distance of the filter ker…
454 For example, under this definition a box filter would be defined as:
466 This is set via `-diff_filter_grid_based_width`.
468 … filter width tensor is most conveniently defined by $\bm{\Delta} = \bm{g}^{-1/2}$ where $\bm g = …
472 The coefficients for that anisotropic scaling are given by `-diff_filter_width_scaling`, denoted he…
490 \zeta = 1 - \exp\left(-\frac{y^+}{A^+}\right)
493 where $y^+$ is the wall-friction scaled wall-distance ($y^+ = y u_\tau / \nu = y/\delta_\nu$), $A^+…
494 … we assume that $\delta_\nu$ is constant across the wall and is defined by `-diff_filter_friction_…
495 $A^+$ is defined by `-diff_filter_damping_constant`.
497 To apply this scalar damping coefficient to the filter width tensor, we construct the wall-damping …
499 The wall-normal filter width is allowed to be damped to a zero filter width.
500 It is currently assumed that the second component of the filter width tensor is in the wall-normal …
513 While we define $\bm{D}\bm{\Delta}$ to be of a certain physical filter width, the actual width of t…
515 To match the "size" of a normal kernel to our differential kernel, we attempt to have them match se…
517 $\beta$ can be set via `-diff_filter_kernel_scaling`.
519 (problem-advection)=
520 ## Advection-Diffusion
522 A simplified version of system {eq}`eq-ns`, only accounting for the transport of total energy, is g…
525 \frac{\partial E}{\partial t} + \nabla \cdot (\bm{u} E ) - \kappa \nabla E = 0 \, ,
526 $$ (eq-advection)
529 In this particular test case, a blob of total energy (defined by a characteristic radius $r_c$) is …
531 - **Rotation**
533 In this case, a uniform circular velocity field transports the blob of total energy.
534 …We have solved {eq}`eq-advection` applying zero energy density $E$, and no-flux for $\bm{u}$ on th…
536 - **Translation**
538 …In this case, a background wind with a constant rectilinear velocity field, enters the domain and …
540 …nditions, a prescribed $E_{wind}$ is applied weakly on the inflow boundaries such that the weak fo…
547 …The weak form boundary integral in {eq}`eq-weak-vector-ns` for outflow boundary conditions is defi…
553 (problem-euler-vortex)=
557 Three-dimensional Euler equations, which are simplified and nondimensionalized version of system {e…
565 $$ (eq-euler)
570 …2 \pi} \, e^{0.5(1-r^2)} \, (-\bar{y}, \, \bar{x}) \, , \\ \delta T &= - \frac{(\gamma-1) \, \epsi…
573 where $(\bar{x}, \, \bar{y}) = (x-x_c, \, y-y_c)$, $(x_c, \, y_c)$ represents the center of the dom…
576 (problem-shock-tube)=
580 …m{cite}`sodshocktubewiki`), a canonical test case for discontinuity capturing in one dimension. Fo…
582 …g a modified version of the $YZ\beta$ operator described in {cite}`tezduyar2007yzb`. This disconti…
600 $\beta$ is a tuning parameter set between 1 (smoother shocks) and 2 (sharper shocks. The parameter …
603 h_{SHOCK} = 2 \left( C_{YZB} \,|\, \bm p \,|\, \right)^{-1}
614 (problem-density-current)=
617 …d to test non-reflecting/Riemann boundary conditions. It's primarily intended for Euler equations,…
619 The problem has a perturbed initial condition and lets it evolve in time. The initial condition con…
623 \rho &= \rho_\infty\left(1+A\exp\left(\frac{-(\bar{x}^2 + \bar{y}^2)}{2\sigma^2}\right)\right) \\
625 E &= \frac{p_\infty}{\gamma -1}\left(1+A\exp\left(\frac{-(\bar{x}^2 + \bar{y}^2)}{2\sigma^2}\right)…
629 where $A$ and $\sigma$ are the amplitude and width of the perturbation, respectively, and $(\bar{x}…
630 The simulation produces a strong acoustic wave and leaves behind a cold thermal bubble that advects…
632 …ing an HLL (Harten, Lax, van Leer) Riemann solver {cite}`toro2009` (option `-freestream_riemann hl…
633 …g a more sophisticated Riemann solver such as HLLC {cite}`toro2009` (option `-freestream_riemann h…
635 ## Vortex Shedding - Flow past Cylinder
637 A cylinder with diameter $D=1$ is centered at $(0,0)$ in a computational domain $-4.5 \leq x \leq 1…
638 We solve this as a 3D problem with (default) one element in the $z$ direction.
640 The viscosity is 0.01 and thermal conductivity is 14.34 to maintain a Prandtl number of 0.71, which…
641 …subjected to freestream boundary conditions at the inflow (left) and Riemann-type outflow on the r…
642 A symmetry (adiabatic free slip) condition is imposed at the top and bottom boundaries $(y = \pm 4.…
643 The cylinder wall is an adiabatic (no heat flux) no-slip boundary condition.
644 As we evolve in time, eddies appear past the cylinder leading to a vortex shedding known as the vor…
663 …, we solve the full Navier-Stokes equations {eq}`eq-ns`, for which a cold air bubble (of radius $r…
667 \begin{aligned} \rho &= \frac{P_0}{( c_p - c_v)\theta(\bm{x},t)} \pi(\bm{x},t)^{\frac{c_v}{ c_p - c…
671 For this problem, we have used no-slip and non-penetration boundary conditions for $\bm{u}$, and no…
675 A compressible channel flow. Analytical solution given in
678 $$ u_1 = u_{\max} \left [ 1 - \left ( \frac{x_2}{H}\right)^2 \right] \quad \quad u_2 = u_3 = 0$$
679 $$T = T_w \left [ 1 + \frac{Pr \hat{E}c}{3} \left \{1 - \left(\frac{x_2}{H}\right)^4 \right \} \ri…
680 $$p = p_0 - \frac{2\rho_0 u_{\max}^2 x_1}{Re_H H}$$
682 where $H$ is the channel half-height, $u_{\max}$ is the center velocity, $T_w$ is the temperature a…
684 Boundary conditions are periodic in the streamwise direction, and no-slip and non-penetration bound…
685 The flow is driven by a body force determined analytically from the fluid properties and setup para…
689 ### Laminar Boundary Layer - Blasius
691 Simulation of a laminar boundary layer flow, with the inflow being prescribed
692 by a [Blasius similarity
696 allowed to float and temperature is set constant. At the outlet, a user-set
698 flux terms use interior solution values). The wall is a no-slip,
699 no-penetration, no-heat flux condition. The top of the domain is treated as an
700 outflow and is tilted at a downward angle to ensure that flow is always exiting
705 Simulating a turbulent boundary layer without modeling the turbulence requires
707 into the simulations either by allowing a laminar boundary layer naturally
709 latter approach has been taken here, specifically using a *synthetic turbulence
715 {cite}`shurSTG2014`. Below follows a re-description of the formulation to match
716 the present notation, and then a description of the implementation and usage.
727 \bm{\hat{x}}^n &= \left[(x - U_0 t)\max(2\kappa_{\min}/\kappa^n, 0.1) , y, z \right]^T
745 The set of wavemode frequencies is defined by a geometric distribution:
748 \kappa^n = \kappa_{\min} (1 + \alpha)^{n-1} \ , \quad \forall n=1, 2, ... , N
751 The wavemode amplitudes $q^n$ are defined by a model energy spectrum $E(\kappa)$:
754 …appa^n}{\sum^N_{n=1} E(\kappa^n)\Delta \kappa^n} \ ,\quad \Delta \kappa^n = \kappa^n - \kappa^{n-1}
760 f_\eta = \exp \left[-(12\kappa /\kappa_\eta)^2 \right], \quad
761 f_\mathrm{cut} = \exp \left( - \left [ \frac{4\max(\kappa-0.9\kappa_\mathrm{cut}, 0)}{\kappa_\mathr…
765 (\nu^3/\varepsilon)^{-1/4}$ with $\nu$ the kinematic viscosity and
767 effective cutoff frequency of the mesh (viewing the mesh as a filter on
776 or temperature using the the `-weakT` flag.
804 lt --Calc-->ke --Calc-->kn
805 y --Calc-->ke
811 u0 --Copy--> u0C[U0]
817 ubar --Copy--> ubarC;
818 y --Copy--> yC;
819 lt --Copy--> ltC;
820 eps --Copy--> epsC;
822 rand --Copy--> randC;
823 rand --> N --Calc--> kn;
824 Rij --Calc--> Cij[C_ij]
827 This is done once at runtime. The spatially-varying terms are then evaluated at
828 each quadrature point on-the-fly, either by interpolation (for $l_t$,
831 The `STGInflow.dat` file is a table of values at given distances from the wall.
832 These values are then interpolated to a physical location (node or quadrature
838 where each `[ ]` item is a number in scientific notation (ie. `3.1415E0`), and `sclr_1` and
853 | ----------------- | -------- | -------------- | --------- |
869 …d by {cite}`shurSTG2014`, but is implemented here as a ramped volumetric forcing term, similar to …
873 S(\bm{q}) = -\sigma(\bm{x})\left.\frac{\partial \bm{q}}{\partial \bm{Y}}\right\rvert_{\bm{q}} \bm{Y…
876 …- P_\mathrm{ref}, \bm{0}, 0]^T$, and $\sigma(\bm{x})$ is a linear ramp starting at `-idl_start` wi…
877 The damping is defined in terms of a pressure-primitive anomaly $\bm Y'$ converted to conservative …
878 $P_\mathrm{ref}$ has a default value equal to `-reference_pressure` flag, with an optional flag `-i…
882 …yer example has custom meshing features to better resolve the flow when using a generated box mesh.
883 …y the nodal layout of the default, equispaced box mesh and are enabled via `-mesh_transform platem…
884 One of those is tilting the top of the domain, allowing for it to be a outflow boundary condition.
885 The angle of this tilt is controlled by `-platemesh_top_angle`.
889 specifying the node locations via a file. Algorithmically, a base node
892 are placed such that `-platemesh_Ndelta` elements are within
893 `-platemesh_refine_height` of the wall. They are placed such that the element
894 height matches a geometric growth ratio defined by `-platemesh_growth`. The
895 remaining elements are then distributed from `-platemesh_refine_height` to the
898 Alternatively, a file may be specified containing the locations of each node.
903 specified via `-platemesh_y_node_locs_path`. If this flag is given an empty
906 ## Taylor-Green Vortex
908 This problem is really just an initial condition, the [Taylor-Green Vortex](https://en.wikipedia.or…
913 v &= -V_0 \cos(\hat x) \sin(\hat y) \sin(\hat z) \\
921 This coordinate modification is done to transform a given grid onto a domain of $x,y,z \in [0, 2\pi…
923 This initial condition is traditionally given for the incompressible Navier-Stokes equations.
924 …e reference state is selected using the `-reference_{velocity,pressure,temperature}` flags (Euclid…