An insight into the biomaterials used in craniofacial tissue engineering inclusive of regenerative dentistry

1.   Introduction

Fractional calculus, the study of derivatives and integrals of non-integer orders, has increasingly become a powerful tool in modeling physical phenomena that defy conventional descriptions under ideal conditions ^[51]. It addresses the limitations inherent in classical models by providing a more accurate representation of systems with anomalous properties. Despite a relatively slow historical progression, the past decade has seen a renaissance, with fractional-order models advancing alongside traditional integer-order counterparts.

These fractional-order models, encompassing both ordinary and partial differential equations with non-integer order terms have been instrumental in advancing our understanding of a diverse array of physical phenomena. They have been successfully applied in fields ranging from signal processing and control theory to more complex systems in diffusion processes, thermodynamics, biophysics, and blood-flow dynamics, as well as in electrodynamic, electrochemical, and electromagnetic theories, not to mention their significant roles in continuum and statistical mechanics, along with dynamical systems, see References^[5,12,48,56].

Substantial progress in the numerical approximation of fractional-order models have been made, enhancing the precision and applicability of these models in computational simulations. Beyond the traditional finite difference schemes for fractional equations ^[11,54], noteworthy advancements have been made in the Finite Element Methods (FEM). For instance, the work by Acosta and Borthagaray ^[2] shed light on regularity results for the analytic solution of the fractional Poisson problem and delineated convergence rates for FEM approximations, with these findings further substantiated by computational examples in ^[1]. Additionally, Ervin and Roop ^[32] contributed to this field by employing the FEM approach to steady-state fractional advection dispersion equations,

where $D$ represents the first-order spatial derivative, $_{0}D_{x}^{-\beta}$, $_{x}D_{1}^{-\beta}$ represent the left and right fractional integral operators with $0\leq\beta < 1$ and $0\leq p$, $q\leq 1$ satisfying $p+q = 1$. The FEM with Galerkin framework was used in ^[33] to solve fractional diffusion wave equations. Jin et al. described the Petrov-Galerkin FEM for the fractional convection-diffusion equations ^[35], where they considered the 1D fractional boundary value problem

where $f\in L^{2}(D)$ and $_{0}D_{x}^{\alpha}u$ denotes Caputo fractional derivative of order $\alpha\in({3}/{2}, 2)$.

A comprehensive review of numerical approximations for fractional models is available in ^[30] and the references cited therein. Recently, focussing on the polygonal approximation techniques, Zhang et al. ^[58] studied a local projection stabilization Virtual Element Method (VEM) for the time-fractional Burger's equation with higher Reynold's number. Recent attempts with the VEM approximation for time-fractional Partial Differential Equations (PDEs) ^[31,57,58] highlight the need to further develop the studies to more general nonlinear problems.

Our research extends the VEM to such a time-fractional nonlinear partial differential equation, providing both theoretical analysis and validation through numerical experiments. In this work, we consider the 2D time-fractional nonlinear convection-diffusion equation

where $\Omega\subset\mathbb{R}^2$ denotes a bounded polygon region with boundary $\partial\Omega$ and the order of time derivative varies such that $0 < \alpha < 1$. The convective field is represented by ${ \boldsymbol{b}(\boldsymbol{x}) = (b_{1}(\boldsymbol{x}), b_{2}(\boldsymbol{x}))^T}$, the function $f(u)$ encapsulates the nonlinear term or terms of interest, while $g(\boldsymbol{x}, t)$ is a forcing term. To show theoretical estimates, suitable regularities of the terms in (1.1) will be discussed in Section 3. According to ^[38], notation ${}^{c}_{0}D^{\alpha}_{t}$ denotes the Caputo fractional derivative of order $\alpha$ and is defined as,

