2.1. Time integration

The Caputo fractional time derivative in (2.1)–(2.2) is discretized using the difference approximation and studied by Podlubny [Reference Podlubny35]. Suppose the time interval [0, T] is discretized uniformly into n subintervals; we define ${t_k = k \Delta t, \,\, k = 0, 1,\ldots ,n}$ , where $\Delta t = T/n$ is the time step. Let $u(\mathbf {x},\,t_k)$ be the exact value of a function $u(\mathbf {x}, t)$ at time step $t_k$ . Then, the fractional derivative can be approximated as

$$ \begin{align*} \frac{\partial^\alpha u(\mathbf{x}, t)}{\partial t^\alpha}\bigg\vert_{t = t_{k+1}} &\approx \frac{1}{\Gamma(1-\alpha)} \sum^k_{j=0} \frac{u(\mathbf{x},\,t_{j+1})-u(\mathbf{x},\,t_j)}{\Delta t} \int^{(j+1)\Delta t}_{j\Delta t} \frac{d t'}{(t_{k+1}-t')^\alpha} \nonumber \\ &=\frac{1}{\Gamma(1-\alpha)} \sum^k_{j=0} \frac{u(\mathbf{x},\,t_{j+1})-u(\mathbf{x},\,t_j)}{\Delta t} \int^{(k-j+1)\Delta t}_{(k-j)\Delta t} t^{\prime -\alpha} \,d t' \nonumber \\ &=\frac{1}{\Gamma(1-\alpha)} \sum^k_{j=0} \frac{u(\mathbf{x},\,t_{k+1-j})-u(\mathbf{x},\,t_{k-j})}{\Delta t} \int^{(j+1)\Delta t}_{j\Delta t} t^{\prime -\alpha} \,d t' \nonumber \\ &=\frac{(\Delta t)^{1-\alpha}}{\Gamma(2-\alpha)} \sum^k_{j=0} \frac{u(\mathbf{x},\,t_{k+1-j})-u(\mathbf{x},\,t_{k-j})}{\Delta t} r^\alpha_j, \end{align*} $$

where the weight, $r^\alpha _j=[ (j+1)^{1-\alpha } - j^{1-\alpha } ]$ . Following an earlier work by the authors that used standard integer-order advection–diffusion-reaction (ADR) equations [Reference Singh, Bansal, Kaur and Sircar41], we retain the “IMEX method philosophy” by explicitly discretizing the advection and the reaction terms (that is, the “ $K_1\nabla \cdot u(\mathbf {x}, t)$ ” term and “ $f(\mathbf { x}, t)$ ” term, respectively, in (2.1)), while the diffusive term (that is, the “ $K_2(\mathbf {x}, t) \nabla ^2 u(\mathbf {x}, t)$ ” term in (2.1)) is treated semi-implicitly as follows:

$$ \begin{align*} K_2(\mathbf{x}, t) \nabla^2 u(\mathbf{x}, t) \approx \theta K_2(\mathbf{x}, t) \nabla^2 u(\mathbf{x}, t)|_{t_k} + (1-\theta) K_2(\mathbf{x}, t) \nabla^2 u(\mathbf{x}, t)|_{t_{k-1}}. \end{align*} $$

This time integration method (referred to as the $\theta $ -method in the subsequent discussion) generalizes the computationally explicit $L1$ -method by Brunner et al. [Reference Brunner, Ling and Yamamoto10] as well as the fully implicit method recently proposed by Jannelli [Reference Jannelli23].

2.2. Adaptive time stepping

When the simulation time is long, the size of the “memory” in the fractional derivative approximation becomes enormously large. However, according to the “fading memory property” [Reference Diethelm and Freed16], for long times, the solution of the FADR systems changes slower than the standard integer-order ADR processes. Hence, it makes sense to employ a large time step at longer times. Let $\widetilde {u(\cdot ,\, t_k)}$ be the numerical approximation for $u(\cdot ,\, t_k)$ . To detect this change, we define a measure between the numerical solutions of two consecutive time steps by

$$ \begin{align*} \Delta u_{t_k} = \frac{\| \widetilde{u(\cdot,\, t_k)} - \widetilde{u(\cdot,\, t_{k-1})} \|_{l^2}}{|| \widetilde{u(\cdot,\, t_{k-1}) ||_{l^2}}}, \quad k = 1, \ldots, n. \end{align*} $$

For some user-defined relaxation parameter $\delta $ , if $\Delta u_{t_k} < \delta $ , then the time spacing is geometrically increased (that is, $\Delta t \rightarrow 2\Delta t$ ), up to some prefixed value $\Delta t_{\text {max}}$ .

2.3. Spatial approximation

The gradients and Laplacian in (2.1) are spatially approximated using the upwind difference scheme (UDS) [Reference Rai and Moin36] and the second-order central difference scheme (CDS), respectively. Recall that the CDS approximation of the convective terms in (2.1) does not model the flow-physics accurately [Reference Singh, Bansal, Kaur and Sircar41]. Furthermore, a first-order upwind scheme is too diffusive, thereby necessitating the use of higher order upwind schemes. For example,

$$ \begin{align*} K_1(\mathbf{x}, t) \frac{\partial u}{\partial x} \approx {K_1}_{ij} \bigg( \frac{u^k_{i+1,j}-u^k_{i-1,j}}{2\Delta x} \bigg) + q(K^+_1 u^-_x + K^-_1 u^+_x), \end{align*} $$

where

