FlowDeck: Our Method

Overview

FlowDeck is a computational fluid dynamics tool (CFD). Every simulation solves the equations of fluid motion, and incorporates additional physics where needed (heat transport, species concentration, turbulence equations) for your geometry and boundary conditions.

The flow is solved using a coupled Lattice Boltzmann and finite difference solver, using the same cartesian grid. The methodology and assumptions employed are outlined below. Each product feature is validated with reference data.

Flow Solver: Lattice Boltzmann

The flow field is computed using the Lattice Boltzmann method. Instead of solving the Navier–Stokes equations directly, it tracks how fluid particles move and interact on a grid. At each time step, particles collide and move to neighbouring grid points, allowing velocity, pressure and stress to be calculated. Through this process, the Lattice Boltzmann Method reproduces the Navier-Stokes equations at a macroscopic scale, and is well suited to model unsteady fluid flow, complex geometry, and use parallel computing (GPU).

Detail

The update equation evolves particle distribution functions $f_i$, which stream between neighbouring grid points along a fixed set of velocities $\mathbf{c}_i$ and then "collide" via $\Omega_i$:

$$f_i(\mathbf{x} + \mathbf{c}_i \Delta t,\; t + \Delta t) = f_i(\mathbf{x}, t) + \Omega_i(\mathbf{x}, t) + F_i(\mathbf{x}, t)$$

Pressure, velocity and stress are functions of the distribution function $f_i$.

In the low-Mach number limit, a multi-scale expansion of this update equation recovers the incompressible Navier–Stokes equations, given an appropriate collision function. The update is local and explicit, and formally second order in space and time. More details about the LBM method can be found here [1].

Grid and Level Transfer

The simulation grid is a cartesian cubic grid in 3 dimensions. Grid cell widths are refined by a factor of 2 (halved) near areas of interest. Grid time step scales with the cell widths each refined level.

Neighbouring grid levels overlap via a band of transfer ghost grid nodes. Information flowing from one side of the interface is used to populate these ghost nodes and ensure the smooth transfer of momentum, stress and pressure values across the interface. This is second order accurate in space.

Flow Boundary Conditions

Boundary conditions are applied on the boundary of every simulation domain. These conditions ensure the behaviour of the flow simulation is well-defined. A boundary condition value or profile is able to vary in space and time.

Tranport Solver: Finite Difference

Additional physics such as temperature, turbulence transport and pollutant concentration can be added to the flow field by solving a corresponding advection-diffusion (scalar transport) equation. At each step in the simulation, the flow and transport solvers are coupled. The scalar variable $c$ is advected using the flow velocity, and the flow may experiences an effect from the scalar $c$.

A transported scalar $c$ obeys the advection–diffusion equation:

$$\frac{\partial c}{\partial t} + \mathbf{u}\cdot\nabla c = \nabla\cdot(D\,\nabla c) + Q$$

with $D$ the diffusivity and $Q$ the source terms of the physics being modelled.

Discretisation

The equation is solved using a Finite Difference scheme on the same grid as the flow, and coupled to ensure correct feedback between the solvers. The update is explicit, with the same time step as the flow. Due to this, there is a maximum lattice diffisivity (Courant number) which enforces the solver stability.

• Advection uses the Lax–Wendroff MacCormack scheme, second order in space and time [6].
• Diffusion and any Source terms use second-order central differences.

Boundary conditions

Coupled Solver Features

Turbulence Modelling

FlowDeck can model turbulence in three ways: with a Reynolds-averaged (RANS) model, by resolving it directly down to the grid scale, or by combining the two in a hybrid DES model. The treatment is chosen to suit each simulation.

Reynolds-Averaged: Spalart–Allmaras

Reynolds-averaged (RANS) models approximate the time averaged influence of turbulent structures in the flow. The selected RANS model is Spalart Allmaras (SA), which solves a single equation transport model for the variable $\tilde\nu$:

$$\frac{\partial \tilde\nu}{\partial t} + \mathbf{u}\cdot\nabla\tilde\nu = RHS$$