Background material on the VEM. The VEM is used for approximation of the solution of PDEs on general polygonal and polyhedral meshes. It is an extension of the finite element method and shares some similarities with it. Like the FEM, VEM also formulates the problem in a weak form, which involves multiplying the PDE by a test function and integrating by parts over the domain. The weak form allows for the use of piecewise smooth approximations and leads to a system of linear equations that can be solved numerically. Differently from FEM, the VEM introduces the concept of "virtual element spaces", which are sets of functions that approximate the exact solution of the PDEs at both the local (elementwise) and global level. Such functions are implicitly defined as the solution to a local PDE problem, and are never explicitly used to compute the coefficients of the global matrix. The global spaces are defined on polygonal or polyhedral meshes and are constructed to be compatible with the underlying mesh structure and possess some degree of global regularity. In depth examinations of the parallelism between the VEM and finite elements for polygonal/polyhedral meshes are provided in ^[25,45], and its relationship with BEM-influenced FEM approaches in ^[24]. The VEM is also a variational reinterpretation of the Mimetic Finite Difference (MFD) approach, specifically its nodal variant introduced in ^{[15,17,22,42]}. An extensive discussion on MFD is found in the comprehensive review article ^[41] and the book ^[20]. A major advantage of the VEM is its ability to handle arbitrary polygonal and polyhedral meshes. Unlike the "traditional" FEM, which typically requires triangular or quadrilateral elements in 2D and tetrahedral or hexahedral elements in 3D, the VEM can handle meshes composed of elements with any number of sides with an upper bound from theoretical level. This flexibility is particularly useful when dealing with complex geometries. The VEM is designed to provide stable and accurate solutions for a wide range of problems. The method incorporates stabilization techniques to handle challenges associated with irregular meshes and distorted elements. Additionally, the VEM can achieve high-order accuracy by incorporating polynomial approximation spaces of higher degree in the virtual element space definition. Initially introduced to address elliptic problems such as the Poisson equation, cf. ^[6,19], the VEM was later applied to linear and nonlinear diffusion, convection-diffusion, and convection-reaction-diffusion problems ^[3,4,10,14] the Stokes equations ^{[13,16,43,44]}, the polyharmonic equation ^[9], and many other models. Meanwhile, the nonconforming formulation for diffusion problems was proposed in ^[29] as the finite element reformulation of ^[40] and later extended to general elliptic problems ^[26], Stokes problem ^[23], and the biharmonic equation ^[8,59]. A recent contributed book documents the state of the art of this methodology, cf. ^[7].

Paper's outline. The structure of this paper is as follows. Section 2 outlines some preliminary and established results, the model problem with weak formulation framework is given in Section 3. Section 4 presents the discrete formulation, well-posedness and convergence analysis of a fully discrete VEM. Section 5 offers numerical examples that substantiate our theoretical findings. Section 6 concludes the paper with a summary of our work.

2.   Definitions and notations

The concept of fractional derivatives has been explored and delineated extensively within the mathematical literature. In this section, we recall some definitions and background results that will be useful in the paper. For a comprehensive exposition of these concepts, the reader is directed to refer ^[49,50].

2.1.   Riemann-Liouville fractional derivatives

Let $\alpha$ be a positive real number, $m$ a non-negative integer, and $f(x)$ a univariate function defined on $x\in[a, b]$. The left and right hand Riemann-Liouville (R-L) fractional derivatives of $f(x)$ of order $\alpha$, with $m-1 < \alpha\leq m$, are defined as

where $\Gamma$ denotes, as usual, the Gamma function ^[53].

2.2.   Caputo fractional derivatives

At the end of the sixties, Caputo introduced the formulation for fractional derivatives, which jointly with Minardi was then applied in their research on viscoelasticity theory ^[27]. Caputo's definition of the fractional derivative is given by

where $I^{m-\alpha}$ is a corresponding integral of $m-\alpha$ fractional order. This derivative can be interpreted as:

2.3.   Notation

On a domain $\Omega$ and $v\in L^2(\Omega)$, the $L^2$ inner product is defined by $\left < v, v\right > = \underset{\Omega}{\int}v(x)v(x) dx$ and the associated $L^2(\Omega)$ norm is defined by $||v|| = (\underset{\Omega}{\int}|v(x)|^2dx)^{1/2}$. For any non-negative integer $n$, $H^n(\Omega)$ denotes a Sobolev space with the related norm

and seminorm $|\cdot|_{n}$. Let $C_{0}^{\infty}(\Omega)$ be a space of infinitely differentiable functions with compact support in $\Omega$ and $H_{0}^{m}(\Omega)$ is a closure of $C_{0}^{\infty}(\Omega)$ with respect to norm $||\cdot||_{m}$. Finally, we denote by $\mathbb{P}_{k}(\Omega)$ the space of all polynomials over $\Omega$ of degree up to $k$. Also, $C$ is independent of mesh sizes, denotes a generic constant that may vary at different occurrences.

3.   Model problem

We consider the following time-fractional nonlinear convection-diffusion-reaction equation:

Here, $0 < \alpha < 1$, $^{c}_{0}D^{\alpha}_{t}u$ denotes the $\alpha$-th order Caputo fractional derivative of $u$, and $g(\boldsymbol{x}, t)\in L^2(\Omega)$ is a forcing/load term. For the purpose of theoretical analysis, we make the assumptions:

(H1) Let $f: \mathbb{R}\to\mathbb{R}$ be a Lipschitz continuous function, that is, there exists $L > 0$ such that,

(H2) We also suppose that $\nabla\cdot\boldsymbol{b}\in L^{\infty}(\Omega)$ and there exists $\mu_{0}\geq 0$ such that for almost every $\boldsymbol{x}\in \Omega$,

Then, the following wellposedness theorem states the existence and uniqueness of the solution to problem (3.1). The proof of these results can be found in ^[36].

