Simplified method for simulation of incompressible viscous flows inspired by the lattice Boltzmann method

The lattice Boltzmann method (LBM) has gained increasing popularity in incompressible viscous flow simulations, but it uses many distribution functions (far more than the flow variables) and is often memory demanding. This disadvantage was overcome by a recent approach that solves the more actual macroscopic equations obtained through Taylor series expansion analysis of the lattice Boltzmann equations [Lu et al., J. Comput. Phys. 415, 109546 (2020)]. The key is to keep some small additional terms (SATs) to stabilize the numerical solution of the weakly compressible Navier-Stokes equations. However, there are many SATs that complicate the implementation of their method. Based on some analyses and numerous tests, we ultimately pinpoint two essential ingredients for stable simulations: (1) suitable density (pressure) diffusion added to the continuity equation and (2) proper numerical dissipation related to the velocity divergence added to the momentum equations. Then we propose a simplified method that is not only easier to implement but noticeably faster than the original method and the LBM. It contains much simpler SATs that only involve the density (pressure) derivatives, and it requires no intermediate steps or variables. As well, it is extended for thermal flows with small temperature variations and for two-phase flows with uniform density and viscosity. Several test cases, including some two-phase problems under two-dimensional, axisymmetric, and three-dimensional geometries, are presented to demonstrate its capability. This work may help pave the way for the simplest simulation of incompressible viscous flows on collocated grids based on the artificial compressibility methodology.

Figure 1

Evolutions of (a) the error in the horizontal velocity component ${\text{Er}}_u$, (b) the error in the pressure ${\text{Er}}_p$, (c) the total kinetic energy $E_k$, and (d) the total enstrophy $\mathrm{\Omega }$ in the simulation of the 2D Taylor-Green vortex. The Reynolds number is $\text{Re}=20$ and the shared simulation parameters are $N_L=50$ and $N_t=1000$ ($c=20$). The solid line is by the present method, the dashed line is by the original MAMEs and the dash-dot-dot line is by the LBM using MRT. In (c), (d) the insets give the plots of the same data with the axis for $E_k$ and $\mathrm{\Omega }$ in log scale, and the analytical predictions are given by the line with filled squares.

Figure 2

Stable and unstable regions in the $N_L\text{−}c$ plane for different methods to simulate the Taylor-Green vortex problem at $\text{Re}=20$. Note the grid size ${\delta }_x=1/N_L$ and the time step ${\delta }_t=1/N_t=1/\left(c{N}_L\right)$. The three lines approximately represent the critical conditions for the present method (solid), the original MAMEs (dashed) and the LBM using MRT (dash-dot-dot). The stable regions are above the respective lines, and the unstable regions are below them. The filled circles represent the stable cases, and the empty circles represent the unstable cases for the present method (for clarity, the specific data for the other two methods are not plotted). The tests were performed for $0\le t\le 4$. In the stable runs the errors do not increase abruptly. For the LBM runs, the MRT parameters follow those in Fig. 1 of [31], and the errors are quite large when $c$ is too small (even though the runs do not blow up).

Figure 3

Evolutions of (a) the total kinetic energy $E_k$ (b) the error in the horizontal velocity component ${\text{Er}}_u$. The solid line was obtained with both the SATs in Eqs. (10) and (11), the dashed line was obtained with the SATs in the momentum equations [Eq. (11)] only, the dash-dot-dotted line was obtained with the SAT in the continuity equation [Eq. (11)] only, and the long-dashed line was obtained without any SAT. The shared simulation parameters are $N_L=50$ and $N_t=1000$ ($c=20$).

Figure 4

Contour plots of (a), (c) the pressure field and (b), (d) the vorticity field at $t=2.0$ in the simulation of the 2D Taylor-Green vortex at $\text{Re}=20$. The simulation parameters are $N_L=50$ and $N_t=1000$ ($c=20$). The black solid line is by the present method, and the red dashed line is the analytical solution. The results in (a), (b) were obtained with the SATs added in both the continuity and momentum equations, whereas those in (c), (d) were obtained with the SATs in the momentum equations only.

Figure 5

Evolutions of the rate of energy dissipation $\varepsilon \left(t\right)$ at (a) $\text{Re}=100$ and (b) $\text{Re}=200$ for the 3D Taylor-Green vortex. The solid line is by the present SMAMEs simulation, the dashed line is by the original MAMEs, and the dash-dot-dot line is by the 3D LBM (D3Q19) using the weighted MRT collision model [32]. The filled circles are by the spectral methods from [36]. The shared simulation parameters for the SMAMEs, MAMEs, and LBM are $N_L\approx 20.372$ (${128}^3$ grid points) and $N_t\approx 407.437$ ($c=20$).