$$ \begin{align*} \begin{array}{lll} &K^-_1 = \text{min}({K_1}_{ij}, 0), &K^+_1 = \text{max}({K_1}_{ij}, 0), \nonumber\\[6pt] &u^-_x = \dfrac{u^k_{i-2,j} - 3u^k_{i-1,j} + 3u^k_{i,j} - u^k_{i+1,j}}{3 \Delta x}, & u^+_x = \dfrac{u^k_{i-1,j} - 3u^k_{i,j} + 3u^k_{i+1,j} - u^k_{i+2,j}}{3 \Delta x}. \end{array} \end{align*} $$

Here the parameter $q=0.5$ represents the third-order accurate upwind formula (UD3) that is used for the interior stencil points. The use of ghost points is avoided by setting $q = 0$ for grid points immediately adjacent to the boundary, leading to a second-order accurate method (UD2) at these points. Since the focus of the present work is on the development of the time integration method, we retain the same spatial approximation outlined above in all the subsequent test cases.

3.1. Time-asymptotic stability analysis

First, we show that the solution of the linear 1D FADR equation, coupled with periodic boundary conditions and discretized using the numerical method outlined in Sections 2.1–2.3, is time-asymptotically stable for a limited range of CFL and Péclet numbers [Reference Sengupta, Sengupta, Puttam and Vajjala39]. The discretization of the linear 1D FADR equation, using the $\theta $ -method, is as follows:

(3.1) $$ \begin{align} &\frac{(\Delta t)^{-\alpha}}{\Gamma{(2-\alpha)}}\bigg\{ \sum_{j=1}^n (j^{(1-\alpha)}-(j-1)^{(1-\alpha)}) (u_i^{n-j+1}-u_i^{n-j})\bigg\} \nonumber\\ &\qquad + c \bigg(\frac{u_{i+1}^{n-1}-u_{i-1}^{n-1}}{2\Delta x}\bigg) + qc \bigg(\frac{u_{i-2}^{n-1}-3u_{i-1}^{n-1}+3u_{i}^{n-1}-u_{i+1}^{n-1}}{3\Delta x}\bigg) \nonumber \\ &\quad = \gamma \theta \frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{(\Delta x)^2} + \gamma(1-\theta)\frac{u_{i+1}^{n-1}-2u_{i}^{n-1} + u_{i-1}^{n-1}}{(\Delta x)^2} + \lambda u_{i}^{n-1}, \end{align} $$

where $n / i$ denotes the temporal/spatial index, and we set the parameter $q=0.5$ . We introduce the following nondimensional parameters: the fractional CFL number $N_c$ , fractional Péclet number $Pe$ and fractional Damköhler number $Da$ , as

$$ \begin{align*} N_c = \frac{c (\Delta t)^\alpha}{\Delta x}, \quad Pe = \frac{\gamma (\Delta t)^\alpha}{(\Delta x)^2}, \quad Da = \frac{\lambda (\Delta x)}{c}. \end{align*} $$

In the subsequent discussion, we drop the prefix “fractional” in these nondimensional parameters. Rearranging the terms in (3.1), we arrive at the following equation:

(3.2) $$ \begin{align} &\{1\!+\! 2 Pe \theta \Gamma{(2\!-\! \alpha)}\} u_i^n \! -\! Pe \theta \Gamma{(2\!-\!\alpha)}u_{i+1}^n \! -\! Pe\theta\Gamma{(2\!-\! \alpha)}u_{i-1}^n \nonumber\\ &\quad= \{1-qN_c\Gamma{(2\!-\! \alpha)} -2 Pe(1\!-\!\theta)\Gamma{(2\!-\!\alpha)} + Da N_c\Gamma{(2-\alpha)}\}u_i^{n-1} \nonumber\\ &\qquad + \bigg\{-N_c\frac{\Gamma{(2-\alpha)}}{2}+q N_c\frac{\Gamma{(2-\alpha)}}{3} + Pe(1\!-\!\theta)\Gamma{(2\!-\!\alpha)}\bigg\}u_{i+1}^{n-1} \nonumber\\ &\qquad + \bigg\{N_c\frac{\Gamma{(2\!-\!\alpha)}}{2} + q N_c\Gamma{(2\!-\!\alpha)} \! +\! Pe (1-\theta)\Gamma{(2-\alpha)} \bigg\}u_{i-1}^{n-1} \nonumber \\ &\qquad - q N_c\frac{\Gamma{(2-\alpha)}}{3}u_{i-2}^{n-1} + \sum_{j=2}^n \{j^{(1-\alpha)}-(j-1)^{(1-\alpha)}\} (u_i^{n-j+1} - u_i^{n-j}). \end{align} $$

If $\widetilde {u_i^n}$ is the (numerically) approximate solution of (3.2), then the round-off error, $\varepsilon _i^n = u_i^n - \widetilde {u_i^n}\quad (i = 0, \ldots , M)$ , identically satisfies the same equation. Due to periodic boundary conditions,

$$ \begin{align*} \varepsilon_0^n=\varepsilon_M^n, \quad n = 1,2,\ldots,T. \end{align*} $$

Assume the round-off error has the form

$$ \begin{align*} \varepsilon_i^n = \xi_n e^{Ik(ih)}, \end{align*} $$

where $\xi _n = |\varepsilon _i^n|,\, I=\sqrt {-1},\, h=\Delta x$ and $k = {2\pi l}/{L}$ (l, index and L, spatial domain length). We substitute the above relation in the equation for the round-off error to arrive at the form

(3.3) $$ \begin{align} \mu_1 \xi_n = \mu_2 \xi_{n-1} +\sum_{j=2}^n r^\alpha_j (\xi_{n-j+1}-\xi_{n-j}), \end{align} $$

where