Theorem 1 (Wellposedness). Under Assumption (H1), problem (3.1) admits a unique solution $u$ for $0 < \alpha < 1$, such that,

where $C$ is a strictly positive, real constant (independent of $u$).

Li et al. ^[39] established a discrete fractional Grönwall type inequality for Caputo type fractional derivative. Taking the motivation from this work, we here establish a discrete fractional Grönwall type inequality (Grunwald-Letnikov approximation) to R-L type fractional derivative. To do so, we consider the relation between Caputo and R-L fractional derivatives expressed as

Since $u(\boldsymbol{x}, 0) = 0$, we find that $^{c}_{0}D^{\alpha}_{t}u(\cdot, t) = {^{R}_{0}D^{\alpha}_{t}}u(\cdot, t)$. As our initial model has the Caputo type derivative, it is important to take zero initial conditions to use the Caputo and R-L type interchangeably because they are equal when initial conditions are zero. It is noteworthy that the methodology (Grünwald equality) for a time-fractional derivative used is specific for the R-L type ^[35,37]. The reason to start with the Caputo-type derivative is because of its versatile and meaningful interpretation of physical phenomena and allows for clear formulations of initial and boundary conditions. Furthermore, the literature indicates that there are no inherent limitations on the initial conditions, as the Grünwald approximation can be extended to the Caputo-type derivative with inhomogeneous initial values. This extension is achieved by adding appropriate correction terms, as referenced in the work by Scherer et al. ^[52]. Now, we write the model problem (3.1) in terms of the R-L fractional derivative reads as

The weak/variational formulation of the model fractional elliptic problem (3.2) is:

Find $u(t)\in H_{0}^{1}(\Omega)$ for almost all $t\in (0, T)$ such that

where the bilinear forms $m(\cdot, \cdot)$, $a(\cdot, \cdot)$ and $b(\cdot, \cdot)$ are defined as:

and $\left < \cdot, \cdot\right >$ is the duality product in $L^2(\Omega)$. The wellposedness of the weak form can be proved along similar lines of the work by Jin et al.^[35,36] for the time-fractional nonlinear subdiffusion equation.

4.   Virtual element discretization

Let $\{{\mathcal{T}_{h}}\}_{h > 0}$ be a family of polygonal meshes covering $\Omega$, each labelled by the mesh size parameter $h$ which is the maximum of the element diameters $h_E$ associated with the mesh ${\mathcal{T}_{h}}$, i.e., $h: = \underset{E\in{\mathcal{T}_{h}}}{\text{max}}\; h_{E}$. For mesh regularity, we assume there exists $\rho > 0$ such that for every $E\in{\mathcal{T}_{h}}$ :

● $E$ is star shaped for a ball of radius $\geq\rho h_{E}$,

● distance between any two vertices of $E$ is $\geq\rho h_{E}$.

Referring to ^[6], for $k\geq\mathbb{N}$, the local virtual element space $V_{h}^{E}$ of order $k$ is defined as,

where $\partial E$ is the polygonal boundary of $E$ and $e\subset\partial E$ denotes a generic edge.

Next, we define some useful polynomial projection operators. We consider the elliptic projection operator ${\Pi_{k}^{\nabla, E}:H^{1}(E)\to\mathbb{P}_{k}(E)}$, which is the orthogonal projector onto the space of polynomials of degree $k$ with respect to the $H^1$ semi-norm. Formally, for every $v\in H^{1}(E)$, the $k$-degree polynomial ${\Pi_{k}^{\nabla, E}v}$ is the solution to the variational problem

where the second condition is needed to fix the gradient kernel. We also define the $L^2$-orthogonal projection operator ${\Pi_{k}^{0, E}:H^{1}(E)\to\mathbb{P}_{k}(E)}$, such that for any $v\in H^1(E)$, the polynomial ${\Pi_{k}^{0, E}v}$ satisfies:

Using the elliptic projection operator, we define the "enhanced" local virtual element space denoted $W_{h}^{E}$ by

where $\mathbb{P}_k/\mathbb{P}_{k-2}(E)$ denotes the quotient space of equivalence classes of polynomials, where two polynomials are equivalent if their difference is a polynomial of degree at most $k-2$ (roughly speaking, in the implementation, we can consider the union of the linear spaces of homogeneous polynomials of degree exactly $k$ and $k-1$). Accordingly, we define the global virtual element space of order $k$ as:

A function $w_h\in W_{h}^{E}$ is uniquely described by the following sets of values that we can take as the local Degrees of Freedom (DoFs):

$\textbf{(D1)}$ For $k\geq1$, the vertex values $w_h(V_i)$, where $V_i$ are the vertices of element $E$;

$\textbf{(D2)}$ For $k > 1$, the values of $w_h$ at the $(k-1)$ internal Gauss-Lobatto quadrature points on each edge $e$ subset of $\partial E$;

