The Efficient Finite Element Methods for Time-Fractional Oldroyd-B Fluid Model Involving Two Caputo Derivatives

: In this paper, we consider the numerical schemes for a time-fractional Oldroyd-B fluid model involving the Caputo derivative. We propose two efficient finite element methods by applying the convolution quadrature in time generated by the backward Euler and the second-order backward difference methods. Error estimates in terms of data regularity are established for both the semidiscrete and fully discrete schemes. Numerical examples for two-dimensional problems further confirm the robustness of the schemes with first- and second-order accurate in time.


Introduction
In this paper, we investigate efficient numerical schemes for the time-fractional Oldroyd-B fluid model involving the Caputo derivative in time. Let Ω ⊂ R n (n = 1, 2, 3) be a bounded convex domain with a smooth boundary ∂ , T > 0 be a fixed time. The equation is stated as follows.
with a homogeneous Dirichlet boundary condition u(x, t) = 0 on ∂ × (0, T] and the initial value conditions u(x, 0) = v(x) and ∂ t u(x, 0) = b(x). Here, the parameters μ, a 1 and a 3 are some fixed positive constants, the source term f and the initial value conditions v and b are given functions. The fractional orders γ ∈ (0, 1] and α ∈ (1,2]. The fractional derivative C D ν 0,t with ν ∈ (n − 1, n] is in Caputo sense defined by where the Riemann-Liouville (R-L) integral RL D −ν 0,t with ν > 0 is given by In order to obtain the desired convergence order, many efficient schemes often rely on some requisite regularity on the solution. Sometimes such assumption on the solution are too strong and make those numerical schemes less efficient or even useless in practical application.
In this paper, we aim to propose the robust and efficient numerical schemes for the fractional model (1) by applying the popular convolution quadrature method. The convolution quadrature method for solving fractional model has received considerable attenction [16,17]. The convolution quadrature method developed by Lubich [18] can provide a flexible framework for constructing some reliable numerical methods. So in this paper, we deduce two fully discrete schemes for (1) based on the linear finite element method in space and convolution quadrature in time generated by the backward Euler and the second-order backward difference methods. By the strategy developed in [19], we investigate the error estimates with respect to data regularity (Theorems 4.3 and 4.5).
We find that only smooth data of v fulfill the requirement for the second-order spatial accuracy of the proposed schemes, while the nonsmooth data (that is v ∈ L 2 ( )) would deteriorate the spatial accuracy. However, both smooth and nonsmooth data v can ensure the desired temporal convergence orders in the fully discrete schemes. We make some discussions on such situation (Remarks 4.1, 4.3, and 4.5) and numerically confirm our conclusions (Tabs. 7 and 8).
To the best of our knowledge, it seems that such results have not been found in the existing literatures yet.
The rest of this paper is organized as follows. In Section 2, we introduce some preliminaries which are need in the numerical study. Two fully discrete schemes using the Galerkin finite element method in space and convolution quadrature in time are developed in Section 3. In Section 4, the error analysis of the schemes are investigated. In Section 5, numerical examples for two-dimensional problems are provided to support our theoretical results. The conclusions are presented in Section 6.
Throughout this paper, the notation c is denoted as a constant which may vary at different occurrences, but is always independent of the mesh size h and the time step size τ .

Preliminary
We first list some useful notations in this part.
(2) Laplace transform and a related property. The Laplace transform of f with respect to t is denoted aŝ One can readily derive that (3) Riemann-Liouville derivative and its relationship with Caputo derivative. For n − 1 < ν < n, The definition of R-L derivative is given by The relationships between Caputo and R-L derivatives are as follows.
We also need the notations from the operational calculus [21].
Let K(z) be a complex valued or operator valued function that is analytic in a sector θ where θ ∈ π 2 , π . Besides, the function K(z) is bounded by for some real numbers λ and M. Then K(z) is the Laplace transform of a distribution k on the real line, which vanishes for t < 0, has its singular support empty or concentrated at t = 0, and which is an analytic function for t > 0 ( [18,21] for more details). By the inversion Laplace transform for K(z) with t > 0, we get Here the contour lies in θ , and parallel to its boundary and oriented with increasing imaginary part. We denote that θ ,δ = {z ∈ C : |z| = δ, | arg z| ≤ θ} ∪ z ∈ C : z = ρe ±iθ , ρ ≥ δ .
We define K(∂ t ) as the operator of convolution with the kernel k : K(∂ t )g = k * g. Here ∂ t is the time differentiation and g is suitably smooth. We divide the time interval [0, T] into a uniform grid with a grid point t k = kτ . Here the time step size τ = T/N with a positive integer N. The convolution quadrature K(∂ τ )g of K(∂ t )g at t = t n is given by where the quadrature weights ω k (τ ) are determined by the generating function Here ϑ is the quotient of the generating polynomials of linear multistep method [18]. In this paper, we focus on the backward Euler (BE) and the second-order backward difference (SBD) methods, for which ϑ(ζ ) = 1 − ζ and ϑ(ζ ) = (1 − ζ ) + (1 − ζ ) 2 /2, respectively.
The convolution quadrature has the following error estimates [22].
Lemma 2.1. Let K(z) be analytic in θ and (3) hold. Then for g(t) = ct σ −1 , the convolution quadrature based on BE (p = 1) or SBD (p = 2) satisfies Finally, we state some basic properties for the functions h (z) = 1/ z + a 1 z α and g (z) = z + a 1 z α / μ 1 + a 3 z γ with z ∈ θ by the following lemma.