(3.4) $$ \begin{align} \mu_1 &= 1+2Pe\theta\Gamma(2-\alpha)(1-\cos(kh)), \nonumber \\ \mu_2 &= 1 - \bigg \{ N_cq \bigg(1 - \frac{2}{3}\cos(kh)\bigg) + 2 Pe(1 - \theta)(1 - \cos(kh)) - DaN_c + IN_c\sin(kh) \nonumber \\ & \qquad\quad - \frac{N_c q}{3} (2e^{-Ikh} - e^{-2Ikh}) \bigg \}\Gamma(2 - \alpha), \nonumber \\ r^\alpha_j &= j^{(1-\alpha)}-(j-1)^{(1-\alpha)}. \end{align} $$

Next, we prove the result that the finite difference scheme in (3.3)–(3.4) is time- asymptotically stable in the following theorem.

Proof. It suffices for us to show that the time-amplification factor, $\xi _n$ in (3.3), obeys the inequality,

$$ \begin{align*} \xi_n \le \xi_{n-1} \le \xi_{n-2} \le \cdots \le \xi_1 \le \xi_0. \end{align*} $$

• • For $n=1$ in (3.3)–(3.4), we observe the $\mu _1 \ge 1$ while $\mu _2 \le 1$ for sufficiently small $N_c$ and $Pe$ , thereby leading us towards the conclusion that $\xi _1 \le \xi _0$ .

• • Using the induction hypothesis, we assume that $\xi _{n-1} \le \xi _{n-2} \le \cdots \le \xi _1 \le \xi _0$ . Next we show that $\xi _n \le \xi _{n-1}$ .

• • Since $\mu _1 \ge 1$ and $\mu _2 \le 1$ , (3.3) can be replaced with the following inequality:

  (3.5) $$ \begin{align} \xi_n \le \xi_{n-1} +\sum_{j=2}^n r^\alpha_j (\xi_{n-j+1}-\xi_{n-j}) \le \xi_{n-1}. \end{align} $$
  The second inequality in (3.5) occurs because:

  1. (i) $r^\alpha _j$ is positive;

  2. (ii) $r^\alpha _j> r^\alpha _{j+1}$ ;

  3. (iii) $\xi _{n-j+1}-\xi _{n-j} \le 0$ (via the induction hypothesis).

  This completes the proof.

We emphasize that while the $\theta $ -method is not unconditionally stable, it is time-asymptotically stable for a restricted range of $N_c$ and $Pe$ .

3.2. Spectral analysis

Although the $\theta $ -method is time-asymptotically stable, the presence of dispersion errors (through negative group velocity and large phase speed errors) would invalidate the long time integration results [Reference Sengupta, Sircar and Dipankar40]. Hence, we couple the result in Section 3.1 along with the toolset of spectral analysis to deduce the relevant range of the parameters, $N_c, Pe$ and $Da$ , for an accurate representation of the numerical solution of the 1D FADR equation.

Using the spectral (Fourier) representation of the approximate solution of (3.2), we have $\widetilde {u_i^n} = \xi ^{\prime }_n e^{I(ikh)}\,\,\,(I=\sqrt {-1},\, k={2\pi l}{L})$ . We define the numerical time-amplification factor

$$ \begin{align*}G_{\text{num}} = \frac{\widetilde{u_i^n}}{\widetilde{u_i^{n-1}}} = \frac{\xi^{\prime}_n}{\xi^{\prime}_{n-1}}.\end{align*} $$

Dividing (3.2) by $\widetilde {u_i^{n-1}}$ , we arrive at the equation governing $G_{\text {num}}$ ,

(3.6) $$ \begin{align} C_0 G_{\text{num}}^n + C_1 G_{\text{num}}^{n-1}+ \sum_{j=2}^n r_j^\alpha (G_{\text{num}}^{n-j+1}-G_{\text{num}}^{n-j}) = 0, \end{align} $$

where the coefficients,

$$ \begin{align*} C_0 &= 1+2Pe\theta \Gamma{(2-\alpha)}(1-\cos(kh)), \nonumber\\ C_1 &= -1+qN_c\Gamma(2-\alpha)\bigg(1-\frac{4}{3}\cos(kh)+\frac{1}{3}\cos(2kh)\bigg) \nonumber \\ &\quad +2Pe(1-\theta)\Gamma(2-\alpha)(1-\cos(kh)) -Da N_c\Gamma(2\!-\! \alpha) \nonumber \\ &\quad +IN_c\Gamma(2-\alpha)\bigg(\!\sin(kh)+\frac{2}{3}q\sin(kh)-\frac{q}{3}\sin(2kh)\bigg). \end{align*} $$

We remark that the order “n” of the polynomial in (3.6) is fixed at $n=75$ in the subsequent analysis, since our numerical studies have shown that an increase in the polynomial order by one shifts the contours of the spectral variables by less than $0.001\%$ .

Next, transforming the exact solution of the linear 1D FADR equation, using the Fourier–Laplace transform as

$$ \begin{align*}u(x, t) = \iint\hat{U}(k, \omega) e^{\mathit{I}(kx - \omega t)}, \end{align*} $$

we arrive at the exact dispersion relation,

$$ \begin{align*} \omega_{\text{exact}} = \omega = I(\lambda-\gamma k^2-c I k)^{{1}/{\alpha}}. \end{align*} $$

Similarly, we have a corresponding numerical dispersion relation for the approximate equation (3.2),

$$ \begin{align*} \omega_{\text{num}} = I(\lambda^N-\gamma^N k^2-c^N I k)^{{1}/{\alpha}}, \end{align*} $$

where the superscript “N” denotes the corresponding numerical values of the parameters, and the subscripts “exact” and “num” indicate the exact (analytical) and the numerical value of the variable, respectively. Since