$\textbf{(D3)}$ For $k > 1$, the polynomial moments up to order $k-2$,

The DoFs of the global virtual element space of order $k$ are given by collecting the local DoFs of the elemental spaces $W^{E}_{h}$ from all elements $E\in{\mathcal{T}_{h}}$ that are glued with $C^{0}$ continuity. The unisolvence of these degrees of freedom is proved in ^[6]. A crucial property of definition (4.1) is that the projections ${\Pi_{k}^{\nabla, E}}w_h$ and ${\Pi_{k}^{0, E}}w_h$ are computable using only the DoFs (D1)–(D3) of $w_h\in W_h^E$.

The semi-discrete VEM approximation of the weak formulation (3.3) is given as:

Find $u_{h}(t)\in W_h$, for almost all $t\in (0, T)$ such that

Here, the global discrete bilinear forms $a_{h}(\cdot, \cdot):W_{h}\times W_{h}\to \mathbb{R}, \; b_{h}(\cdot, \cdot):W_{h}\times W_{h}\to \mathbb{R}\; \text{and}\; m_{h}(\cdot, \cdot):W_{h}\times W_{h}\to \mathbb{R}$ are defined as:

where the $computable$ discrete local bilinear forms $a_{h}^{E}(\cdot, \cdot):W_{h}^{E}\times W_{h}^{E}\to\mathbb{R}$, $b_{h}^{E}(\cdot, \cdot):W_{h}^{E}\times W_{h}^{E}\to\mathbb{R}$, and $m_{h}^{E}(\cdot, \cdot):W_{h}^{E}\times W_{h}^{E}\to \mathbb{R}$ are defined over every $E\in\mathcal{T}_h$ as follows:

where the bilinear forms $s_{a}^{E}:{V_{h}^{E}}\times {V_{h}^{E}}\to \mathbb{R}$ and $s_{m}^{E}:{V_{h}^{E}}\times {V_{h}^{E}}\to \mathbb{R}$ are a locally admissible stabilizing terms that are symmetric positive definite and satisfy

for some positive constants $c_{0}, c_{1}, d_{0}\; \text{and}\; d_{1}$ independent of $h$ and $E$.

These stabilization terms ensure the local stability of the corresponding elemental bilinear forms $a_{h}^{E}(\cdot, \cdot)$ and $m_{h}^{E}(\cdot, \cdot)$ (a formal definition follows below). For the computational purpose, we consider the "dofi-dofi" stabilization, see ^[46]. For the sake of completeness, we refer to various specific choices of the $computable$ definitions for the stabilizers available in the literature ^{[18,19,21,28]}.

We recall that the bilinear forms $a_{h}^{E}(\cdot, \cdot)$ and $m_{h}^{E}(\cdot, \cdot)$ possess the fundamental properties of polynomial consistency and stability, that are useful for theoretical estimation.

$\textbf{Polynomial consistency:}$ $\forall\; E\in{\mathcal{T}_{h}}$ and for any $w_h\in W_{h}^{E}$, it holds that

$\textbf{Stability:}$ There exist two pairs of constants independent of $h$, say, $(\alpha_{*}, \alpha^{*})$ and $(\gamma_{*}, \gamma^{*})$ with $0 < \alpha_{*}\leq\alpha^{*}$ and $0 < \gamma_{*}\leq\gamma^{*}$, such that for all $w_h\in W_{h}^{E}$, the two following equivalence holds:

The conditions above ensure that the non-polynomial parts $s^E_a(\cdot, \cdot)$ and $s^E_m(\cdot, \cdot)$ scale as the polynomial parts of $a^E_h(\cdot, \cdot)$ and $m^E_h(\cdot, \cdot)$, respectively.

To compute the nonlinear forcing term, we use the $L^{2}$-orthogonal projector ${\Pi_{k}^{0, E}}$ previously defined. For each element $E$, we first define $f_{h}(u_{h})$ as

The orthogonality of ${\Pi_{k}^{0, E}}$ implies that

Hence, for $f(u_{h})\in L^{2}(\Omega)$ and $u_h\in W_h$, we state that

by introducing a global $L^2$-orthogonal projection operator ${\Pi_{k}^{0}}:L^2(\Omega)\to\mathbb{P}_{k}({\mathcal{T}_{h}})$ onto the piecewise polynomial space of degree $k$ built on the mesh partition ${\mathcal{T}_{h}}$. The global projection operator is such that

We outlined that this approximation of the forcing term is computable since the projection ${\Pi_{k}^{0, E}}u_{h}$ is computable from the degrees of freedom of $u_h$ in the enhanced space $W_h^E$.

4.1.   Fully-discrete VEM

