Prediction of the moments in advection-diffusion lattice Boltzmann method. I. Truncation dispersion, skewness, and kurtosis

The effect of the heterogeneity in the soil structure or the nonuniformity of the velocity field on the modeled resident time distribution (RTD) and breakthrough curves is quantified by their moments. While the first moment provides the effective velocity, the second moment is related to the longitudinal dispersion coefficient ($k_T$) in the developed Taylor regime; the third and fourth moments are characterized by their normalized values skewness ($\text{Sk}$) and kurtosis ($\text{Ku}$), respectively. The purpose of this investigation is to examine the role of the truncation corrections of the numerical scheme in $k_T$, $\text{Sk}$, and $\text{Ku}$ because of their interference with the second moment, in the form of the numerical dispersion, and in the higher-order moments, by their definition. Our symbolic procedure is based on the recently proposed extended method of moments (EMM). Originally, the EMM restores any-order physical moments of the RTD or averaged distributions assuming that the solute concentration obeys the advection-diffusion equation in multidimensional steady-state velocity field, in streamwise-periodic heterogeneous structure. In our work, the EMM is generalized to the fourth-order-accurate apparent mass-conservation equation in two- and three-dimensional duct flows. The method looks for the solution of the transport equation as the product of a long harmonic wave and a spatially periodic oscillating component; the moments of the given numerical scheme are derived from a chain of the steady-state fourth-order equations at a single cell. This mathematical technique is exemplified for the truncation terms of the two-relaxation-time lattice Boltzmann scheme, using plug and parabolic flow in straight channel and cylindrical capillary with the d2Q9 and d3Q15 discrete velocity sets as simple but illustrative examples. The derived symbolic dependencies can be readily extended for advection by another, Newtonian or non-Newtonian, flow profile in any-shape open-tabular conduits. It is established that the truncation errors in the three transport coefficients $k_T$, $\text{Sk}$, and $\text{Ku}$ decay with the second-order accuracy. While the physical values of the three transport coefficients are set by Péclet number, their truncation corrections additionally depend on the two adjustable relaxation rates and the two adjustable equilibrium weight families which independently determine the convective and diffusion discretization stencils. We identify flow- and dimension-independent optimal strategies for adjustable parameters and confront them to stability requirements. Through specific choices of two relaxation rates and weights, we expect our results be directly applicable to forward-time central differences and leap-frog central-convective Du Fort–Frankel–diffusion schemes. In straight channel, a quasi-exact validation of the truncation predictions through the numerical moments becomes possible thanks to the specular-forward no-flux boundary rule. In the staircase description of a cylindrical capillary, we account for the spurious boundary-layer diffusion and dispersion because of the tangential constraint of the bounce-back no-flux boundary rule.

Figure 1

The first diagram illustrates the three centered average concentration profiles in bimodal channel (composed of open layer and porous layer of porosity $\varphi ={10}^{-1}$) produced by (i) the averaged $\text{ADE}$ (red, dotted line), (ii) the $\text{d2Q5}$ $\text{TRT}$ scheme (blue, solid line), and (iii) the reconstruction from the predicted by the $\text{EMM}$ values (black, dashed line): $k_T=0.94$, ${\text{Sk}}_{☆}=-0.84\sqrt{\text{Pe}}$, and ${\text{Ku}}_{☆}=1.49\text{Pe}$ for $\text{Pe}=13.6$. The second diagram displays similar comparison for the Darcy flow in stratified periodic soil of porosity contrast 50 where $k_T=19.57$, ${\text{Sk}}_{☆}=2.37\sqrt{\text{Pe}}$, and ${\text{Ku}}_{☆}=7.82\text{Pe}$ for $\text{Pe}=9.07$. In the last diagram, the reconstruction procedure involves the four next moments and achieves the entire agreement between the $\text{TRT}$ and $\text{EMM}$. These figures are borrowed [15]: Figs. 24(b), 14(a), and 15(a) there, respectively.

Figure 2