$$ \begin{align*} G_{\text{num} }= \frac{u(k, t+\Delta t)}{u(k, t)} = \frac{e^{I(kx-\omega_{\text{num}} (t+\Delta t))}}{e^{I(kx-\omega_{\text{num}} t)}} = e^{-I\omega_{\text{num}}\Delta t} = e^{I \beta}, \end{align*} $$

and the numerical phase speed is given by $c_{\text {num}} = {\omega _{\text {num}}}/{k}$ ,

$$ \begin{align*} c_{\text{num}} = -\frac{\beta}{k \Delta t} = -\frac{1}{k \Delta t}\tan^{-1} \bigg( \frac{(G_{\text{num}})_{Im}}{(G_{\text{num}})_{Re}} \bigg), \end{align*} $$

where the subscript “ $Im$ ” and “ $Re$ ” denote the imaginary and real values, respectively. Using the expression for the exact phase speed, $c_{\text {exact}} = {\omega _{\text {exact}}}/{k}$ and the non-dimensional parameters described in (3.2), we find the ratio of the phase speeds,

$$ \begin{align*} \frac{c_{\text{num}}}{c_{\text{exact}}} = \frac{-I\beta}{(DaN_c-(kh)^2Pe-I(kh)N_c)^{{1}/{\alpha}}}. \end{align*} $$

Finally, the expression for the numerical and the exact group velocities $V_g$ are given by

$$ \begin{align*} \begin{aligned} {V_g}_{\text{num}}& = \bigg[\frac{\partial}{\partial k} (\omega_{\text{num}})\bigg]_{Re} = \bigg[\frac{\partial}{\partial k} \bigg(\frac{\beta}{\Delta t}\bigg)\bigg]_{Re} = \bigg[\frac{h}{\Delta t}\frac{d\beta}{d(kh)}\bigg]_{Re} , \nonumber\\[.2cm] {V_g}_{\text{exact}} &= \bigg[\frac{\partial}{\partial k} (\omega_{\text{exact}})\bigg]_{Re} = \bigg[\frac{\partial}{\partial k} (I(\lambda\!-\!\gamma k^2 \!-\! cI k)^{{1}/{\alpha}})\bigg]_{Re} =\bigg[\frac{1}{\alpha} (c\!-\!2k\gamma I)(-I\omega)^{1\!-\!\alpha}\bigg]_{Re}. \end{aligned} \end{align*} $$

Hence,

$$ \begin{align*} \frac{(V_g)_{\text{num}}}{(V_g)_{\text{exact}}} = \bigg( \frac{\alpha}{N_c(r^{1-\alpha}(\Delta t)^{1-\alpha}\cos{(1-\alpha)\phi}){+2 kh (r^{1-\alpha}(\Delta t)^{1-\alpha}\sin{(1-\alpha)\phi})Pe}} \bigg) \bigg(\frac{d\beta}{d(kh)}\bigg)\bigg|_{Re}, \end{align*} $$

where $\omega = r e^{I (\phi +{\pi }/{2})}$ .

The contour plots for the group velocity ratio, $V_g = [{V_{g, {\text {num}}}}/{V_{g, \text {exact}}}]_{Re}$ , for Péclet numbers, $Pe = 0.001, 0.01$ and $1.0$ , and Damköhler numbers, $Da = -0.01, 0.0$ and $0.01$ , are presented for two values of $\theta $ , namely, $\theta =0.5$ (in Figure 1) and $\theta =1.0$ (in Figure 2) in the $(N_c, kh)$ -plane. The corresponding contour plots for the absolute phase speed error, $\Delta c = |1 - {c_{\text {num}}}/{c_{\text {exact}}}|$ , are shown in Figures 3 and 4 for $\theta =0.5$ and $\theta =1.0$ , respectively. Regions of favourable spectral properties ( $V_g> 0$ and $\Delta c \le 0.1$ ) are highlighted with grey colour.

Figure 1 Group velocity ratio contours in the $N_c$ – $k_h$ plane ( $N_c$ plotted on the x-axis), $V_g$ , at $\theta =0.5$ , $\alpha =0.9$ and (a) $Pe=0.001,~Da=-0.01$ , (b) $Pe=0.001,~Da=0.0$ , (c) $Pe=0.001,~Da=0.01$ , (d) $Pe=0.01,~Da=-0.01$ , (e) $Pe=0.01,~Da=0.0$ , (f) $Pe=0.01,~Da=0.01$ , (g) $Pe=1.0,~Da=-0.01$ , (h) $Pe=1.0,~Da=0.0$ and (i) $Pe=1.0,~Da=0.01$ .

Figure 2 Group velocity ratio contours in the $N_c$ – $k_h$ plane ( $N_c$ plotted on the x-axis), $V_g$ , at $\theta =1.0$ , $\alpha =0.9$ and (a) $Pe=0.001,~Da=-0.01$ , (b) $Pe=0.001,~Da=0.0$ , (c) $Pe=0.001,~Da=0.01$ , (d) $Pe=0.01,~Da=-0.01$ , (e) $Pe=0.01,~Da=0.0$ , (f) $Pe=0.01,~Da=0.01$ , (g) $Pe=1.0,~Da=-0.01$ , (h) $Pe=1.0,~Da=0.0$ and (i) $Pe=1.0,~Da=0.01$ .

Figure 3 Absolute phase speed error contours in the $N_c$ – $k_h$ plane ( $N_c$ plotted on the x-axis), $\Delta c$ at $\theta =0.5$ , $\alpha =0.9$ and (a) $Pe=0.001,~Da=-0.01$ , (b) $Pe=0.001,~Da=0.0$ , (c) $Pe=0.001,~Da=0.01$ , (d) $Pe=0.01,~Da=-0.01$ , (e) $Pe=0.01,~Da=0.0$ , (f) $Pe=0.01,~Da=0.01$ , (g) $Pe=1.0,~Da=-0.01$ , (h) $Pe=1.0,~Da=0.0$ and (i) $Pe=1.0,~Da=0.01$ .