To develop a fully discrete VEM, the next step is to approximate the fractional derivative in the temporal direction. Let $0 = t_0 < t_1 < \ldots < t_K = T$ be a uniform partition of the interval $(0, T]$ with time step size $\tau = T/K$ for some positive integer $K$. By employing the Grunwald-Letnikov approximation, we effectively calculate the R-L type fractional derivative as obtained in the work of Kumar et al. ^[37]. This computation yields

Here, the weights are given by

where $\boldsymbol{R_n}$ satisfies the estimate $||\boldsymbol{R_n}||\leq c\tau$, and the approximation (4.2) is also $O(\tau)$ accurate.

Denote

where $0\leq \theta\leq 1$. For any arbitrary choice of $\theta\in [0, 1]$ in the discrete scheme, the numerical results show no significant effect on the accuracy, independent of the fractional-order $\alpha\in (0, 1)$. Hence, for convenience and reduced simplified notation, we fix $\theta\, = \, \alpha$ in the following fully discrete VEM formulation:

Find $U_{h}^{n}\in W_h$, $n = 1, 2, \ldots, K$, such that

Using the definition of operator $^{R}D^{\alpha}_{t_n}$ from (4.2), we rewrite the equation above as:

This discrete formulation provides a system of algebraic equations that we solve to compute the virtual element approximation.

To prove the well-posedness of the fully discrete VEM (i.e., the uniqueness and existence of the virtual element approximation), we need the following proposition, which results from an application of the Brouwer's fixed point theorem, see, e.g., ^[55] for the details.

Proposition 1. Let $\mathcal{H}$ be a finite-dimensional Hilbert space with inner product $\left < \cdot, \cdot\right >$ and norm $||\cdot||$ and let $S:\mathcal{H}\to\mathcal{H}$ be a continuous map such that,

Then, there exists $w\in \mathcal{H}$ such that,

Applying Proposition 1, we can derive the wellposedness result that we state in the next theorem.

Theorem 2 (Wellposedness). Let $U_{h}^{0}, \cdots, U_{h}^{n-1}$ be the first $n$ virtual element fields solving the fully discrete VEM for some given initial and boundary condition. Then, a unique solution $U_{h}^{n}$ of the fully discrete VEM exists.

Proof. We rewrite Eq (4.3) as

and multiply each term by $\left(1-\frac{\alpha}{2}\right)$ to obtain

Equations (4.4) and (4.5) are equivalent, and their solutions coincide. Then, we define the continuous operator $G:W_h\to W_h$ such that

for $W^{n, \alpha}$ and $V$ in $W_h$. Setting $V = W^{n, \alpha}$, we find that

From Assumption (H1), we have that

from which it follows that

By using the Cauchy-Schwarz inequality in (4.7), Eq (4.8) and noting that $w_{j}^{(\alpha)} < 0$ for $1\leq j\leq n$, we find that

Since, $\tau^{\alpha}||W^{n, \alpha}|| > 0$, we rewrite (4.9) as

Then, $m(G(W^{n, \alpha}), W^{n, \alpha}) > 0$, if and only if it holds that

Choosing $\tau^{\alpha} < \frac{1}{\left(1-\frac{\alpha}{2}\right)a}$, there exists a function $W^{n, \alpha}$, such that

This inequality implies that $m(G(W^{n, \alpha}), W^{n, \alpha}) > 0$. Thus, for $||W^{n, \alpha}|| = \rho$, we have that

and Proposition 1 implies the existence of the virtual element approximation.

To prove the uniqueness of a virtual element function $U_{h}^{n, \alpha}$ solving problem (4.3), assume that $U_{h, 1}^{n, \alpha}$ and $U_{h, 2}^{n, \alpha}$ are two solutions. To ease the notation, we denote $U_1 = U_{h, 1}^{n, \alpha}$ and $U_2 = U_{h, 2}^{n, \alpha}$. Then, from (4.5) it follows that

Setting $v_h = U_1-U_2 = r$ in (4.10) and using Assumption (H1), we find that

and choosing $\tau^{\alpha} < \frac{1}{\left(1-\frac{\alpha}{2}\right)L}$ sufficiently small, with $L > 0$ being the Lipschitz constant, we find that $||r^2||\leq 0$, which implies the uniqueness of the solution. □

4.2.   Convergence and a priori error bounds

From the theoretical results that we prove in this subsection, some important aspects of the VEM and the fully discrete scheme (4.3) are achieved. These points are given as follows:

● The efficiency with which the VEM is combined with a finite difference scheme to form a fully discrete scheme.

● The dependence of force/load term and reaction term on both the spatial and temporal domains and extended VEM to such terms nested with both domain discretizations.

● The optimal convergence order for a fully discrete scheme is achieved with theoretical analysis and supported with numerical results on several sets of meshes, including non-convex meshes.

To prove the convergence of the VEM, we derive suitable a priori error bounds. To this end, we need the fractional Grönwall type inequality discussed here. The Grönwall weights $w_{i}^{(\alpha)}$ can be computed as