This figure compares the averaged concentration profiles of solute advected by the parabolic velocity field in cylindrical capillary of $R=10$ (l.u.) when $\text{Pe}\approx 9.5$ (left diagram) and $\text{Pe}\approx 95$ (right diagram). The profiles are restored (dashed-dotted line, black) from the predicted values for $\mathcal{U}$, $k_T^{\left(c\right)}$, ${\text{Sk}}_{☆}$, and ${\text{Ku}}_{☆}$ in Table 2. They are compared to Gaussian distribution (solid line, magenta) where ${\text{Sk}}_{☆}={\text{Ku}}_{☆}=0$. The profiles are monitored at $t^{\prime }=D_0/R^{2}t\in [\frac{2\sqrt{3}}{9},\frac{\sqrt{3}}{3},\frac{\sqrt{3}}{2},\frac{2\sqrt{3}}{3}]$.

Figure 3

The first diagram displays together two distribution profiles computed by the $\text{d3Q15}\text{BB}$ with $\mathrm{\Lambda }=\frac{1}{4}$ (black, dotted-dashed line) and $\mathrm{\Lambda }=\frac{1}{12}$ (blue, dotted line) in the same simulations as in Fig. 2 but at earlier time $t^{\prime }=\frac{\sqrt{3}}{6}$. The $\text{d3Q15}\text{BB}$ is applied with $t_c^{\left(a\right)}=t_c^{\left(u\right)}=\frac{1}{2}$, $t_c^{\left(m\right)}=\frac{1}{3}$, $\mathrm{\Lambda }=\frac{1}{4}$, $\text{Pe}\approx 95$. The Gaussian distribution ($k_T=k_T^{\left(c\right)}$, ${\text{Sk}}_{☆}={\text{Ku}}_{☆}=0$) is plotted by solid line (magenta). Three next diagrams plot the relative differences of the numerical dispersion, skewness, and kurtosis with respect to their reference values from Table 1, with $\mathrm{\Lambda }=\frac{1}{4}$ (black, dotted-dashed line) and $\mathrm{\Lambda }=\frac{1}{12}$ (blue, dotted line).

Figure 4

Pure diffusion in straight channel. First diagram plots predicted solution (33) for ${\text{Ku}}_{☆}=\text{Ku}\left(t\right)t=\frac{6A_3}{D_0^2}$ versus $\mathrm{\Lambda }$ with four parameter sets $\{c_e,{\left({\mathrm{\Lambda }}^{-}\right)}^2\}$: (i) $\{\frac{1}{3},\frac{1}{12}\}$ (solid line, red), (ii) $\{\frac{1}{3},\frac{1}{3}\}$ (dotted-dashed line, blue), (iii) $\{\frac{1}{6},\frac{1}{12}\}$ (dotted line, black), and (iv) $\{\frac{1}{6},\frac{1}{3}\}$ (dashed line, magenta). This prediction is $\left\{t_q^{\left(m\right)}\right\}$ and $H$ independent. In two next diagrams, the horizontal lines plot this prediction in cases (a) and (d) where $D_0=\frac{1}{3}\sqrt{\frac{1}{12}}$, when $\mathrm{\Lambda }=\{\frac{1}{12},\frac{1}{6},\frac{1}{4}\}$ [dotted line (black), dotted-dashed line (blue), solid line (red)], and compare it with the numerical results (“filled” symbols). The predicted solution is quasi-exact.

Figure 5

Poiseuille flow in straight channel by $\text{d2Q9}$ scheme. This figure displays the predicted differences in Eq. (23) for truncation result ${\text{Sk}}_{☆}^{\left(\mathrm{tr}\right)}$ and ${\text{Ku}}_{☆}^{\left(\mathrm{tr}\right)}$ relative to the physical values ${\text{Sk}}_{☆}$ and ${\text{Ku}}_{☆}$ from Table 1. This prediction is derived with Eqs. (48) and it is plotted in the limit $c_e\to 0$ for $\mathrm{\Lambda }=\{\frac{1}{12},\frac{1}{6},\frac{1}{4}\}$ [dotted line (black), dotted-dashed line (blue), solid line (red)]. Top row: $t_c^{\left(a\right)}=t_c^{\left(m\right)}=\frac{1}{2}$. Bottom row: $t_c^{\left(a\right)}=t_c^{\left(m\right)}=\frac{1}{3}$. The two first diagrams are plotted for $H=20$, when $\text{Pe}\in [20,3\times {10}^2]$. The two last diagrams are displayed for $\text{Pe}={10}^2$, when $H\in [10,40]$.