Figure 4 Absolute phase speed error contours in the $N_c$ – $k_h$ plane ( $N_c$ plotted on the x-axis), $\Delta c$ at $\theta =1.0$ , $\alpha =0.9$ and (a) $Pe=0.001,~Da=-0.01$ , (b) $Pe=0.001,~Da=0.0$ , (c) $Pe=0.001,~Da=0.01$ , (d) $Pe=0.01,~Da=-0.01$ , (e) $Pe=0.01,~Da=0.0$ , (f) $Pe=0.01,~Da=0.01$ , (g) $Pe=1.0,~Da=-0.01$ , (h) $Pe=1.0,~Da=0.0$ and (i) $Pe=1.0,~Da=0.01$ .

We outline the spectral properties, namely, the group velocity ratio and the absolute phase speed error, for a specific case of the fractional order, $\alpha = 0.9$ . In general, we find that in the limit of vanishingly small values of $kh$ , the numerical method has favourable spectral properties (that is, $V_g> 0$ and $\Delta c \le 0.1$ in the limit, $kh \rightarrow 0$ ). The region of positive group velocities ( $V_g> 0$ ) indicate a region where the numerical solution travels in the correct direction and the numerical instabilities in the form of q-waves are avoided [Reference Singh, Bansal, Kaur and Sircar41]. Figures 1 and 2 indicate that this favourable region is restricted to smaller values of $kh$ with larger values of $Pe$ as well as with smaller values of $\theta $ . While the former observation can be attributed to the fact that a stiff diffusive term (that is, larger “ $K_2(\mathbf {x}, t) \nabla ^2 u(\mathbf {x}, t)$ ” term in (2.1)) can be correctly approximated with smaller grid-size, h; the latter observation can be explained through a numerical stabilization due to the implicit treatment of fast time scale term (in this case, the diffusive term). Both of these observations are in congruence with the corresponding analysis of the integer order ADR equations [Reference Singh, Bansal, Kaur and Sircar41]. Similarly, the absolute phase error contours highlight a spectrally favourable region ( $\Delta c \le 0.1$ , equivalently the phase errors are restricted to 10% of the exact phase speed, see Figures 3 and 4) at lower values of $Pe$ and at lower values of $\theta $ , indicating that the phase errors in this numerical method can be reduced by introducing finer grids [Reference Jin, Li and Zhou24, Reference Jin and Li25].

We summarize our discussion by indicating two sources of error that are particularly perplexing in the direct numerical simulation (DNS) of FADR equations. The first source of error is the existence of q-waves for those set of numerical parameters for which the spatiotemporal discretization is time-asymptotically stable. In such a situation, the q-waves do not attenuate and have to be eliminated by deploying an explicit filter [Reference Visbal and Gaitonde54]. Another aspect of spectral error is related to Gibbs’ phenomenon that occurs as a consequence of sharp discontinuity in the solution and that causes fictitious oscillations, a problem that can be remedied using high accuracy dispersion relation preserving schemes (which has at least shown promise in integer order partial differential equations) [Reference Sengupta, Sircar and Dipankar40].

5.1. Initial and boundary conditions

Initial conditions: a rectilinear coordinate system is used with $x, y$ denoting the channel flow direction and the transverse direction, respectively. The origin of this coordinate system is chosen at the left end of the lower wall of the channel. The mean flow is assumed to be a plane Poiseuille flow with its variation entirely in the transverse direction, namely,

$$ \begin{align*} \mathbf{U}_0 = (y - y^2) \mathbf{e_x}, \end{align*} $$

where $\mathbf {e_x}$ is the unit vector along the x-direction. The mean flow, $\mathbf {U}_0$ , defines the mean vorticity, $\Omega _0 = 2y-1$ , and the mean streamfunction, $\psi _0 = {y^2}/{2} - {y^3}/{3}$ . The base state polymer conformation tensor, $\mathcal {C}_0 = [{\mathcal {C}_0}_{i j}]$ , is given by

$$ \begin{align*} \begin{aligned} {\mathcal{C}_0}_{11} &= 1, \nonumber \\ {\mathcal{C}_0}_{12} &= {\mathcal{C}_0}_{21} = We (1 - 2y), \nonumber \\ {\mathcal{C}_0}_{22} &= 2 We^2 (1 - 2y)^2 + 1. \end{aligned} \end{align*} $$

The initial condition comprises the mean flow superposed with a perturbed unstable mode, as follows:

(5.2) $$ \begin{align} \Omega|_{t=0} &= \Omega_0 + \epsilon \Omega_1 = \Omega_0 + \epsilon \mathcal{R}\{\widetilde{\Omega}(y)|_{t=0} e^{I k x} \}, \nonumber \\ \psi|_{t=0} &= \psi_0 + \epsilon \psi_1 = \psi_0 + \epsilon \mathcal{R}\{\widetilde{\psi}(y)|_{t=0} e^{I k x} \}, \nonumber\\ \mathcal{C}|_{t=0} &= \mathcal{C}_0 + \epsilon \mathcal{C}_1 = \mathcal{C}_0 + \epsilon \mathcal{R} \{\widetilde{\mathcal{C}_1}(y)|_{t=0} e^{I k x} \}, \end{align} $$

where $(\Omega _1, \psi _1, \mathcal {C}_1)$ are the perturbations that are Fourier transformed in the x-direction. Here, $\mathcal {R}\{ \}$ denotes the real part of the complex valued function. The equations governing the initial conditions in (5.2) are listed in Appendix A., whose solution is found using a standard MATLAB boundary value solver.