Let $g_{n}^{(\alpha)} = \sum_{i = 1}^{n}w_{i}^{(\alpha)}$, so that

Note that $w_{i}^{(\alpha)}$ satisfy

Therefore, $g_{i-1}^{(\alpha)} > g_{i}^{(\alpha)}\; \text{for}\; i = 1, \ldots, n$. Using such a definition of $g_{n}^{(\alpha)}$, we rewrite the R-L derivative expansion as

Since $u(\boldsymbol{x}, t_{0}) = 0$, we have that

where $\delta u(\boldsymbol{x}, t_{i}) = u(\boldsymbol{x}, t_{i})-u(\boldsymbol{x}, t_{i-1})$ for every $i = 1, \ldots, n$.

Furthermore, we consider the following technical lemmas that we will use in the next derivations. Their proof can be found in ^[39].

Lemma 1. Consider the sequence $\{\phi_n\}$ defined by

Then, $\{\phi_n\}$ satisfies

and the following inequalities:

for any integer $k\geq1$.

Lemma 2. Let $\big\{c^n:\, n\geq 0\big\}$ and $\big\{d^n:\, n\geq 0\big\}$ be two non-negative sequences and $\lambda_1$ and $\lambda_2$ be two non-negative constants. Then, for

a positive constant $\tau^*$ exists such that for $\tau\leq\tau^*$ we have that

where $E_{\alpha}(z) = \sum_{j = 0}^{\infty}\frac{z^j}{\Gamma(1+\alpha j)}$ is the Mittag-Leffler function ^[34] and $\lambda = \lambda_1+\lambda_2$.

Lemma 3. For any sequence $\big\{i^k\big\}_{k = 0}^{K}\subset W_h$, the following inequality holds

Proof. We rewrite Eq (4.16) as

Then, use the fact that, $\left(1-\frac{\alpha}{2}\right)w_{1}^{(\alpha)}+\frac{\alpha}{2}w_{0}^{(\alpha)} < 0$, and $w_{j}^{(\alpha)} < 0\; \forall j\geq 1$, and we obtain

which proves inequality (4.16). □

Now, we derive the a priori error bound for the virtual element approximation $U_{h}^{n}$, which we state in the following theorems.

Theorem 3 (A priori bound of discrete VEM). Let $U_{h}^{n}$ be the solution to the semi-discrete VEM. Then, there exists a positive constant $\tau^*$ such that for $\tau\leq\tau^*$ the virtual element field $U_{h}^{n}$ satisfies

where $C$ is a positive constant independent of $h$ and $\tau$.

Proof. The virtual element solution field $U_{h}^{n}$ must satisfy the variational condition

Setting $v_h = U_{h}^{n, \alpha}$, we obtain

Then, we use Assumption (H1) to find

where $C = \text{max}\{1, \frac{L}{2}\}$, $L$ being the Lipschitz constant. We apply the standard inequality $(a+b)^{2}\leq 2(a^{2}+b^{2})$, $a, b\geq0$, and we have that

Using the results of Lemma 3 yields

which implies that

Finally, in view of Lemma 2, we find a positive constant $\tau^*$ such that $\tau\leq \tau^*$ and $||U_{h}^{n}||^{2}\leq C$, which proves the theorem's assertion. Here $C$ is a positive constant independent of the mesh size $h$ and the time step size $\tau$, and it represents more specifically an upper bound on semi discrete VEM solution $U_{h}^{n}$. It depends on the assumption of Lipschitz continuity of function $f$ and Lipschitz constant $L$. □

The last result of this section is about the convergence of the fully discrete VEM, which we prove in the next theorem by deriving an a priori error estimate.

Theorem 4 (Convergence of the fully discrete VEM). Let $u$ and $U_{h}^{n}$ be the solution to problem (3.1) and the fully discrete VEM (4.3) under Assumption (H1) respectively, with convex regularity of all the elements $E\in\mathcal{T}_h$. Then, there exists a positive constant $\tau^*$ such that

for $\tau\leq\tau^{*}$, where $C$ is a positive constant independent of $\tau$ and $h$.

Proof. Let $\Pi_{h}$ be the projection operator that satisfies

According to Cangiani et al. ^[26] and the references therein, there exists a positive constant $C$ independent of $h$ such that

The regularity of the domain at $j = 0$ is implied from the convex regularity of the local elements. Now, using this projection operator, we rewrite the approximation error as

with the obvious definition of $\rho_{h}^{n}$ and $\Theta_{h}^{n}$.

Assuming that

for any arbitrary $v_h\in W_h$, we find that $\Theta_{h}^{n}$ satisfies

Then, by using (4.3) and (4.19) in (4.22), we obtain that

The weak form of (3.1) implies that