Figure 6

Contour plots of the velocity components (a), (c) $u$ and (b), (d) $w$ at $t=3.5$ and $\text{Re}=100$ in the plane $z=\pi /4$ for the 3D Taylor-Green vortex. The results by the SMAMEs are in the upper row whereas those by the original MAMEs are in the lower row. The shared simulation parameters for the SMAMEs and MAMEs are $N_L\approx 20.372$ (${128}^3$ grid points) and $N_t\approx 407.437$ ($c=20$).

Figure 7

Velocity profiles along the centerlines (a) $u\left(y\right)$ at $x=0.5$ and (b) $v\left(x\right)$ at $y=0.5$ for the 2D cases; (c) $u\left(z\right)$ at $x=0.5$ and (d) $w\left(x\right)$ at $z=0.5$ in the plane $y=0.5$ for the 3D cases. The reference data for the 2D cases are from [38], and those for the 3D cases are from [39].

Figure 8

Contour plots of the stream function $\mathrm{\Psi }$ (a), (d), the vorticity $-\omega$ (b), (e) and the pressure deviation $\mathrm{\Delta }p$ (c), (f) at steady state for the 2D driven cavity at $\text{Re}=1000$ (upper row) and 5000 (lower row) obtained by the present SMAMEs simulations with $N_L=256$ and $c=50$. The blue lines represent the wall boundaries.

Figure 9

Evolutions of (a) the total enstrophy and (b) the total kinetic energy for the doubly periodic shear layer at $\text{Re}=10000$. The present simulation parameters are $N_L=320$, $N_t=16000$ ($c=50$). For the EDAC results from [14], $N_L=512$ and $N_t=5120$ ($c=10$). Note that ”Pspect” denotes the results obtained by a pseudospectral solve on a $768\times 768$ grid [14, 41].

Figure 10

The vorticity at $t=1.0$ for the doubly periodic shear layer at $\text{Re}=10000$ (a) by the present simulation using SMAMEs and (b) by using the original MAMEs. The simulation parameters are $N_L=320$, $N_t=16000$ ($c=50$).

Figure 11

The pressure at $t=1.0$ for the doubly periodic shear layer at $\text{Re}=10000$ (a) by the present simulation using SMAMEs and (b) by using the original MAMEs. The simulation parameters are the same as in Fig. 10.

Figure 12

Contour plots of the stream function for the 2D natural convection in a square cavity in steady state at different Rayleigh numbers: (a) $\text{Ra}={10}^3$, (b) $\text{Ra}={10}^4$, (c) $\text{Ra}={10}^5$, (d) $\text{Ra}={10}^6$, and $\text{Pr}=0.71$. The shared parameter is $c=20$ and other simulation parameters are $N_L=100$ for $\text{Ra}={10}^3$, $N_L=150$ for $\text{Ra}={10}^4$ and ${10}^5$, and $N_L=200$ for $\text{Ra}={10}^6$. The blue lines represent the wall boundaries.

Figure 13

Contour plots of the temperature for the 2D natural convection in a square cavity in steady state at different Rayleigh numbers: (a) $\text{Ra}={10}^3$, (b) $\text{Ra}={10}^4$, (c) $\text{Ra}={10}^5$, (d) $\text{Ra}={10}^6$, and $\text{Pr}=0.71$. The simulation parameters are the same as in Fig. 12.

Figure 14

Evolutions of the (scaled) interface displacement at the left boundary for the capillary wave at $\text{Re}=1000$. In both (a) and (b) the dashed line is the analytical prediction by Eq. (39). In (a) the solid line is by the present method, and the dash-dotted line is by the MRT-LBM. In (b) the solid line is obtained at $N_L=64$, and the dash-dotted line, $N_L=128$. The shared simulation parameters in (a) are $N_L=64$, $N_t=384$ ($c=6$), $\text{Cn}=0.0625$, $\text{Pe}=2\times {10}^4$. For the fine mesh solution in (b), $N_t=3072$ ($c=24$), $\text{Cn}=0.03125$, $\text{Pe}=2\times {10}^4$.

Figure 15