Boundary conditions: to imitate an infinitely long channel, periodic boundary conditions are assumed at the flow inlet and outlet. No-slip (that is, $u = v = 0$ ) and zero tangential conditions (that is, ${\partial u}/{\partial x} = {\partial v}/{\partial x} = 0$ ) are imposed on the lower wall ( $y = 0$ ) and the upper wall ( $y = 1.0$ ) of the channel, respectively. Further, the incompressibility constraint provides an additional condition on the walls: ${\partial v}/{\partial y} = 0$ . Since the flow is parallel to the channel walls, the walls may be treated as streamline. Thus, the streamfunction value, $\psi $ , on the wall is set as a constant. That constant (which may be different on the lower and the upper wall) is found from the no-slip condition. The zero tangential condition implies that all tangential derivatives of streamfunction vanish on the wall. Thus, the boundary condition for vorticity is found from the Poisson equation (5.1b),

$$ \begin{align*} \frac{\partial^2 \psi}{\partial y^2}|_{\text{boundary}} = -\Omega_{\text{boundary}}, \end{align*} $$

where the subscript “boundary” denotes the value of the variable at the boundary. Finally, the boundary conditions for the elastic stress tensor are constructed from (5.1c), coupled with the no-slip and zero tangential conditions, as follows:

$$ \begin{align*} \begin{aligned} & \mathcal{C}_{11} = 1, \nonumber \\ & \frac{\mathcal{C}_{12}}{We} + \frac{\partial^\alpha \mathcal{C}_{12}}{\partial t^\alpha} - \frac{\partial u}{\partial y} = 0, \nonumber \\ & \frac{\mathcal{C}_{22}-1}{We} + \frac{\partial^\alpha \mathcal{C}_{22}}{\partial t^\alpha} - 2 \frac{\partial u}{\partial y} \mathcal{C}_{12} = 0. \end{aligned} \end{align*} $$

5.2. Algorithmic details

The size of the domain is chosen to be $(x, y) \in [0,\,5] \times [0,\,1]$ . The domain is discretized using $76 \times 51$ points, such that the discrete points are equally spaced at $\Delta x = {5}/{75}$ and $\Delta y ={1}/{50}$ in the x- and the y-directions, respectively. The variable in the $\theta $ -method is fixed at $\theta =1.0$ . The minimum and the maximum time steps are chosen as $\Delta t_{\text {min}}=10^{-3}$ and $\Delta t_{\text {max}}=1.6 \times 10^{-2}$ , respectively. The parameters in the initial conditions in (5.2) are set at $\epsilon = 0.1, k = 0.1$ . The time step in the simulation is increased adaptively using the technique outlined in Section 2.2. The Poisson equation (5.1b) is iteratively solved using the Gauss–Seidel iteration technique [Reference Atkinson3]. Assuming $(N+1)$ (or $(M+1)$ ) points in the x- (or y-) direction, the size of the coefficient matrix is $N(M-1) \times N(M-1)$ . Hence, due to size constraints, the numerical inversion of this coefficient matrix imposes severe memory restrictions. Instead, we note that the coefficient matrix has the following block tridiagonal structure with $(M-1)$ blocks:

(5.3) $$ \begin{align} \begin{pmatrix} D_1 & D_2 & 0 & \cdots & 0 \\ D_2 & \ddots & \ddots & & \vdots \\ 0 & \ddots & \ddots & \ddots& 0\\ \vdots & & \ddots & \ddots & D_2\\ 0 & \cdots& 0 & D_2 & D_1 \end{pmatrix} \end{align} $$

where $D_1$ ( $D_2$ ) are $N \times N$ tridiagonal (diagonal) matrices, such that $D_1$ = Tridiag( $r_1$ , R, $r_1$ ), $D_2 = \text {Diag}~(r_2)$ , $r_1 = (\Delta x)^{-2}, r_2 = (\Delta y)^{-2}$ and $R = -2(r_1 + r_2)$ . Since the exact location of the nonzero entries is known, the invocation of the entire coefficient matrix in the Guass–Siedel iteration is avoided. The solution is updated by using the nonzero entries of each row on a case-by-case basis. The row diagonal dominance of the matrix in (5.3) ensures that the iteration converges in a finite number of steps. Similarly, the vorticity equation (5.1a) involves another Laplacian term and its discretization involves a coefficient matrix identical to the form in (5.3) with

$$ \begin{align*}r_1 = -\nu \theta (dx)^{-2}, \quad r_2 = -\nu \theta (dy)^{-2}, \quad R = -2(r_1 + r_2) + \frac{Re (\Delta t)^{-\alpha}}{\Gamma(2-\alpha)}.\end{align*} $$

Thus, the Laplacian term in (5.1a) is treated identically as above.

5.3. Numerical evolution of the conformation tensor equations

The advection term in the conformation tensor equation (5.1c) is second-order accurate in space with the exception of a few points where it is only first-order accurate. This change in order is due to the special treatment of the advection term in the conformation tensor that is needed to avoid the Hadamard instabilities resulting in the loss of the symmetric, positive definite (SPD) property of the conformation tensor [Reference Sureshkumar, Beris and Handler52]. The special treatment is an efficient transformation of the slope-limiting approach of Vaithianathan et al. [Reference Vaithianathan, Robert, Brasseur and Collins53], which requires the evaluation of three schemes: forward, backward and centred stencils, to approximate the advective flux at each grid point and at each time step. The resultant scheme that maximizes the eigenvalues of the conformation tensor at the grid points is then chosen. Further, following Dubief et al. [Reference Dubief, Terrapon, White, Shaqfeh, Moin and Lele18], we ensure that the polymer does not exceed its maximum extensibility by using a semi-implicit approach for the relaxation term in (5.1c).