This is solved via the finite difference solver scheme on the same grid as the flow (see x). The variable $\tilde\nu$ is used to calculate an "eddy viscosity" term $\nu_t$, which represents the additional momentum mixing caused by the turbulent structures. The total viscosity experienced locally at each grid cell by the fluid is now:

$$\nu_{total} = \nu + \nu_t$$

The SA model is computationally efficient and commonly used for aerodynamic flows, particularly turbulent wall boundary layers. The particular model used is the Negative SA model which increases numerical stability in the turbulence solver. The full equations and coefficients can be found here [8].

Scale-resolving: The Cumulant Operator

Time-resolved simulations resolve the larger turbulent structures on the grid and use no turbulence model. Turbulent structures smaller than the grid are removed, hence increasingly small grids are needed as the strutures become small. Conventional large-eddy simulations can approximate these 'sub-grid' structures by adding an explicit sub-grid model, typically an eddy viscosity model such as Smagorinsky:

$$\nu_t = (C_s dx)^2\,|S|$$

Where $S$ is the local momentum shear of the flow, $C_s$ a constant. These types of models add dissipation to any resolved structures and require highly accurate shear calculations to reduce this error, which is computationally expensive.

Instead, the optimised cumulant collision operator is used to model the turbulence [3]. The method cannot resolve any sub-grid structures, but instead reduces the diffusion experienced by the well-resolved structures. The optimised parameter set for the cumulant collision operator enables fourth order diffusion calculations, adding very little damping as would be otherwise experienced with LES.

Individual cumulants of the distribution function $f_i$ may grow large if the flow is under-resolved by the grid. Limits on these cumulants are added as part of this cumulant collision operator scheme, providing dissapation only where cumulants grow large and would otherwise be unstable. The operator therefore acts as its own sub-grid model, without adding any additional local viscosity.

Detached Eddy Simulation

Turbulent flows can be separated into two main regions: the turbulent boundary layer near walls, and the turbulent wake regions away from the wall. To accurately model the flow behaviour in these two regions, a Hybrid DES approach is used.

Hybrid DES is employed via a blend function [7] to apply a Reynolds Averaged (RANS) model in the boundary layer, near wall regions, and then to resolve the larger turbulent structures in the turbulent wake using the accurate cumulant collision method in the Lattice Boltzmann code.

The blend function is necessary as time-averaged RANS models can severely underpredict the unsteady vortex shedding that occurs in turbulent wakes. Turbulent eddy structures are also extremely small in boundary layers, meaning resolving these explicitly using a very small grid is computationally infeasible.

Wall functions

Within a turbulent boundary layer, flow velocity changes rapidly near the wall. To resolve such extremely high shear, grid cells much smaller than reasonable are normally required. Instead, flow in this boundary layer region can be estimated by well established descriptions of this velocity profile (wall functions) without the need for explicit modelling.

From the wall function, velocity, stress and eddy viscosity can be calculated near the wall and applied as a velocity boundary condition and dirichlet boundary condition respectively. This allows the flow to experience the influence of the near wall region of the turbulent boundary, without needing it to be explicitly modelled.

Detail

The wall functions used in the code are equilibrium functions, which do not consider the effect of pressure gradients or the mixing of the flow in their calculation of the velocity profile. Due to this, away from the wall, the wall functions predict a region of constant shear which extends indefinitely named the "log layer". The near wall region modelled by the wall function includes 1 to 2 grid cells from near the wall, and will assume these cells are within this "log-layer" region.

In reality, pressure gradients and flow advection become significant away from the wall and closer to the edge of the boundary layer. In very thin boundary layer simulations, the nearest grid cells from the wall may be too far away to be considered in the log layer, and for the wall function to be accurate. This may also be the case in boundary layer flows experiencing pressure gradients, where the shape of the desired velocity profile is different to the equilibrium wall function.

Equilibrium wall functions are chosen as they computationally efficient, and accurate for a large range of simulation scenarios. The employed wall functions are:

• "Spalart Allmaras": used for smooth wall boundary layers [8].
• "Rough Wall Logarithmic": used for rough wall boundary layers, with a roughness height $z_0$ (eg: Atmospheric Boundary Layer) [9].

Heat

To include heat effects to the flow, the scalar transport equation solves fir the temperature value $T$. The diffusivity $D$ is the thermal fluid diffusivity, and $Q$ is any variable heat source.