Fully Discrete Schemes
In this section, we propose two fully discrete finite element methods by using the standard Galerkin finite element method in space and convolution quadrature in time.

Semidiscrete Galerkin Scheme in Space
A partition of the domain is denoted by T h in which h is the maximal length of the sides of the triangulation T h . Then the continuous piecewise linear finite element space V h is given by The L 2 ( ) orthogonal projection P h : and the Ritz projection R h : From the definitions of P h and R h , one can see that they are stable in L 2 and H 1 0 , respectively. We remark that these properties would be used frequently and may not be mentioned explicitly in the error estimates. We also need the following approximation properties of these two operators P h and R h , one may also refer to Theorem 1.1 and Lemma 1.1 in the book [20].
Lemma 3.1. The following approximation properties of the operators P h and R h hold: The semidiscrete Galerkin scheme for (1) reads: with the initial value conditions Here v h and b h are proper approximations to the functions v and b, respectively. The bilinear form a(u, v) is given by (∇u h , ∇χ ).
By introducing the discrete Laplacian Δ h : V h → V h with the definition: and letting A h = −Δ h , we can rewrite the semidiscrete scheme (8) as follows.
where f h (t) = P h f (t).

Fully Discrete Schemes with Convolution Quadrature in Time
In this part we apply convolution quadrature based on BE and SBD formulas in time to obtain two fully discrete schemes.
Using the relationship between Caputo and R-L derivatives, we derive from the semidiscrete scheme (9) that on both sides of the semidiscrete scheme (10), we obtain with t > 0. Combining the convolution quadrature with the associativity of convolution, we have the approximation of u h (t n ) at t n = nτ with U n h by Remark 3.1. Actually, the term ∂ α τ tb h in Eq. (13) may also be replaced with the exact expression . Next, we present a robust numerical scheme which can maintain the second-order accuracy for the scheme (12). Since the semidiscrete scheme (11) can be rewritten as where Letting 1 τ = (0, 3/2, 1, 1, . . .), and noting 1 τ = ∂ τ RL D −1 0,t , we obtain the SBD scheme as follows: Find U n h for n ≥ 1 such that that is,

Error Estimates
In this section, we derive the error estimates for the proposed numerical schemes by using the strategy developed in [19] which is originated in [21]. For the convenience of discussion, we always assume a 1 = 1 when necessary.

Error Estimates for the Semidiscrete Scheme
We first consider the homogeneous problem with f = 0 and derive the corresponding integral representation of the solution.
An application of the Laplace transform to (1) yields wherê and the functions h(z) and g(z) are given by Eq. (7).
So the solution u(t) can be represented by Similarly, the solution u h (t) to Eq. (9) can be represented by the following: wherê From the structures ofÊ (z) andÊ h (z), we let F h (z) = (g(z)I + A) −1 − (g(z)I + A h ) −1 P h for notation conveience. The operator F h (z) has the following property which plays a key role in the error estimate [9]. Lemma 4.1. The following estimate holds: Now we are ready to present the error estimates for the semidiscrete problem.  (1) with v ∈Ḣ 2 (Ω) and f = 0, and set u h be the solution of Eq. (9) with v h = R h v. Then, the following estimates hold: For the first term I, by the identity and noting that A h R h = P h A, we can obtain By Lemmas 3.1 and 4.1, we obtain Note that the last two terms II and III can be recast as Thus, for the second term II, by Lemma 4.1, we have Putting the bounds of I, II, and III together, we complete the proof for the case (a) with L 2 -estimates.
For the second case b ∈ L 2 ( ), by the error estimates in Eq. (20) for the first case (a), one just need to bound the third term III in Eq. (20) with P h replacing R h . By Lemma 4.1, we derive that So the L 2 -estimates are proved completed. A similar argument yields the H 1 -estimates for both cases (a) and (b).
By Lemma 3.1, we conclude that the above error estimate is not O(h 2 ) for the case v ∈ L 2 ( ) due to the appearance of the term v − v h . Now we consider the inhomogeneous problem with f = 0. By applying the Laplace transform, the solution u(t) can be represented by where By a similar argument, the solution u h (t) to Eq. (9) is represented by where Subtracting Eq. (21) from Eq. (22), we get The operatorG h has the following error estimate.
Proof. By Lemma 4.1, we readily obtain A similar argument also yields The proof is thus completed.
We are in position to state the error estimate for the inhomogeneous problem.
Proof. Using Lemma 4.2 and Hölder's inequality for the error e h (t) in Eq. (23), we can readily derive the desired result. One may refer to Theorem 3.4 in [10] for similar discussions. Thus the proof is completed.

