Presets/cfd

Introduction

cfd.xml is the main preset to be loaded in order to perform Fluid Dynamics simulations.

Governing equations

The following governing equations are approximately solved with this preset (called conservation of mass equation, momentum equation, and equation of state):

${\displaystyle{\frac{d\rho}{dt}}=-\rho{\mathrm{div}}\left({\mathbf{u}}\right),}$

${\displaystyle{\frac{d{\mathbf{u}}}{dt}}={\mathbf{g}}+{\frac{{\mathrm{div}}\left({\mathbb{T}}\right)}{\rho}},}$

$p=p\left(\rho\right);$

where ${\displaystyle\rho}$ and ${\displaystyle{\mathbf{u}}}$ are the density and the velocity respectively, ${\displaystyle{\mathbf{g}}}$ is an external volumetric force field, $p$ is the pressure, and ${\displaystyle{\mathbb{T}}}$ is the stress tensor for a Newtonian fluid:

${\displaystyle{\mathbb{T}}=\left(-p+\lambda{\mathrm{tr}}\left({\mathbb{D}}\right)\right){\mathbb{1}}+2\mu{\mathbb{D}},}$

with ${\displaystyle\lambda}$ and ${\displaystyle\mu}$ the viscous coefficients and ${\displaystyle{\mathbb{D}}}$ the rate of strain tensor:

${\displaystyle{\mathbb{D}}={\frac{\nabla{\mathbf{u}}+\nabla{\mathbf{u}}^{T}}{2}}.}$

As it has been stated in the 3^rd governing equation (formerly Equation Of State, or EOS), we need to stablish a relation between the density and the pressure:

$p={\frac{c_{s}^{2}\rho_{0}}{\gamma}}\left(\left({\frac{\rho}{\rho_{0}}}\right)^{\gamma}+1\right),$

where $c_{s}$ is the sound speed (at least 10 times the maximum speed expected during the simulation), ${\displaystyle\rho_{0}}$ the density of reference of the fluid, and ${\displaystyle\gamma}$ an incompressibility factor (1 is recommended, but it can be increased to get better incompressibility properties).

Numerical modelling

Conservation of mass equation

As usual in SPH the fluid domain is discretized in a set of particles, such that the aproximate solution of the conservation of mass equation can be written as:

${\displaystyle\left\langle{\frac{d\rho_{a}}{dt}}\right\rangle={\frac{1}{\gamma_{a}}}\left[\sum_{{b\in\Omega}}\left({\mathbf{u}}_{b}-{\mathbf{u}}_{a}\right)\cdot\nabla W\left({\mathbf{r}}_{b}-{\mathbf{r}}_{a}\right)m_{b}+\sum_{{b\in\partial\Omega}}\rho_{b}\left({\mathbf{u}}_{b}-{\mathbf{u}}_{a}\right)\cdot{\mathbf{n}}_{b}W\left({\mathbf{r}}_{b}-{\mathbf{r}}_{a}\right)S_{b}\right],}$

where $a$ and $b$ are particle indexes, ${\displaystyle{\mathbf{r}}_{a}}$ and ${\displaystyle{\mathbf{r}}_{b}}$ their positions, $m_{b}$ the mass of the particle $b$ , ${\displaystyle{\mathbf{n}}_{b}}$ and $S_{b}$ the outward normal and the surface of the boundary element $b$ , and ${\displaystyle\gamma_{b}}$ a renormalization factor:

${\displaystyle\gamma_{a}=\sum_{{b\in\Omega}}W\left({\mathbf{r}}_{b}-{\mathbf{r}}_{a}\right){\frac{m_{b}}{\rho_{b}}},}$

Momentum equation

The momentum equation aproximation can be separated in 2 differential operators aproximation:

${\displaystyle\left\langle{\frac{d{\mathbf{u}}_{a}}{dt}}\right\rangle=-\left\langle{\frac{\nabla p_{a}}{\rho_{a}}}\right\rangle+\mu\left\langle\Delta{\mathbf{u}}_{a}\right\rangle+{\mathbf{g}},}$

such that:

${\displaystyle\left\langle{\frac{\nabla p_{a}}{\rho_{a}}}\right\rangle={\frac{1}{\gamma_{a}}}\left[\sum_{{b\in\Omega}}\left({\frac{p_{a}}{\rho_{a}^{2}}}+{\frac{p_{b}}{\rho_{b}^{2}}}\right)\cdot\nabla W\left({\mathbf{r}}_{b}-{\mathbf{r}}_{a}\right)m_{b}+\sum_{{b\in\partial\Omega}}\rho_{b}\left({\frac{p_{a}}{\rho_{a}^{2}}}+{\frac{p_{b}}{\rho_{b}^{2}}}\right)\cdot{\mathbf{n}}_{b}W\left({\mathbf{r}}_{b}-{\mathbf{r}}_{a}\right)S_{b}\right],}$

${\displaystyle\left\langle\Delta{\mathbf{u}}_{a}\right\rangle={\frac{1}{\gamma_{a}}}\sum_{{b\in\Omega}}K{\frac{\left({\mathbf{u}}_{b}-{\mathbf{u}}_{a}\right)\cdot\left({\mathbf{r}}_{b}-{\mathbf{r}}_{a}\right)}{\rho_{a}\rho_{b}\left|{\mathbf{r}}_{b}-{\mathbf{r}}_{a}\right|^{2}}}\nabla W\left({\mathbf{r}}_{b}-{\mathbf{r}}_{a}\right)m_{b};}$

where the constant $K=6,8,10$ for 1D, 2D and 3D respectively.

Time integration

A pseudo second-order time integration scheme is used in AQUAgpusph, composed by a predictor and a corrector stage.

The predictor, executed at the start of the time step, can be written as follows:

${\displaystyle\rho_{{n+1/2}}=\rho_{{n}}+\Delta t_{{n}}\,{\frac{d\rho_{{n}}}{dt}}}$

${\displaystyle{\mathbf{u}}_{{n+1/2}}={\mathbf{u}}_{{n}}+\Delta t_{{n}}\,{\frac{d{\mathbf{u}}_{{n}}}{dt}}}$

${\displaystyle{\mathbf{r}}_{{n+1/2}}={\mathbf{r}}_{{n}}+\Delta t_{{n}}\,{\mathbf{u}}_{{n}}+{\frac{\Delta t_{{n}}^{2}}{2}}\,{\frac{d{\mathbf{u}}_{{n}}}{dt}}}$

Regarding the corrector stage, computed after the derivatives evaluation described before, can be written as:

${\displaystyle\rho_{{n+1}}=\rho_{{n+1/2}}+{\frac{\Delta t_{{n+1}}^{2}}{2}}\left({\frac{d\rho_{{n+1}}}{dt}}-{\frac{d\rho_{{n}}}{dt}}\right)}$

${\displaystyle{\mathbf{u}}_{{n+1}}={\mathbf{u}}_{{n+1/2}}+{\frac{\Delta t_{{n+1}}^{2}}{2}}\left({\frac{d{\mathbf{u}}_{{n+1}}}{dt}}-{\frac{d{\mathbf{u}}_{{n}}}{dt}}\right)}$

Both, predictor and corrector are just applied to the fluid particles.

Usage

To use it just use the following include tag:

<Include file="/usr/share/aquagpusph/resources/Presets/cfd.xml" />

Where /usr/share/aquagpusph/ should be conveniently modified. For backguard compatibility you can include /usr/share/aquagpusph/resources/OpenCLMain.xml instead, which is just redirecting to cfd.xml.

Remember that this preset should be the first included.