Resulting schemes

We discretize the mentioned equations as follows. Note that we define several symbols for notational convenience.

The advective terms are implemented in tri-diagonal-like ways to insist the advective operators are skew-symmetric (c.f., Verstappen and Veldman, J. Comput. Phys. (187), 2003). Since the diffusive terms are given by the vector Laplacian of the velocity vector, the discrete Laplace operators are pre-computed.

Incompressibility constraint

We define this relation for each cell center \(\left( \ccidx{i}, \ccidx{j}, \ccidx{k} \right)\):

\[\ddiv{1} + \ddiv{2} + \ddiv{3} = 0.\]

The left-hand side is monitored to confirm if the maximum divergence of the flow field is sufficiently small (\(\approx 0\)):

src/logging/max_divergence.c
69  double divmax = 0.;
70  BEGIN
71    const double hx_xm = HXXF(i  );
72    const double hx_xp = HXXF(i+1);
73    const double jd_xm = JDXF(i  );
74    const double jd_x0 = JDXC(i  );
75    const double jd_xp = JDXF(i+1);
76#if NDIMS == 2
77    const double ux_xm = UX(i  , j  );
78    const double ux_xp = UX(i+1, j  );
79    const double uy_ym = UY(i  , j  );
80    const double uy_yp = UY(i  , j+1);
81#else
82    const double ux_xm = UX(i  , j  , k  );
83    const double ux_xp = UX(i+1, j  , k  );
84    const double uy_ym = UY(i  , j  , k  );
85    const double uy_yp = UY(i  , j+1, k  );
86    const double uz_zm = UZ(i  , j  , k  );
87    const double uz_zp = UZ(i  , j  , k+1);
88#endif
89    const double div = 1. / jd_x0 * (
90        - jd_xm / hx_xm * ux_xm + jd_xp / hx_xp * ux_xp
91        - jd_x0 / hy    * uy_ym + jd_x0 / hy    * uy_yp
92#if NDIMS == 3
93        - jd_x0 / hz    * uz_zm + jd_x0 / hz    * uz_zp
94#endif
95    );
96    // check maximum
97    divmax = fmax(divmax, fabs(div));
98  END

Momentum balance

The following relations are defined at the corresponding cell faces.

\[ \begin{align}\begin{aligned}\pder{\vel{1}}{t} = & - \dmomadv{1}{1} - \dmomadv{2}{1} - \dmomadv{3}{1}\\& - \dmompre{1}\\& + \dmomdif{1}{1} + \dmomdif{2}{1} + \dmomdif{3}{1}\\& + \dmombuo.\end{aligned}\end{align} \]
\[ \begin{align}\begin{aligned}\pder{\vel{2}}{t} = & - \dmomadv{1}{2} - \dmomadv{2}{2} - \dmomadv{3}{2}\\& - \dmompre{2}\\& + \dmomdif{1}{2} + \dmomdif{2}{2} + \dmomdif{3}{2}.\end{aligned}\end{align} \]
\[ \begin{align}\begin{aligned}\pder{\vel{3}}{t} = & - \dmomadv{1}{3} - \dmomadv{2}{3} - \dmomadv{3}{3}\\& - \dmompre{3}\\& + \dmomdif{1}{3} + \dmomdif{2}{3} + \dmomdif{3}{3}.\end{aligned}\end{align} \]

The implementations are elaborated below.

Internal energy balance

We define this relation for each cell center \(\left( \ccidx{i}, \ccidx{j}, \ccidx{k} \right)\):

\[ \begin{align}\begin{aligned}\pder{T}{t} = & - \dtempadv{1} - \dtempadv{2} - \dtempadv{3}\\& + \dtempdif{1} + \dtempdif{2} + \dtempdif{3}.\end{aligned}\end{align} \]

The implementations are elaborated below.