Using (4.24) in (4.23) yields

Setting $v_{h} = \Theta_{h}^{n, \alpha}$ in (4.25) and using Cauchy-Schwarz inequality, we obtain

We add and subtract $u^{n, \alpha}$ and $\Pi_hu^n$ to the argument of $||u^{n}-U_{h}^{n, \alpha}||$; then, we apply the triangular inequality twice and use the definition of $\rho_{h}^{n}$ and $\Theta_{h}^{n}$ from (4.21), and we obtain

We also note that

and

Using (4.27)–(4.29) in (4.26), we obtain

which gives

From (4.31), we can get

where $C^{*} = 3L+1$. Furthermore, we also can write

Applying Lemma 2, there exists a positive constant $\tau^*$ such that $||\Theta_{h}^{n}||^{2}\leq C{({\tau}+h^{k+1})}^{2}$ for $\tau\leq\tau^{*}$, which implies that $||\Theta_{h}^{n}||\leq C{({\tau}+h^{k+1})}$. The theorem's assertion follows on applying the triangular inequality along with (4.20) and the definition of $\Theta_{h}^{n}$ from (4.21). □

4.3.   Implementation details

A short note on the implementation is discussed in this subsection. We solve the resulting system of nonlinear algebraic equations from (4.3) by incorporating Newton's method.

Let $N$ be the dimension and $\{\phi_i\}_{ i = 1}^N$ be the canonical basis for the global virtual element space $W_h$. For some $\beta_{i}^{n}\in \mathbb{R}, \; i = 1, 2, \ldots, N$, we can write the solution of $U_{h}^{n}\in W_h$ of (4.3) as

After defining $\boldsymbol{\beta}^{n}: = (\beta_{1}^{n}, \beta_{2}^{n}, \ldots, \beta_{N}^{n})^T$, using the value of $U_{h}^{n}$ from (4.33) in (4.3), we get the following nonlinear algebraic equation

where

Using the Newton's method in (4.34), we obtain the matrix system

where $\boldsymbol{H} = (H_1, H_2, \ldots, H_N)^T$ and the entries of the $(N\times N)$ Jacobian matrix $\boldsymbol{J}$ are given by

where $1\leq i, l\leq N$.

5.   Numerical experiments

In this section, we assess the performance of our approximation method by solving problem (3.1) on three different test cases and by using a sequence of refined uniform square, regular Voronoi meshes, and non-convex meshes (sample mesh configurations given in Figure 1). There are a few reasons to choose the different sets of meshes, including the non-convex meshes.

● To check the efficiency of the constructed VEM as it is mentioned that the VEM is a generalization of FEM over polygonal meshes.

● The use of square meshes was to provide an explanation that VEM can adjust the triangular or quadrilateral elements as well and can provide the accuracy of classical FEMs on these meshes.

● The use of non-convex meshes provides the explanation of the point that in VEM, the basis functions are constructed virtually as compared to the implicit definition of basis function in classical FEMs.

To this end, we measure the approximation errors as the difference between the exact solution $u$, the $L^2$ orthogonal projection $\Pi_{k}^{0}u_{h}$, and the elliptic projection $\Pi_{k}^{\nabla}u_{h}$ of the virtual element approximation $u_h$ through the formulas

By comparing the errors at two subsequent mesh refinements, we compute the convergence rate. From the theoretical results of Section 3, we expect to see the optimal convergence rates $e_{h, 0}^2 = \mathcal{O}(h^{k+1})$ and $e_{h, 1}^2 = \mathcal{O}(h^{k})$, assuming that the exact solution is smooth enough. Also, the error values in spatial directions for all examples are calculated at $T = 1$ and constant time steps with $K = 100$ that is $\tau = 1/100$ wherein time-step $\tau$ obeys the $CFL$ condition $\tau\approx \mathcal{O}(h^{k+1})$.

The purpose of introducing VEM is to take care of the two-dimensional spatial domain, where it is efficiently combined with a finite difference approach for the temporal direction, and we get a fully discrete scheme. In all the example cases, we have taken zero initial condition then only we are able to use the Grunwald-Letnikov approximation for Caputo type derivatives because with zero initial value, the Caputo type coincides with R-L type, and we can use the Grunwald approximation designed for R-L type derivatives. However, it is important to mention that in Section 5.2, we have taken a case of a non-smooth analytical solution in the temporal direction, which affects the convergence optimality of the temporal direction.

5.1.   Example 1

Consider the two-dimensional time-fractional nonlinear convection-diffusion equation (3.1), which models the anomalous processes, that is, the processes governed by random walks against the ideal Brownian motion ^[47] by applying the VEM scheme (4.3).

We let $\Omega = (0, 1)\times(0, 1)$, $J = [0, 1]$ with $\boldsymbol{b} = (1, 1)^{T}$ with the nonlinear term $f(u) = u+u^2$ and the RHS function $g(\boldsymbol{x}, t)$ given by

The virtual element solution and error plots at $\alpha = 0.8$, over the family of uniform square meshes and regular Voronoi meshes are shown in Figures 2 and 3 for $k = 1$, and likewise, for $k = 2$ in Figures 4 and 5. The $L^2$ norm and $H^1$ semi-norm error obtained along with the convergence rates for linear order VEM are presented in Tables 1 and 2 for the uniform square meshes and in Tables 3 and 4 for the regular Voronoi meshes, wherein Table 3 presents the error values in spatial direction with respect to different time step lengths and for different $\alpha$ values. We report the corresponding results for $k = 2$ in Tables 5–8. In addition, Table 9 demonstrates that the VEM (4.3) is of global order $\mathcal{O}(\tau)$ at time $T = 1$ and $T = 0.1$. The error and convergence in temporal direction are calculated by refining $\tau$ within a range of $\frac{1}{2^n}, n = 1, \ldots, 4$ for fixed mesh size $h = 10^{-3}$, wherein convergence rate in the time direction is calculated by $rate = log2\frac{(e(\tau_{n})}{(e(\tau_{n+1}))}.$

5.2.   Example 2

Consider the nonlinear term $f(u) = u+u^3$ with $\boldsymbol{b} = (1, 1)^{T}$ and define the function $g(\boldsymbol{x}, t)$ such that the exact solution is given by,

over domain $\Omega = (0, 1)\times(0, 1)$ and time interval $J = [0, 1]$. Figure 6 and 7 show the virtual element solution and error plots at $\alpha = 0.6$, over the uniform square and regular Voronoi meshes for $k = 1$. Figure 8 and 9 show the virtual element solution and error plots at $\alpha = 0.8$, over uniform square and regular Voronoi meshes for $k = 2$. We report the approximation errors measured using the $L^2$ norm and $H^1$ semi-norm and the convergence rates for $k = 1$ when applying the VEM for $k = 1$ on the uniform square meshes in Tables 11 and 12, and on the regular Voronoi meshes in Tables 13 and 14, wherein Table 13 presents the error values in spatial direction with respect to different time step lengths and for different $\alpha$ values. We report the results for $k = 2$ in Tables 15–18. This example is different in a sense, as in we have assumed the analytical solution in time direction as $t^\alpha$. The optimal order of convergence, that is $\mathcal{O}(\tau)$, is achieved. The error table 10 provides the corresponding values.

5.3.   Example 3

In this problem, we choose a variable velocity field $\boldsymbol{b} = (y, 0)^{T}$, the nonlinear term as $f(u) = u^3$ and set $g(\boldsymbol{x}, t)$ from the non-smooth exact solution $u = tx^{1.6}$ over unit square domain. We report the $L^2$ norm and $H^1$ semi-norm and convergence rates over uniform square meshes using VEM of order $k = 1$, in Tables 19 and 20 respectively, while over non-convex mesh configuration results are reported in Tables 21 and 22. Figure 10 shows the VEM solution and an error plot at $\alpha = 0.8$ for non-convex meshes. Due to the weak regularity of the exact solution, an optimal convergence rate is not achieved for higher-order VEM ($k\geq 2$).

Tables 3 and 13 for Examples 1 and 2 respectively show how the error's concur using VEM of order $k = 1$ over regular Voronoi mesh configurations. For a suitable range of $h$ (0.25-0.01) and $\tau$ (0.05-0.01) we see with increasing time steps/decreasing temporal mesh size the error values have very slight changes that infers, as we keep on increasing time steps with constant spatial mesh sizes the error values have very slight changes.

6.   Conclusions

Time-fractional models are of great significance in process modeling, necessitating the development of efficient numerical techniques for their approximation. This paper provides a comprehensive study of the Virtual Element Method in the context of time-fractional nonlinear convection-diffusion equations. We introduce the VEM, its formulation, and a thorough discretization scheme for the aforementioned equations. Theoretical error estimates are rigorously derived and subsequently validated through numerical examples. This work represents an extension of the VEM approach to address time-fractional PDEs, offering a detailed analysis of both theoretical and computational aspects. In the future, we plan to explore the application of VEM to various fractional models and equations while also comparing their efficiency against existing numerical methods. Our primary focus will be on approximating real-time application problems, including anomalous physical processes, impulse, and fluid transport in complex geometries, and the impact of turbulence on design considerations.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory partially supported the work of GM under project number 20220129ER. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). G. Manzini is a member of the Gruppo Nazionale Calcolo Scientifico-Istituto Nazionale di Alta Matematica. Also, Vellore Institute of Technology, Vellore supported the work of M. Chandru under a SEED grant (Sanction Order No. SG20230081).

Conflict of interest

The authors declare that they have no conflicts of interest.