Error Estimates for the Fully Discrete Schemes
In this part, we derive the L 2 -error estimates for the fully discrete schemes (13) and (16).

Error Estimates for the BE Method
We first consider the homogeneous problem with f = 0.
Since the spatial error u(t) − u h (t) has been discussed in the previous part, we focus on the temporal error u h (t n ) − U n h . We state the error estimates as below. Lemma 4.3. Let u h and U n h be the solutions of Eq. (9) and Eq. (13), respectively. Particularly, v ∈Ḣ 2 (Ω) , f = 0, and v h = R h v. Then the following error estimates hold: Proof. By Eqs. (11) and (12), we have and So it follows from Lemma 2.1 (with σ = 1, p = 1, and the suitable chosen parameters λ) that where the estimate for the term R h v in the last inequality has used the following two facts: H 1 0 (Ω)-stability of R h and the equivalence between ∇v and v 1 in H 1 0 (Ω). For b ∈ L 2 ( ) and b h = P h b, by the bound of F 2 (z) and Lemma 2.1 (with σ = 2, p = 1, and some suitable parameters λ), the second term II can be derived that This together the bound I yields the case (b). The case (a) can be derived similarly, so the proof is completed.

Remark 4.2. When proving the error estimate in terms of smooth data inḢ 2 (Ω), we only can obtain the terms with v h and b h instead of those with A h v h and A h b h , and then bound them by Av and Ab in view of v h = R h v and b h = R h b. One may wonder whether this result can be improved to A h v h and A h b h
where we can naturally have Av and Ab by using the identity A h R h = P h A and the L 2 -stability of P h . We find that such situation may not happen unless the term ∂ t u in Eq. (1) vanishes.

Remark 4.3. By the equalitŷ
one can improve the above result with respect to the initial data v ∈ L 2 ( ) and v h = P h v, that is Indeed, one may consider the term μa 3Êh (z) z γ A h = : F 12 (z) in Eq. (24). By noting that a use of Lemma 2.1 (with σ = 1, p = 1, and the suitable chosen parameters λ) yields In view of Lemma 4.3, we readily obtain the desired results.
We now summarize the results for the error estimates of the fully discrete BE scheme (13) in the next theorem.
Next we state the error estimates for the inhomogeneous problem with f = 0 but v = b = 0.
Proof. It suffices to bound U n h − u h (t n ) . From (11) and (12), we have By the Taylor expansion For the first term I, noting that F (z) ≤ c|h (z) |, we apply Lemma 2.2 and Lemma 2.1 (with σ = 1, p = 1, and some suitable parameters λ) to derive that The second term II has the following estimate: Combining the above estimates I and II with Theorem 4.2, we complete the proof.

Error Estimate for the SBD Method
We now present the error estimates for the homogeneous problem in the next lemma.  (14) and (15), one has From Lemma 2.2, we have and So for the case b ∈ L 2 and b h = P h b, by Lemma 2.1 (with σ = 2, p = 2, and the suitable chosen parameters λ), we derive that The error estimates for the case b ∈Ḣ 2 (Ω) and b h = R h v is similar with that in (b), so the proof is completed.
It seems that the convergent order may maintain the first/second-order temporal accuracy for the semidiscrete scheme (14) by using the BE/SBD scheme, even though the initial data is nonsmooth, i.e., v ∈ L 2 ( ), while one need to impose the condition a 3 = 0 in order to obtain second-order spatial accuracy in the spatial error estimate of the semidiscrete scheme, see Remark 4.1 for more details.
Now we summarize the error estimates of the SBD scheme (16) for the homogeneous problem.
Finally, we present the following error estimate for the inhomogeneous problem with v = b = 0. Theorem 4.6. Let u be the solution of the problem (1) with v = b = 0 and f ∈ L∞(0, T; L 2 ( )), and let U n h be the solution of (16) with (s) ds)).
Proof. Similar to the proof of Theorem 4.4, we just need to bound e n h = u h (t n ) − U n h . Using the Taylor expansion of f h (t) at t = 0: f h (t) = f h (0) + tf h (0) + t * f h , we rewrite the expressions of u h (t n ) and U n h in (14) and (15) respectively as Noting that for ∀z ∈ θ , one has So for the first two terms I and II, by Lemma 2.1 (with σ = 2, p = 2, and the suitable chosen parameters λ), we obtain For the third term III, we similarly have Putting the bounds together we thus finish the proof.

