Physics related keywords

Compressible flow

Shock capturing

The shock-capturing module helps stabilize cases with discontinuous solutions, and may also improve the results of under-resolved turbulent cases. It is built on top of a Sensors module that detects problematic flow regions, classifying them according to the value of the sensor, $s$, mapped into the interval $a \in [0,1]$,

$$a = \left\{\begin{array}{ll} 0, & \text{if } s \leq s_0 - \Delta s / 2, \\ \frac{1}{2}\left[1+\sin\left(\frac{s-s_0}{\Delta s}\right)\right], & \text{if } s_0 - \Delta s / 2 < s < s_0 + \Delta s / 2, \\ 1, & \text{elsewhere}. \end{array}\right.$$

The values of $s_0 = (s_1 + s_2)/2$ and$\Delta s = s_2 - s_1$ depend on the sensor thresholds$s_1$ and$s_2$.

At the moment, flow regions where $a \leq 0$ are considered smooth and no stabilization algorithm can be imposed there. In the central region of the sensor, with $0 < a < 1$, the methods shown in the next table can be used and even scaled with the sensor value, so that their intensity increases in elements with more instabilities. Finally, the higher part of the sensor range can implement a different method from the table; however, the intensity is set to the maximum this time.

All the methods implemented introduce artificial dissipation into the equations, which can be filtered with an SVV kernel to reduce the negative impact on the accuracy of the solution. Its intensity is controlled with the parameters$\mu$ (similar to the viscosity of the Navier-Stokes equations) and$\alpha$ (scaling of the density-regularization term of the Guermond-Popov flux), which can be set as constants or coupled to the value of the sensor or to a Smagorinsky formulation.

Spectral vanishing viscosity

The introduction of an SVV-filtered artificial flux helps dissipate high-frequency oscillations. The baseline viscous flux can be chosen as the Navier-Stokes viscous flux or the flux developed by Guermond and Popov. In any case, this flux is expressed in a modal base where it is filtered by any of the following three filter kernels:

• power: $\hat{F}^{\text{1D}}_i = (i/N)^P$,
• sharp: $\hat{F}^{\text{1D}}_i = 0$ if $i<P$, $\hat{F}^{\text{1D}}_i = 1$ elsewhere,
• exponential: $\hat{F}^{\text{1D}}_i = 0$ if $i \leq P$, $\hat{F}^{\text{1D}}_i=\exp\left(-\frac{(i-N)^2}{(i-P)^2}\right)$ elsewhere.

The extension to three dimensions allows the introduction of two types of kernels based on the one-dimensional ones:

• high-pass: $\hat{F}^{\text{H}}_{ijk} = \hat{F}^{\text{1D}}_i \hat{F}^{\text{1D}}_j \hat{F}^{\text{1D}}_k$,
• low-pass: $\hat{F}^{\text{L}}_{ijk} = 1 - \left(1-\hat{F}^{\text{1D}}_i\right)\left(1-\hat{F}^{\text{1D}}_j\right)\left(1-\hat{F}^{\text{1D}}_k\right)$,

being the low-pass one more dissipative and, thus, more suited to supersonic cases. The high-pass filter, on the other hand, works better as part of the SVV-LES framework for turbulent cases.

The cutoff parameter $P$ can be set as "automatic", which uses a sensor to differentiate troubled elements from smooth regions. The stabilisation strategy then depends on the region:

• smooth regions: $P=4$, $\mu=\mu_2$, $\alpha=\alpha_2$,
• shocks: $P=4$, $\mu=\mu_1$, $\alpha=\alpha_1$.

In addition to this, the viscosity $\mu_1$ can be set to "Smagorinsky" to use the implemented SVV-LES approach. In this case, the $\mu=\mu_{\text{LES}}$ viscosity is computed following a Smagorinsky formulation with $C_s=\mu_2$ and the viscosity parameters do not depend on the region anymore,

$$\mu = C_s^2 \Delta^2|S|^2, \quad \alpha = \alpha_1.$$

Acoustic

The Ffowcs Williams and Hawkings (FWH) acoustic analogy is implemented as a complement to the compressible NS solver. It can run both during the execution of the NS (in-time) or as at a post-process step. The version implemented includes both the solid and permeable surface variations, but both of them for a static body subjected to a constant external flow, i.e. a wind tunnel case scenario. The specifications for the FWH are divided in two parts: the general definitions (including the surfaces) and the observers definitions. The former is detailed in the table below, while the latter are defined in a block section, similar to the monitors (see monitors):

#define acoustic observer 1
   name     = SomeName
   position = [0.d0, 0.d0, 0.d0]
#end
# end

To run the in-time computation, the observers must be defined in the control file. Beware that adding an additional observer will require to run the simulation again. To use the post-process computation, the solution on the surface must be saved at a regular time. Beware that it will need more storage. To run the post-process calculation the horses.tools binary is used, with a control file similar to the one use for the NS simulation (without monitors), and adding the keywords ''tool type'' and ''acoustic files pattern'', as explained in the table below:

Incopressible navier-Stokes

Among the various incompressible Navier-Stokes models, HORSES3D uses an artificial compressibility formulation, which converts the elliptic problem into a hyperbolic system of equations, at the expense of a non divergence–free velocity field. However, it allows one to avoid constructing an approximation that satisfies the inf–sup condition. This methodology is well suited for use as a fluid flow engine for interface–tracking multiphase flow models, as it allows the density to vary spatially.

The artificial compressibility system of equations is:

$$\begin{equation} \rho_{t} + \vec{\nabla}\cdot\left(\rho\vec{u}\right)=0 , \label{eq:incomressible:continuity} \end{equation}$$

$$\begin{equation} \left(\rho\vec{u}\right)_t+\vec{\nabla}\cdot\left(\rho \vec{u}\vec{u}\right) = -\vec{\nabla}p + \vec{\nabla}\cdot\left(\frac{1}{\mathrm{Re}}\left(\vec{\nabla}\vec{u} + \vec{\nabla}\vec{u}^{T}\right)\right)+\frac{1}{\mathrm{Fr^{2}}}\rho\vec{e}_{g}, \label{eq:incomressible:momentum} \end{equation}$$

$$\begin{equation} p_t + \rho_0 c_0^2 \vec{\nabla}\cdot\vec{u} = 0, \label{eq:incomressible:ACM} \end{equation}$$

The factor $\rho_0$ is computed as $max(\rho_1/\rho_1,\rho_2/\rho_1)$.

This solver is run with the binary horses3d.ins.

The incompressible Navier Stokes solver has two modes: with 1 fluid and with 2 fluids. The required keywords are listed below.

Table 1: Keywords for incompressible Navier-Stokes solver. Mode with 1 fluid

Table 2: Keywords for incompressible Navier-Stokes solver. Mode with 2 fluids

Multiphase

The multiphase flow solver implemented in HORSES3D is constructed by a combination of the diffuse interface model of Cahn–Hilliard with the incompressible Navier–Stokes equations with variable density and artificial compressibility. This model is entropy stable and guarantees phase conservation with an accurate representation of surface tension effects. The modified entropy-stable version approximates:

$$\begin{equation} c_t + \vec{\nabla}\cdot\left(c\vec{u}\right) = M_0 \vec{\nabla}^2 \mu, \label{eq:governing:cahn--hilliard} \end{equation}$$ % $$\begin{equation} \sqrt{\rho}\left(\sqrt{\rho}\vec{u}\right)_t+\vec{\nabla}\cdot\left(\frac{1}{2}\rho \vec{u}\vec{u}\right) +\frac{1}{2}\rho\vec{u}\cdot\vec{\nabla}\vec{u}+c\vec{\nabla}\mu = -\vec{\nabla}p + \vec{\nabla}\cdot\left(\eta\left(\vec{\nabla}\vec{u} + \vec{\nabla}\vec{u}^{T}\right)\right)+\rho\vec{g}, \label{eq:governing:momentum-skewsymmetric-sqrtRho} \end{equation}$$ % $$\begin{equation} p_t + \rho_0 c_0^2 \vec{\nabla}\cdot\vec{u} = 0, \label{eq:governing:ACM} \end{equation}$$ % where $c$ is the phase field parameter, $M_0$ is the mobility, $\mu$ is the chemical potential, $\eta$ is the viscosity and $c_0$ is the artificial speed of sound. The factor $\rho_0$ is computed as $max(\rho_1/\rho_1,\rho_2/\rho_1)$. Mobility $M_0$ is computed from the control file parameters chemical characteristic time $t_{CH}$, interface width $\epsilon$ and interface tension $\sigma$ with the formula $M_0 = L_{ref}^2 \epsilon /(t_{CH} \sigma)$.

The term $M_0 \vec{\nabla}^2 \mu$ can be implicity integrated to reduce the stiffnes of the problem with the keyword time integration = IMEX. This is only recomended if the value of $M_0$ is very high so that the time step of the explicit scheme is very small.

This solver is run with the binary horses3d.mu.

Particles

Horses3d includes a two-way coupled Lagrangian solver. Particles are tracked along their trajectories, according to the simplified particle equation of motion, where only contributions from Stokes drag and gravity are retained, $$\begin{equation} \label{eq:part_motion} \frac{d y_i}{dt} = u_i, \quad \frac{d u_i}{dt} = \frac{v_i - u_i}{\tau_p} + g_i, \end{equation}$$ where $u_i$ and $y_i$ are the \emph{ith} components of velocity and position of the particle, respectively. Furthermore, $v_i$ accounts for the continuous velocity of the fluid at the position of the particle. We consider spherical Stokesian particles, so their mass and aerodynamic response time are $m_p = \rho_p \pi D_p^3/6$ and $\tau_p = \rho_p D_p^2 / 18\mu$, respectively, $\rho_p$ being the particle density and $D_p$ the particle diameter.

Each particle is considered to be subject to a radiative heat flux $I_o$. The carrier phase is transparent to radiation, whereas the incident radiative flux on each particle is completely absorbed. Because we focus on relatively small volume fractions, the fluid-particle medium is considered to be optically thin. Under these hypotheses, the direction of the radiation is inconsequential, and each particle receives the same radiative heat flux, and its temperature $T_p$ is governed by $$\begin{equation} \label{eq:part_energy} \frac{d}{dt} (m_p c_{V,p} T_p) = \frac{\pi D_p^2}{4} I_o - \pi D_p^2 h (T_p-T), \end{equation}$$ where $c_{V,p}$ is the specific heat of the particle, which is assumed to be constant with respect to temperature. $T_p$ is the particle temperature and $h$ is the convective heat transfer coefficient, which for a Stokesian particle can be calculated from the Nusselt number $Nu = hD_p/k = 2$.

In practical simulations, integrating the trajectory of every particle is too expensive. Therefore, particles are agglomerated into parcels, each of them accounting for many particles with the same physical properties, position, velocity, and temperature. The evolution of the parcels is tracked with the same set of equations presented for the particles.

The two-way coupling means that fluid flow is modified because of the presence of particles. Therefore, the Navier-Stokes equations are enriched with the following source terms: % $$\begin{equation} \boldsymbol{S} = \beta\left[\begin{array}{c} 0 \\ \sum_{n=1}^{N_p} \frac{m_p}{\tau_p} (u_{1,n}-v_1)\delta(\mathbf{x} - \mathbf{y}_n) \\ \sum_{n=1}^{N_p} \frac{m_p}{\tau_p} (u_{2,n}-v_2)\delta(\mathbf{x} - \mathbf{y}_n) \\ \sum_{n=1}^{N_p} \frac{m_p}{\tau_p} (u_{3,n}-v_3)\delta(\mathbf{x} - \mathbf{y}_n) \\ \sum_{n=1}^{N_p} \pi D_p^2 h (T_{p,n} - T) \delta(\mathbf{x} - \mathbf{y}_n ) \end{array}\right], \label{eq:particlessource} \end{equation}$$ % where $\delta$ is the Dirac delta function, $N_p$ is the number of parcels, $\beta$ is the number of particles per parcel and $u_{i,n}$, $\mathbf{y}_{i,n}$, $T_{p,n}$ are the velocity, spatial coordinates, and temperature of the parcel \emph{nth}. The dimensionless form of the Navier Stokes equations can be seen in the appendix at the end of this document.

Particles are solved in a box domain inside the flow domain. The box is defined with the keywords minimum box'' andmaximum box''. The boundary conditions for the particles are defined with the keyword ``bc box''. Possible options are inflow/outflow, periodic and wall (with perfect rebound).

Complementary Modes

Wall functions

The wall function overwrites the viscous flux on the specified boundaries based on an specific law using a Newman condition. It must be used as a complement of no slip boundary condition. The table below shows the parameters that can be set in the control file. The frictional velocity is calculated using the instantaneous values of the first node (either Gauss or Gauss-Lobatto) of the element neighbour of the face element (at the opposite side of the boundary face). Currently is only supported for the compressible Navier-Stokes solver.

The standard wall function uses the Reichardt law, solving the algebraic non-linear equation using the newton method to obtain the frictional velocity. The ABL function uses the logarithmic atmospheric boundary layer law, using the aerodynamic roughness; the frictional velocity is without using any numerical method.

Tripping

A numerical source term is added to the momentum equations to replicate the effect of a tripping mechanism used commonly in explerimental tests. The forcing is described via the product of two independent functions: one that depends streamwise and vertical directions (space only) and the other one describing the spanwise direction and time (space and time). It can be used for the compressible NS, both LES and RANS. The keywords for the trip options are: