IMPLICIT–EXPLICIT TIME INTEGRATION METHOD FOR FRACTIONAL ADVECTION–DIFFUSION-REACTION EQUATIONS

Abstract

We propose a novel time-asymptotically stable, implicit–explicit, adaptive, time integration method (denoted by the $\theta $-method) for the solution of the fractional advection–diffusion-reaction (FADR) equations. The spectral analysis of the method (involving the group velocity and the phase speed) indicates a region of favourable dispersion for a limited range of Péclet number. The numerical inversion of the coefficient matrix is avoided by exploiting the sparse structure of the matrix in the iterative solver for the Poisson equation. The accuracy and the efficacy of the method is benchmarked using (a) the two-dimensional fractional diffusion equation, originally proposed by researchers earlier, and (b) the incompressible, subdiffusive dynamics of a planar viscoelastic channel flow of the Rouse chain melts (FADR equation with fractional time-derivative of order ) and the Zimm chain solution (). Numerical simulations of the viscoelastic channel flow effectively capture the nonhomogeneous regions of high viscosity at low fluid inertia (or the so-called “spatiotemporal macrostructures”), experimentally observed in the flow-instability transition of subdiffusive flows.

1. Introduction

Fractional partial differential equations (FPDEs) have emerged as a powerful tool for modelling the multiphysics and multiscale processes in numerical simulations ranging from physics [20–22] and biology [30] to quantitative finance [13]. For example, some of the most significant and profoundly published experimental results in random flow environments, such as the cytosol and the plasma membrane of biological cells [38], crowded complex fluids and polymer solutions [32], dense colloidal suspensions [29], and single-file diffusion in colloidal systems [28], are better rationalized within the fractional calculus framework. The growing number of applications of fractional derivatives in various fields of science and engineering indicate that there is a significant demand for better numerical algorithms. However, such algorithms are predominantly designed for one-dimensional (1D) problems [34], due to the severe memory restrictions imposed by these derivatives, a challenge that we alleviate in Section 2.2.

The individual physics or scale components in FPDE have very different properties that are reflected in their discretization, for example, in fractional advection diffusion reaction (FADR) systems [23], the discrete advection has relatively slow (or “nonstiff”) dynamics while the diffusion is a fast (or “stiff”) evolving process. Implicit–explicit (IMEX) integrators have been proposed as an attractive alternative (compared with the fully explicit or fully implicit time integration methods), where one combines the explicit (implicit) integration for the slow (fast) scale [14]. In an IMEX scheme, the system of equations assumes the form [2],

(1.1) $$ \begin{align} \frac{\partial^\alpha \Phi}{\partial t^\alpha} = \mathbf{f}(\Phi) + \eta \mathbf{g}(\Phi), \end{align} $$

where the superscript, $\alpha $ , denotes order of the fractional derivative and $\eta $ is a nonnegative parameter. In (1.1), the terms collected in $\mathbf {f}(\Phi )$ are on a slow time scale. Because they are (possibly) nonlinear, the implementation of a fully implicit integration scheme faces performance challenges, either from a poor performance of iterative solvers or from a complex Jacobian matrix associated with the problem. Therefore, it makes sense to solve this term explicitly. The terms in $\mathbf {g}(\Phi )$ , however, are on a fast time scale, and their explicit solution may require excessively small time steps to maintain numerical stability. Being usually linear in nature, they can be solved implicitly without further complications. The stability of such methods is still bounded by the Courant–Friedrichs–Lewy (CFL) condition [31], but because the fast terms are treated implicitly, these conditions are less strict when compared with a fully explicit scheme of similar formal order of accuracy. A class of such a time-asymptotically stable IMEX method is introduced in Section 2.

We have limited our focus in this work on the development, analysis and applicability of a novel spatiotemporal discretization method for the numerical solution of the 1D and two-dimensional (2D) FADR systems. To the best knowledge of the authors, a detailed analysis of the numerical methods for the FADR equation, in the finite-difference framework, does not exist. Such analysis is available for the fraction-diffusion equation [10] and the advection–diffusion equation [1]. Section 2 describes this method. Using the 1D linear FADR equation, the time-asymptotic stability analysis and the spectral analysis for the method is outlined in Section 3. In Section 4, the numerical method is benchmarked using the test cases for the 2D fraction diffusion equation, which was proposed by Brunner et al. [10]. The numerical results of the subdiffusive dynamics of the viscoelastic channel flow are delineated in Section 5. Section 6 concludes with a brief discussion of the implication of these results in future studies.

2. Methodology

Let $\chi $ be a bounded domain in $\mathbb {R}^2$ with a sufficiently smooth boundary $\partial \chi $ . We present the numerical method for an initial-boundary value problem for the time-dependent FADR equation with fractional time-derivative of order $\alpha \in (0,\,1)$ , as follows:

(2.1) $$ \begin{align} &\frac{\partial^\alpha u(\mathbf{x}, t)}{\partial t^\alpha} + K_1(\mathbf{x}, t) \nabla \cdot u(\mathbf{x}, t) = K_2(\mathbf{x}, t) \nabla^2 u(\mathbf{x}, t) + f(\mathbf{x}, t), \quad \mathbf{x} \in \chi, \quad t \in (0, T), \nonumber \\ & u(\mathbf{x}, 0) = u_0(\mathbf{x}), \quad \mathbf{x} \in \chi, \nonumber \\ & \alpha_1 u(\mathbf{x},t) + \alpha_2 \frac{\partial u(\mathbf{x}, t)}{\partial n} = g_0(\mathbf{x},t), \quad \mathbf{x} \in \partial \chi, \quad t \in (0, T), \end{align} $$

where $({\partial ^\alpha u(\mathbf {x}, t)})/({\partial t^\alpha })$ denotes the Caputo fractional derivative of the variable, $u(\mathbf {{x}},t)$ , of order $\alpha $ [35] with respect to t defined by

(2.2) $$ \begin{align} \frac{\partial^\alpha u(\mathbf{x}, t)}{\partial t^\alpha} = \frac{1}{\Gamma(1-\alpha)} \int^t_0 \frac{d t'}{(t-t')^\alpha} \frac{\partial u(\mathbf{x}, t')}{\partial t'}, \quad 0 < \alpha < 1, \end{align} $$

and the operators $\nabla (\cdot )$ and $\nabla ^2(\cdot )$ in (2.1) are the (integer order) gradient and the Laplacian operator in $\mathbb {R}^2$ , respectively. Here, $K_1(\mathbf {x}, t), K_2(\mathbf {x}, t)$ and $g_0(\mathbf {x}, t)$ are continuous real valued functions of $(\mathbf {x}, t)$ , and $\alpha _1, \alpha _2$ are constants in $\mathbb {R}$ . Additionally, $\Gamma (\cdot )$ is the standard $\Gamma $ -function. The FADR equation is related to the non-Markovian continuous-time random walk and is a model for anomalous diffusional flow-fields such as polymer melts [4–6, 11, 12], flows of liquid crystals [33, 46, 49, 50] as well as biological flows including mucus [42, 45, 47, 48, 51] and cartilage [43]. In the next three sections, we propose the numerical method for the spatiotemporal discretization of (2.1).

2.1. Time integration

The Caputo fractional time derivative in (2.1)–(2.2) is discretized using the difference approximation and studied by Podlubny [35]. 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 [41], 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. [10] as well as the fully implicit method recently proposed by Jannelli [23].

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” [16], 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) [36] 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 [41]. 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. Time-asymptotic stability and spectral analysis: linear 1D FADR equation

The linear 1D FADR equation (2.1) (with $K_1=c, K_2=\gamma , f(x, t) = \lambda u(x, t)$ for ${x \in (-\infty , \infty )}$ , where $c, \gamma $ and $\lambda $ are constants specifying the advection speed, coefficient of diffusion and coefficient of reaction, respectively) serves as a model for the FPDE replicating multiscale processes.

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 [39]. 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.

Theorem 1. The approximate solution to the 1D linear FADR equation using the $\theta $ -method in (3.3) with $\alpha \in (0,\,1)$ , on the finite domain $x \in [-L,\, L]$ with periodic boundary conditions, is time-asymptotically stable for all $t \ge 0$ .

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 [40]. 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 [41]. 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 [41]. 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 [24, 25].

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 [54]. 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) [40].

4. Method validation: 2D fractional diffusion equation

Using the spectrally relevant parameter values discussed in Section 3.2, the $\theta $ -method is verified by solving the 2D fractional diffusion equation ( $K_1 = f = 0$ and $K_2 = 1$ in (2.1)), originally proposed by Brunner et al. [10], for the case of (a) zero Dirichlet boundary conditions $(\alpha _1 = 1, \alpha _2 = 0,\, g_0(\mathbf {x},t) = 0)$ and (b) zero Neumann conditions $(\alpha _1 = 0, \alpha _2 = 1,\, g_0(\mathbf {x},t) = 0)$ , over the square domain, $(x, y) \in \Gamma = [-1,\,1]^2$ . The initial condition for the two test cases can be outlined as follows:

1. (i) Dirichlet case: $u_0(x, y) = \cos ({\pi }/{2} x ) \cos ({\pi }/{2} y )$ ;

2. (ii) Neumann case: $u_0(x, y) = \sin ({\pi }/{2} x ) \sin ({\pi }/{2} y )$ ,

and the exact solution for both cases is given by

$$ \begin{align*} u_{\text{Exact}}(x, y, t) = E_\alpha \big({-}\tfrac{1}{4} \pi^2 t^\alpha \big) u_0(x, y), \end{align*} $$

where $E_\alpha (z) = \sum ^\infty _{k=0} ({z^k}/{\Gamma (\alpha k + 1)}), \alpha> 0$ , is the one-parameter Mittag–Leffler function [15].

The domain, $\Gamma \cup \partial \Gamma $ , is discretized using $51 \times 51$ points and the relative error in $l_2$ -norm, $({\| \widetilde {u} - u_{\text {Exact}}\|_2}/{\| u_{\text {Exact}} \|_2})$ , versus discretization time, $\Delta t$ , is shown in Figure 5. The total computation time is fixed at $T=0.35$ for all cases. We note that slope of the error curves at $\alpha =1.0$ is one (highlighted with solid curves in Figure 5(a),(b)). This observation can be explained since at this value of $\alpha $ , the $\theta $ -method reduces to the standard, integer order, Euler method, which is first-order accurate. Further, note that the slope of the error curves for $\alpha < 1.0$ is sub-linear, indicating that for the fractional-order diffusion equation, the $\theta $ -method has sub-linear order of accuracy.

Figure 5 Relative error for the solution of the 2D fractional diffusion equation at simulation time $T=0.35$ , at $\alpha =1.0, \theta =1.0$ ( $\star $ ); $\alpha =0.9, \theta =1.0$ ( $\square $ ); $\alpha =0.67, \theta =1.0$ ( $\circ $ ); $\alpha =0.5, \theta =1.0$ ( $\triangle $ ); $\alpha =1.0, \theta =0.6$ ( $\ast $ ); $\alpha =0.9, \theta =0.6$ ( $\diamond $ ); $\alpha =0.67, \theta =0.6$ ( $\bullet $ ); $\alpha =0.5, \theta =0.6$ ( $+$ ); and for (a) Dirichlet boundary conditions and (b) Neumann boundary conditions.

5. Numerical simulation: 2D fractional viscoelastic channel flow

Next, we describe the incompressible, subdiffusive dynamics of a planar viscoelastic channel flow for polymer melts and present the numerical solution using the $\theta $ -method described in Section 2. Using the following scales for nondimensionalizing the governing equations: the height of the channel H for length, the timescale T corresponding to maximum base flow velocity (that is, $T = (H / \mathcal {U}_0)^{1/\alpha }$ ) for time and $\rho \mathcal {U}^2_0$ for pressure and stresses (where $\rho , \mathcal {U}_0$ is the density and the velocity scale, respectively), we summarize the model in streamfunction-vorticity formulation as follows:

(5.1a) $$ \begin{align} &Re \bigg[ \frac{\partial^\alpha \Omega}{\partial t^\alpha} + {\mathbf{v}} \cdot \nabla \Omega \bigg] = \nu \nabla^2 \Omega + \frac{(1-\nu)}{We} \nabla \times \nabla \cdot \mathcal{C}, \end{align} $$

(5.1b) $$ \begin{align} & \nabla^2 \psi = -\Omega, \end{align} $$

(5.1c) $$ \begin{align} &\frac{\partial^\alpha \mathcal{C}}{\partial t^\alpha} + {\mathbf{v}}\cdot \nabla \mathcal{C} - (\nabla {\mathbf{v}})^T\mathcal{C} -\mathcal{C} \nabla {\mathbf{v}} = \frac{{\mathbf{I}} - \mathcal{C}}{We}, \end{align} $$

where the variables $t, \psi , \mathbf {v} = (u, v) = ({\partial \psi }/{\partial y}, -{\partial \psi }/{\partial x}),\ \mathcal {C}, \Omega = \nabla \times \mathbf {v}$ denote time, streamfunction, velocity, polymer conformation tensor and vorticity, respectively. The parameters $\eta _0 = \eta _s + \eta _p$ and $\nu = \eta _s / \eta _0$ are the total viscosity and the viscous contribution to the total viscosity of the fluid, respectively. The dimensionless groups characterizing inertia and elasticity are Reynolds number, $Re = {\rho \mathcal {U}_0 H}/{\eta _0}$ , and Weissenberg number, $We = {\lambda ^\alpha \mathcal {U}_0}/{H}$ , respectively. The parameter, $\lambda ^\alpha $ , is the polymer relaxation time. Note that (5.1c) represents a fractional version of the regular Oldroyd-B model for viscoelastic fluids [44], derived in a recently reported work [11]. Two specific cases of the fractional derivative are considered, namely, monomer diffusion in Rouse chain melts ( $\alpha = {1}/{2}$ ) [37] and in Zimm chain solution ( $\alpha ={2}/{3}$ ) [55].

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 [3]. 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 [52]. The special treatment is an efficient transformation of the slope-limiting approach of Vaithianathan et al. [53], 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. [18], 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 [18, 26]. 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 [9]. 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 [7] and (ii) it is equal to the sum of its principal stretches [52]. However, $\text {tr} \mathcal {C}$ does not provide a thorough enough description of polymer deformation. For instance, Beris and Edwards [7] 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 [8], 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 [8] 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 () [37] and in Zimm chain solution () [55] 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 [27]) 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 [44]. 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.1. Elastic stress dominated case: $\nu < 0.5$

Flows with parameter values, $\nu < 0.5$ (or the elastic stress dominated case), are those associated with higher concentration of polymers per unit volume. The contours of the first invariant, $\delta _1$ , for the Zimm’s model and the Rouse model are shown in Figures 7 and 8, respectively, at equal snapshots. At lower values of elastic relaxation (represented by the parameter, $We=15$ ), the flow-macrostructure evolution in the Rouse model is slow, in comparison with the Zimm’s model (that is, comparing the contour values in the first two rows in Figure 7 versus the first two rows in Figure 8). A reversal of this trend occurs at higher values of elasticity (that is, the third and the fourth rows in Figures 7, 8 at $We=20$ ). We attribute this observation due to the following reasoning. The Rouse model represents a“thicker” fluid, or fluids with slower diffusion than the Zimm’s solution, due to the smaller fractional time derivative. Flows with smaller time derivative (or the thicker polymer melt case) are those associated with higher concentration of polymers per unit volume [9]. Physically, the formation of these “spatiotemporal macrostructures” are associated with the entanglement of the polymer chains at microscale [38], leading to localized, nonhomogeneous regions with higher viscosity. Experiments [19] have shown that non-Newtonian fluids with a higher polymer concentration have a greater tendency for the polymer strands to agglomerate (especially at higher elastic relaxation), or the so-called “over-crowding effect” [17].

Further, notice the disappearance of the structures at larger values of $Re$ (that is, notice the time snippets at $Re=1000$ , the second and the fourth rows in Figures 7 and 8). This observation can be attributed to the fact that the macrostructures are “washed out” of the channel at higher flow velocities. In an earlier work, we had indicated the existence of “temporally stable regions at high fluid inertia” through a spatiotemporal linear stability analysis [11]. In this work, we associate temporally stable regions with regions devoid of flow structures.

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 [11], but also highlight the potential of fractional partial differential equations in effectively capturing the flow-instability transition in subdiffusive flows.

6. Conclusion

This investigation addresses the development as well as the time-asymptotic and spectral stability of a novel class of numerical method for the spatiotemporal discretization of FPDE. Section 2 presented the method, including the time integration (Section 2.1) and the spatial discretization (Section 2.3). Using the 1D linear FADR equation, the time-asymptotic stability and spectral analysis were outlined in Section 3. The method was validated using the test cases for the 2D fraction diffusion equation in Section 4. Section 5 described the numerical results of the subdiffusive dynamics of the viscoelastic channel flow. We conclude our discussion with a remark that the focus of this present work was on the development of the numerical method and not on the physical description of the subdiffusive flow dynamics. Hence, a comprehensive study on the mechanics of subdiffusive channel flow, using the numerical method developed in this work, is reported elsewhere ¹ .