5.4. Quantifying viscoelastic macrostructures via non-Euclidean metric

The statistics of polymer forces and torques have been used in previous methods to visualize the dynamics of polymers in diluted liquids [Reference Dubief, Terrapon, White, Shaqfeh, Moin and Lele18, Reference Kim and Sureshkumar26]. However, the conformation tensor, $\mathcal {C}$ , a second-order positive definite tensor that is generated by averaging the dyad formed by the polymer end-to-end vector over all molecular realizations, is a better indicator of the polymer deformation history [Reference Bird, Armstrong and Hassager9]. The trace of $\mathcal {C}$ (referred to as $\text {tr} \mathcal {C}$ from now on) is frequently used in the literature to analyse $\mathcal {C}$ , because (i) it is proportional to the elastic energy in purely Hookean constitutive models of the polymers [Reference Beris and Edwards7] and (ii) it is equal to the sum of its principal stretches [Reference Sureshkumar, Beris and Handler52]. However, $\text {tr} \mathcal {C}$ does not provide a thorough enough description of polymer deformation. For instance, Beris and Edwards [Reference Beris and Edwards7] discovered that an increase in elasticity does not always result in an equivalent change in the mean velocity profile, since the mean stress deficit does not depend on any of the normal components of $\mathcal {C}$ . This illustration emphasizes the significance of simultaneously taking into account all of the elements of $\mathcal {C}$ to fully understand the polymer deformation and its impact on the velocity field. One way to acquire pertinent higher-order statistical descriptions of $\mathcal {C}$ is to use the fluctuating conformation tensor, $\mathcal {C}'$ (obtained by deducting the mean conformation tensor, $\mathcal {\overline {C}}$ , from the instantaneous tensor, $\mathcal {C}$ ). This fluctuating tensor, however, is not guaranteed to be physically realizable, because (i) whenever $\text {tr} \mathcal {C}'\leq 0$ , this implies negative material deformation and this tensor loses positive-definiteness, and (ii) equally probable states of contraction and expansion ( $\text {tr} \mathcal {C}'\in (0, 1)$ ) would be described by fluctuations of very different magnitudes. Because the logarithm of a positive definite matrix is a symmetric matrix and the collection of symmetric matrices forms a vector space [Reference Bhatia8], using $\log \mathcal {C}$ to assess fluctuations in $\mathcal {C}$ is more appealing. The Riemannian manifold of positive definite, polymer conformation tensors forces us to use the following metric, which is constructed as follows (see [Reference Bhatia8] for more information):

(5.4) $$ \begin{align} d(\mathbf{X},\mathbf{Y}) = \bigg[ \sum\limits_{i = 1}^3 ( \log \sigma_i ( \mathbf{X}^{-1} \mathbf{Y} ) )^2 \bigg]^{1/2}, \end{align} $$

where $\sigma _i$ are the eigenvalues of the matrix $\mathbf {X}^{-1} \mathbf {Y}$ . Using this metric, we define three scalar measures used to quantify the fluctuating polymer conformation tensor, which are enumerated as follows.

1. (i) Scalar Invariant $\delta _1$ is the volume ratio of fluctuations, defined as

   (5.5) $$ \begin{align} \delta_1 = \log \bigg( \frac{\det \mathcal{C}}{\det \overline{\mathcal{C}}} \bigg). \end{align} $$
   When $\delta _1 = 0$ , the mean and the instantaneous conformation tensors have the same volume and when $\delta _1$ is negative (positive), the instantaneous conformation tensor has a smaller (larger) volume than the mean volume.

2. (ii) Scalar Invariant $\delta _2$ is the measure of the shortest path between $\overline {\mathcal {C}}$ and $\mathcal {C}$ , defined as

   (5.6) $$ \begin{align} \delta_2 = d^2(\overline{\mathcal{C}}, \mathcal{C}). \end{align} $$
   Using (5.4), we can easily derive that $d^2(\overline {\mathcal {C}}, \mathcal {C}) = d^2(\overline {\mathcal {C}}^{-1}, \mathcal {C}^{-1})$ , which implies that this squared geodesic treats both expansions and compressions identically.

3. (iii) Scalar Invariant $\delta _3$ is the squared geodesic distance between $\mathcal {C}$ and its closest isotropic tensor,

   (5.7) $$ \begin{align} \delta_3 = \delta_2 - \tfrac{1}{3} \delta^2_1. \end{align} $$
   Notice that $\delta _3 = 0$ if and only if $ \delta ^2_1 = 3 \delta _2$ , in which case, $\mathcal {C}$ reduces to an isotropic tensor.

Note that the metric in (5.4) along with the three scalar invariants are designed to bypass the caveat of the regular invariants (such as trace and determinant) of the polymer conformation tensor, namely, the fluctuating conformation tensor fields may not retain positive definiteness in their standard arithmetic decomposition and hence, may not retain their physical meaning.

Figure 6 Contours of instantaneous (a, b) volume ratio, $\delta _1$ , (c, d) shortest distance from the mean, $\delta _2$ , (e, f) anisotropy index, $\delta _3$ , for the Zimm’s model (left column) and the Rouse model (right column) at simulation time, $T=7.15$ . Other parameters are fixed at $\nu = 0.3,~Re=70,~We=20$ .

Figure 7 Contours of instantaneous volume ratio, $\delta _1$ , for the elastic stress dominated Zimm’s model (). Other parameters set at $Re = 70, We = 15$ (first row), $Re = 1000, We = 15$ (second row), $Re = 70, We = 20$ (third row) and $Re = 1000, We = 20$ (fourth row).

Figure 8 Contours of instantaneous volume ratio, $\delta _1$ , for the elastic stress dominated Rouse model (). Other parameters set at $Re = 70, We = 15$ (first row), $Re = 1000, We = 15$ (second row), $Re = 70, We = 20$ (third row) and $Re = 1000, We = 20$ (fourth row).

5.5. Numerical results

Next, the in silico studies for two specific cases of monomer diffusion in Rouse chain melts () [Reference Rouse37] and in Zimm chain solution () [Reference Zimm55] are reported. The Rouse model predicts that the viscoelastic properties of the polymer chain can be described by a generalized Maxwell model, where the elasticity is governed by a single relaxation time, which is independent of the number of Maxwell elements (or the so-called “submolecules”). In contrast, the Zimm’s model predicts the (“shear rate and polymer concentration independent”) viscosity of the polymer solution by calculating the hydrodynamic interaction of flexible polymers (an idea that was originally proposed by Kirkwood [Reference Kirkwood27]) by approximating the chains using a bead–spring setup.

The instantaneous principle invariant of the conformation tensor, $\text {tr} \mathcal {C}$ , as well as the new invariants, (5.5)–(5.7) for a specific set of parameters values, $\nu =0.3, Re=70, We=20$ , and simulation time, $T=7.15$ , are presented for both the Zimm’s case (left column) and the Rouse case (right column) in Figure 6. Figure 6(a) (as well as Figure 6(b)) shows the logarithmic volume ratio, $\delta _1$ . We find that in Figure 6(a), we have predominantly negative values, indicating that the instantaneous volume is smaller than the volume of the mean conformation. Further, we observe regions of very high values of $\delta _1$ interspersed with regions of very low values, especially near the wall. This observation is the result of the slow diffusion of polymers in subdiffusive flows since there is no direct mechanism for smoothening out these “elastic shocks” in the tensor field. Since these elastic shocks originate away from the flow inlet and outflow, we conclude that these structures predominantly form due to the flow interaction with the rigid walls.

The measure, $\delta _1$ , does not distinguish between volume-preserving deformations. For example, $\delta _1$ does not distinguish between $\mathcal {C}$ and $\det (\mathcal {C}) \mathcal {C}_1$ for any tensor $\mathcal {C}_1$ with a unit determinant. In particular, $\delta _1 = 0$ does not imply $\mathcal {C} = \overline {\mathcal {C}}$ . To identify regions where the instantaneous polymer conformation equals the mean conformation and quantify the deviation when it is not, we use the squared geodesic distance away from the mean along the Riemannian manifold, $\delta _2$ (Figures 6(c) and 6(d)). Figure 6(c) indicates that the conformation tensor field, $\mathcal {C}$ , is significantly far away from $\overline {\mathcal {C}}$ , near the wall. This deviation of $\delta _2$ , in the near-wall region, can be explained via the “memory effect”, previously observed in regular Oldroyd-B fluids [Reference Sircar and Bansal44]. Finally, Figure 6(e)–(f) show the instantaneous contours of the anisotropy index, $\delta _3$ . This index shows how close the shape of instantaneous conformation tensor is to the shape of the mean conformation tensor, irrespective of volumetric changes. The visual resemblance of $\delta _2$ and $\delta _3$ suggests that deformations to the mean conformation are largely anisotropic near the wall. Next, we limit our focus on the contours of the first invariant, $\delta _1$ (after noting in Figure 6 that the contours of the other invariants are qualitatively similar) and for two different values of $\nu $ : $\nu =0.3$ (elastic stress dominated case) and $\nu =0.5$ (viscous stress dominated case).

5.5.2. Viscous stress dominated case: $\nu \ge 0.5$

Observe that the “flow structures” for the Zimm’s model in the elastic stress dominated case ( $\nu =0.5$ , Figure 9) are larger in size as well as in magnitude, in comparison with the Rouse model (Figure 10), at equal simulation times. Even within the respective models, we find that the macrostructures are more prominent (both in size and magnitude) at higher values of elastic relaxation (that is, at $We=20$ , third and fourth rows in Figures 9 and 10). Again, note that the structure formation is conspicuously absent at larger values of $Re$ (the second and the fourth rows in Figures 9 and 10).

Figure 9 Contours of instantaneous volume ratio, $\delta _1$ , for the viscous stress dominated Zimm’s model (). Other parameters set at $Re = 70, We = 15$ (first row), $Re = 1000, We = 15$ (second row), $Re = 70, We = 20$ (third row) and $Re = 1000, We = 20$ (fourth row).

Figure 10 Contours of instantaneous volume ratio, $\delta _1$ , for the viscous stress dominated Rouse model (). Other parameters set at $Re = 70, We = 15$ (first row), $Re = 1000, We = 15$ (second row), $Re = 70, We = 20$ (third row) and $Re = 1000, We = 20$ (fourth row).

We summarize our discussion by noting that for the selected values of the parameters, $\nu , \alpha $ , $Re$ and $We$ : (a) the elastic stress dominated case ( $\nu =0.3$ case) is comparatively more unstable than the viscous stress dominated case; (b) the Rouse model is comparatively more unstable than the Zimm’s model at low fluid inertia (or lower values of $Re$ ) and higher elastic relaxation (or higher values of $We$ ); and (c) temporal stability is achieved at higher values of $Re$ , irrespective of the model or the polymer concentration. Thus, our in silico studies not only corroborate the existence of the previously established temporally stable region at high fluid inertia [Reference Chauhan, Bansal and Sircar11], but also highlight the potential of fractional partial differential equations in effectively capturing the flow-instability transition in subdiffusive flows.