Figure 6

Numerical results (symbols) for Poiseuille profile in channel of width $H=20$ are compared to truncation predictions (lines) when velocity weight $t_c^{\left(a\right)}=\frac{1}{2}$ [dashed lines (blue), “circles”] and $t_c^{\left(a\right)}=0$ [dotted line (magenta), “triangles”]. The ${\text{err}}_D$ and $\text{err}\left(\text{Sk}\right)$ in two first columns are independent of the mass weight $\left\{t_q^{\left(m\right)}\right\}$. The $\text{err}\left(\text{Ku}\right)$ is plotted for $t_c^{\left(m\right)}=\frac{1}{2}$ (in the third column) and $t_c^{\left(m\right)}=0$ (in the fourth column). The $\text{d2Q9}\text{SNL}$ is applied for $\text{Pe}\approx 6.3$ with the three parameter sets: top: ${\left({\mathrm{\Lambda }}^{-}\right)}^2=\frac{1}{3}$, $c_e=\frac{1}{4}$; middle: ${\left({\mathrm{\Lambda }}^{-}\right)}^2=\frac{1}{3}$, $c_e=\frac{1}{40}$; bottom: ${\left({\mathrm{\Lambda }}^{-}\right)}^2=\frac{1}{300}$, $c_e=\frac{1}{4}$. The velocity is reduced by a factor of 10 from top diagram to middle and bottom diagrams.

Figure 7

Expt II. Similarly as in Fig. 6 but for $\text{Pe}\approx 9.5$ modeled with ${\mathrm{\Lambda }}^{-}=\sqrt{\frac{1}{12}}$, $c_e=\frac{1}{3}$.

Figure 8

Expt. IV. Similarly as in Figs. 6 and 7 but for $\text{Pe}\approx 95$ modeled with ${\mathrm{\Lambda }}^{-}=\sqrt{\frac{1}{12}}$, $c_e=\frac{1}{30}$.

Figure 9

Expt. IV. Similarly as the second and third diagrams in Fig. 8 but with ${\mathrm{\Lambda }}^{-}=\sqrt{\frac{1}{3}}$, $c_e=\frac{1}{60}$ in the two first diagrams, and ${\mathrm{\Lambda }}^{-}=\frac{1}{8}\times \sqrt{\frac{1}{3}}$, $c_e=\frac{2}{15}$ in the two last diagrams.

Figure 10

Poiseuille flow in cylindrical pipe. This figure displays predicted solution (73) for truncation corrections $\text{err}\left(\text{Sk}\right)$ and $\text{err}\left(\text{Ku}\right)$ in the limit $c_e\to 0$ for three $\mathrm{\Lambda }$ values: $\mathrm{\Lambda }=\frac{1}{12}$ (black, dotted line), $\mathrm{\Lambda }=\frac{1}{6}$ (blue, dotted-dashed line), $\mathrm{\Lambda }=\frac{1}{4}$ (red, solid line). Top row: $t_c^{\left(a\right)}=\frac{1}{2}$, $t_c^{\left(m\right)}=\frac{1}{3}$. Bottom row: $t_c^{\left(a\right)}=t_c^{\left(m\right)}=\frac{1}{3}$. The two first diagrams are plotted for $R=20$, when $\text{Pe}\in [10,2\times {10}^2]$. The two last diagrams are plotted for $\text{Pe}={10}^2$, when $R\in [5,40]$.

Figure 11