$$\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T = \nabla\cdot(D\,\nabla T) + Q$$

The flow is coupled to the temperature by the Boussinesq approximation, creating a bouyancy term to the flow momentum. Fluid is treated as incompressible, where density variations from temperature differences drive the bouyancy force. These density variations are assumed to be small.

$$\mathbf{F} = -\,\mathbf{g}\,\beta\,(T - T_0)$$

$$\rho \approx \rho_0 [1 - \beta(T-T_0)]$$

where $\mathbf{g}$ is gravity unit vector, $T_0$ is a reference temperature, $\beta$ is the thermal expansion coefficient, $\rho$ is the fluid density and $\rho_0$ is the reference density.

Convergence and outputs

Each simulation is performed using explicit time stepping. Each monitored time step, the simulation is checked if convergence to a steady state solution within a tolerance has occured, or if the input number of time steps has been reached.

This convergence tolerance is simply measured as the difference in global velocity magnitudes compared to the previous time step over $N$ cells. For unsteady or time periodic simulations this will never satisfy this criteria.

$$Convergence = \frac{1}{N} \sum_{\text{cells}} \Big|\, |\mathbf{u}^{n+1}| - |\mathbf{u}^{n}| \,\Big| \lt tolerance $$

Outputs

The outputs of the solver are calculated at each grid cell, at the last run timestep. Queried values are calculated using second order interpolation.

• Velocity components ($u_i$)
• Pressure ($p$)
• Eddy Viscosity ratio (${\nu_t}/{\nu}$)
• Scalar Transport variable (eg: $c$, $T$, $\tilde\nu$)

References

1. He, Xiaoyi, and Li-Shi Luo. “Lattice Boltzmann Model for the Incompressible Navier–Stokes Equation.” Journal of Statistical Physics 88, no. 3/4 (1997): 927–44. https://doi.org/10.1023/B:JOSS.0000015179.12689.e4.
2. Geier, Martin, Martin Schönherr, Andrea Pasquali, and Manfred Krafczyk. “The Cumulant Lattice Boltzmann Equation in Three Dimensions: Theory and Validation.” Computers & Mathematics with Applications 70, no. 4 (2015): 507–47. https://doi.org/10.1016/j.camwa.2015.05.001.
3. Geier, Martin, Andrea Pasquali, and Martin Schönherr. “Parametrization of the Cumulant Lattice Boltzmann Method for Fourth Order Accurate Diffusion Part I: Derivation and Validation.” Journal of Computational Physics 348 (November 2017): 862–88. https://doi.org/10.1016/j.jcp.2017.05.040.
4. Ginzburg, Irina, and Frederik Verhaeghe. “Two-Relaxation-Time Lattice Boltzmann Scheme: About Parametrization, Velocity, Pressure and Mixed Boundary Conditions.” Commun. Comput. Phys., 2008.
5. Junk, Michael, and Zhaoxia Yang. “Outflow Boundary Conditions for the Lattice Boltzmann Method.” Progress in Computational Fluid Dynamics, An International Journal 8, no. 1/2/3/4 (2008): 38. https://doi.org/10.1504/PCFD.2008.018077.
6. Machalinska-Murawska, Justyna, and Michał Szydłowski. “Lax-Wendroff and McCormack Schemes for Numerical Simulation of Unsteady Gradually and Rapidly Varied Open Channel Flow.” Archives of Hydro-Engineering and Environmental Mechanics 60, nos. 1–4 (2014): 51–62. https://doi.org/10.2478/heem-2013-0008.
7. Spalart, P. R., S. Deck, M. L. Shur, K. D. Squires, M. Kh. Strelets, and A. Travin. “A New Version of Detached-Eddy Simulation, Resistant to Ambiguous Grid Densities.” Theoretical and Computational Fluid Dynamics 20, no. 3 (2006): 181–95. https://doi.org/10.1007/s00162-006-0015-0.
8. Allmaras, Steven R., Forrester T. Johnson, and Philippe R. Spalart. “Modifi cations and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model.” Computational Fluid Dynamics, 2012.
9. Schlichting & Gersten (2000) p. 530