Numerical Experiments
In this part, we demonstrate numerical examples to verify the theoretical results derived in previous section. In all numerical tests, we fix coefficient μ = a 1 = a 3 = 1, unless otherwise specified. The spatial and temporal convergence rates are computed separately on the domain = (0, 1) 2 and at a fixed final time T. The L 2 -norm errors are measured by U n h − u h (t n ) L 2 (Ω) at t n = T. Example 5.1. (The solution is known). Consider problem (1) with the following source term: The corresponding solution to (1) is u(x, t) = sin(π x) sin(π y)(1 + t p ). Taking p = 2.5 and T = 1, we compute the numerical solution of (1) by using the BE (13) and SBD scheme (16), respectively. The convergence rate O(τ 2 + h 2 ) in the L 2 -norm can be observed from the numerical results in Tabs. 1 and 2, which are in line with the theoretical results.
We compute a reference solution in each case on very refined mesh due to the exact solutions of (1) with the above data are difficult to obtain. The numerical results for cases (a)-(d) with T = 0.1 are demonstrated on Tabs. 3-6, respectively. We observe the O(τ ) and O(τ 2 ) rates for the temporal error for the BE and SBD schemes of (1), respectively. These observations agree well with the convergence theory. We also numerically examine the errors for the nonsmooth initial data v ∈ L 2 ( ) by using the BE and SBD schemes for (1), in which the error estimates are not available in this paper, see Tabs. 7 and 8. We observe that the temporal ones are O(τ ) and O(τ 2 ) with three different combination of α and γ , which verify the discussions in Remarks 4.3 and 4.5. From Tab. 8, one can see that the spatial errors are affected by the parameter a 3 . The smaller value of a 3 , the smaller the spatial errors. This further numerically confirms the conclusion in Remark 4.1.        The computational domain is on a circle = (x, y)|x 2 + y 2 = 1, in which we are more interested. The numerical results are demonstrated in Fig. 1.
For fixed a 3 = 0.01 in Fig. 1, the peak of the solution is located at the origin at the beginning. However as time goes by, say, t = 0.1, the peak of the solution to the Eq. (1) is no longer at the origin, but spreads around the origin. This phenomenon becomes more apparent with increasing time and exhibits a period-like evolution. Compared to the case a 3 = 0.01, the evolution of the solution in case a 3 = 0.1 becomes slower, see the right-hand column in Fig. 1 for the details. From Fig. 1, we can also observe that the solution of the model is not significantly different for various parameters a 3 = 0.01 and 0.1 at the beginning t = 0.01. As time goes by, the effect of the value of parameter a 3 on the model becomes more and more obvious. It seems that the smaller the parameter a 3 , the more the fractional model (1) behaves like the time-fractional diffusionwave equation with damping [13]. This may further validates the claims in [9] that the fractional derivative acting on the Laplacian in the model (1) may be used to capture the viscoelastic behavior of the flow.

Conclusions
In this paper we have proposed two efficient and robust fully discrete schemes for the time-fractional Oldroyd-B fluid equation which involves the Caputo derivative in time. These two numerical schemes are based on Galerkin finite element method in space and convolution quadrature in time. The error estimates for the schemes are investigated with respect to data regularity. Numerical results are provided and they are all well agree with the theoretical results.
We find that the spatial errors of the proposed schemes would be deteriorated if the data v ∈ L 2 ( ) due to the appearance of the term C D γ 0,t acting on the Laplacian in (1), while the temporal errors are not affected. It seems that such situation does not exist in other similar fractional model which involves the R-L derivative in time [9,10]. So this makes the numerical investigation of fractional model considered in this paper more difficult. Nevertheless, we make some comments on such situation and numerically confirm our conclusions, see Remarks 4.1, 4.3, and 4.5, and Tabs. 7 and 8.
Our strategy in this paper can be extended to solve other types of fractional models, for example, the non-Newtonian fluid model [2] ∂ t + a 1C D α 0,t u (x, t) which can be regarded as a more general case of (1). It seems that the numerical study is trivial so we omit the details here. It is noteworthy that our numerical schemes are not very computationally efficient for long time problems due to the nonlocal nature of Caputo derivative and the high dimensionality of the space. Recently, some authors have proposed fast algorithms to resolve such issue, [23][24][25], just to name a few. So combining the numerical schemes here with a fast algorithm is an interesting thing to do and worthy of further study.