Evolutions of (a) the centroid velocity $U_{\mathrm{drop}}$ in the axial direction and (b) the aspect ratio ${\alpha }_{\mathrm{drop}}$ of the drop at $\text{Eo}=144$ and $\text{Oh}=0.0466$ for the falling drop problem by the present method, by the axisymmetric LBM in Ref. [28] and by the finite difference solution of the NSEs and the front-tracking method [46]. The present simulation parameters are $N_L=50$, $N_t=4000$ ($c=80$), $\text{Cn}=0.06$, and $\text{Pe}=1000$.

Figure 16

Snapshots of the drop interfaces (red lines) and the stream traces around the drop (in the frame moving with the drop) at three selected times with an interval $\mathrm{\Delta }t/T_c^{\prime }=3.873$ (a) $t/T_c^{\prime }=3.873$, (b) $t/T_c^{\prime }=7.746$, (c) $t/T_c^{\prime }=11.619$ at $\text{Eo}=144$, $\text{Oh}=0.0466$, $r_{\rho }=1.15$, and $r_{\eta }=1$. The dash-dotted line is the $z$ axis. The simulation parameters are as in Fig. 15.

Figure 17

Evolutions of (a) the droplet height $H_d$ (b) the maximum velocity magnitude ${\left|\mathit{u}\right|}_{\max }$ during the drop dewetting/spreading computed by three different methods. The solid lines are by the present 3D simulations, the filled circles are by the present axisymmetric simulations, the dash-dot-dot lines are by the 3D LBM using the weighted MRT collision model [32] and the dashed lines are by the axisymmetric simulation using the VS formulation [47]. The shared simulation parameters are $N_L=32$, $N_t=320$ ($c=10$), $\text{Cn}=0.1$, and $\text{Pe}=8000$.

Figure 18

Contour plots of the velocity components along the $z$ axis and along the radial direction, $u_z$ and $u_r$, at $t=10$ for axisymmetric drop dewetting on a hydrophobic wall. The physical parameters are $\text{Re}=100$ ($\text{Oh}=0.1$) and $\theta ={135}^{\circ }$. The red lines represent the drop interfaces. The shared simulation parameters are $N_L=32$, $N_t=320$ ($c=10$), $\text{Cn}=0.1$, and $\text{Pe}=8000$. The upper row (a), (b) are the present results using the SMAMEs. The lower row (c), (d) are by the axisymmetric VS simulation in Ref. [47].

Figure 19

Contour plots of the velocity components $u_z$ and $u_r$ at $t=10$ for axisymmetric drop spreading on a hydrophilic wall. The physical parameters are $\text{Re}=100$ ($\text{Oh}=0.1$) and $\theta ={60}^{\circ }$. The red lines represent the drop interfaces. The simulation parameters are the same as those in Fig. 18. The upper row (a), (b) are the present results using the SMAMEs. The lower row (c), (d) are by the axisymmetric VS simulation in Ref. [47].

Figure 20

Evolutions of the $z$ component of the centroid velocity of the drop at $\text{Oh}=0.037$ and 0.119. The shared simulation parameters are $N_L=40$, $N_t=400$ ($c=10$), $\text{Cn}=0.1$, $\text{Pe}=8\times {10}^3$.

Figure 21

Snapshots of the drop shapes (in red solid lines) and the contour plots of the velocity component $u$ (in the middle $x\text{−}z$ plane) and $v$ (in the middle $y\text{−}z$ plane) at $t=1.0$, 2.0, 3.0, 4.0, 5.0, and 6.0 [(a)–(f) by the present SMAMEs simulation and (g)–(l) by the 3D MRT-LBM simulation] at $\text{Oh}=0.037$ for the coalescence induced drop jumping on a nonwetting wall. In each panel, the left part shows the middle $y\text{−}z$ plane at $x=0$ and the right shows the middle $x\text{−}z$ plane at $y=0$. The shared simulation parameters are as in Fig. 20.

Figure 22

Snapshots of the drop shapes (in red solid lines) and the contour plots of the velocity component $w$ in the middle $y\text{−}z$ and $x\text{−}z$ planes at $t=1.0$, 2.0, 3.0, 4.0, 5.0, and 6.0 [(a)–(f) by the present SMAMEs simulation and (g)–(l) by the 3D MRT-LBM simulation) at $\text{Oh}=0.037$ for the coalescence induced drop jumping on a nonwetting wall. In each panel, the left part shows the middle $y\text{−}z$ plane at $x=0$ and the right shows the middle $x\text{−}z$ plane at $y=0$. The shared simulation parameters are as in Fig. 20.