This figure presents numerical results (lines with symbols) for the relative velocity error ${\text{err}}_U(t_c^{\left(a\right)},\mathrm{\Lambda })-{\text{err}}_U^{\left(\mathrm{sum}\right)}$ versus $t_c^{\left(a\right)}$, ${\text{err}}_U={\mathcal{U}}^{\left(\mathrm{num}\right)}/\mathcal{U}-1$ due to the bounce-back retardation of the Poiseuille profile in cylindrical pipe on the diagonal velocity stencil $\left\{t_q^{\left(a\right)}\right\}$ for $R=\{5,10,20\}$ from the left to the right. The results are plotted for uniform distribution $\mathrm{\Lambda }=\{\frac{1}{12},\frac{1}{6},\frac{1}{4}\}$ [dotted line (black), dotted-dashed line (blue), solid line (red)]. The theoretical prediction for ${\text{err}}_U^{\left(bb\right)}$ is constructed in work [51].

Figure 12

This figure displays four different components (in row) of the dispersion correction in parabolic flow prescribed in cylindrical capillary of $R=20$, in four numerical experiments Expts. I–IV (from the top to the bottom), when $\mathrm{\Lambda }=\{\frac{1}{12},\frac{1}{6},\frac{1}{4}\}$ [dotted line (black), dotted-dashed line (blue), solid line (red)]. First diagram: numerical result ${\text{err}}_D^{\left(\mathrm{num}\right)}=D^{\left(\mathrm{num}\right)}/\left[D_0(1+k_T^{\left(c\right)})\right]-1$ versus velocity weight $t_c^{\left(a\right)}$. Second diagram: numerical value of the boundary-layer diffusion correction ${\text{err}}_D^{\left(bb\right)}$ given by Eq. (25). Third diagram: predicted [51] boundary-layer dispersion ${\text{err}}_D^{(bb,U)}$. Fourth diagram: numerical result for ${\text{err}}_D^{\left(\mathrm{bulk}\right)}={\text{err}}_D^{\left(\mathrm{num}\right)}-({\text{err}}_D^{\left(bb\right)}+{\text{err}}_D^{(bb,U)})$ (solid lines) and prediction ${\text{err}}_D^{\left(\mathrm{bulk}\right)}={\text{err}}_D^{\left(\mathrm{tr}\right)}+{\text{err}}_D^{\left(\mathrm{sum}\right)}$ (dashed lines) are plotted together. Further details are given in Table 7; the mass weight $t_c^{\left(m\right)}=\frac{1}{3}$ in all simulations.

Figure 13

This figure demonstrates the effect of the anisotropic (except for $\mathrm{\Lambda }=\frac{1}{6}$) truncation dispersion component in numerical results, when $\mathrm{\Lambda }=\{\frac{1}{12},\frac{1}{6},\frac{1}{4}\}$ [dotted line (black), dotted-dashed line (blue), solid line (red)]. The four diagrams are similar to the last column in Fig. 12 but the (anisotropic) coordinate mass weight $t_c^{\left(m\right)}=\frac{1}{2}$ is applied in computations.

Figure 14

Cylindrical pipe of $R=20$, Expt. IV, $t_c^{\left(m\right)}=\frac{1}{3}$, $t_c^{\left(a\right)}=\frac{1}{2}$, uniform-$\mathrm{\Lambda }$ distribution. The first diagram compares numerical results in the form of Eq. (24) (symbols) to their truncation predictions for ${\omega }_2$ [solid line (red), “lozenges”], ${\omega }_3$ [dashed line (blue), “triangles”], and ${\omega }_4$ [dotted line (black), “square”]. The predictions are obtained with the $\text{EMM}$ using either summation in given geometry with $\mathcal{U}\to {\mathcal{U}}^{\left(\mathrm{sum}\right)}$ (thick lines) or integration over the exact circular shape (the thin lines bisecting at $\mathrm{\Lambda }=\frac{1}{6}$). The second diagram displays numerical and $\text{EMM}$ summation results for ${\text{Sk}}_{☆}$ [dashed line (blue), “triangles”] and ${\text{Ku}}_{☆}$ [dotted line (black), “squares”].

Figure 15

Cylindrical pipe of $R=20$, Expt. IV, $t_c^{\left(m\right)}=\frac{1}{3}$, $t_c^{\left(a\right)}=\frac{1}{2}$. Similarly as in Fig. 15 but increasing $\mathcal{U}$ by factor of 4 at fixed $\text{Pe}$.