A Finite Difference Method for Solving Unsteady Fractional Oldroyd-B Viscoelastic Flow Based on Caputo Derivative

,


Introduction
Recently, non-Newtonian fluid models including fractional derivatives have aroused great interest because of their extensive application in practical engineering problems. Dassios and Baleanu [1] studied singular linear systems of finite difference equations by using the Caputo fractional derivative and its two alternative forms. Alsharif and Elmaboud [2] investigated the physical characteristics of electroosmotic flow by solving the fractional Cattaneo model. Khan et al. [3,4] obtained some analytical solutions for the fractional Oldroyd-B fluid with different initial and boundary conditions. Based on the Fourier sine transform, Nadeem [5] presented exact solutions for the fractional Oldroyd-B fluid in periodic motion. Qi and Xu [6] obtained the analytical solution for the Oldroyd-B viscoelastic flow by using the discrete Laplace transform. However, it is very difficult to carry them out in the field of engineering applications. On the one hand, analytical solutions always appear in complex series form. On the other hand, the Mittag-Leffler function or the Fox-H function is often included in the exact solution. To overcome the above problems, researchers studied the Oldroyd-B viscoelastic flow with other methods. In combination with the numerical inversion of Laplace transforms, Huang et al. [7] analyzed the motion characteristics of the generalized second fluid in a double barrel rheometer. Yao and Zhang [8] set up a kind of implicit difference scheme for the viscoelastic fluid following a generalized Maxwell fractional differential equation. Wang et al. [9] studied the unsteady Poiseuille flow of the fractional Oldroyd-B viscoelastic fluid between two parallel plates and used the Stehfest algorithm to carry out numerical inversion of the Laplace transform. Liu and Zheng [10] used the Laplace transform coupling with the Hankel transform to study unsteady helical flows and obtained exact solutions of a generalized Oldroyd-B fluid between two infinite concentric cylinders. Yu [11] adopted a high-order compact finite difference method to solve the generalized fractional Oldroyd-B fluid model and mainly investigated the corresponding stability and convergence for the finite difference method. However, it should be pointed out that the related work on the numerical analysis is much less compared with that on viscoelastic rheological behaviors and analytical solutions. With the help of a discretization of the Caputo time-fractional derivative and the finite difference method, Zhang et al. [12] obtained an approximate implicit difference method to solve a two-dimensional multiterm timefractional Oldroyd-B equation. Al-Maskari and Karaa [13] considered the numerical approximation of a generalized fractional Oldroyd-B fluid problem with a semidiscrete scheme based on the piecewise linear Galerkin finite element method. In other engineering fields, Raslan et al. [14] present a computational scheme using the bi-finite difference method (Bi-FDM) to solve the hyperbolic telegraph equation in two dimensions. The fluid flow model governed by partial differential equations is reduced to a system of nonlinear ordinary differential equations for solving the calculation in reference [15]. However, no relevant literature has been seen in which the method is applied to fractionalorder models. It shows that the finite difference method has great advantage in solving the fractional Oldroyd-B fluid problem.
Large quantities of experiments and engineering practices have shown that the fractional derivative model can describe the constitutive relation of viscoelastic materials and their rheological properties more accurately than other models. Moreover, traditional constitutive models can be used as special cases of fractional derivative constitutive models. Therefore, the fractional derivative is widely used in the expanding research of the constitutive model. However, the research of applying the fractional-order derivative to the Oldroyd-B model is still relatively rare, especially when the analytical and numerical solutions are given in specific arithmetic cases and compared. The fractionalderivative-type viscoelasticity principal structure relationship established by the theory on fractional derivatives was a newly developed principal structure model nearly two decades ago. Compared to the Riemann-Liouville [16] fractional derivative, Caputo fractional derivatives [14] are more suitable and convenient for numerical analysis. Khan and Rasheed [17] have adopted the Galerkin finite element method to determine the approximate solution that is uniform with the finite difference approximation of the Caputo fractional time derivative. In another article by them, the Caputo fractional derivative is used to simulate the physical properties of the fractional unsteady magnetohydrodynamic flow [18]. Zaid and Dumitru [19] studied and proposed a new system of fractional differential equations in conjunction with the generalized Caputo fractional derivative.
Erturk et al. [20] developed a mathematical analysis for a corneal-shaped model in combination with the Caputo fractional derivative. These show that the Caputo fractional derivative is widely used in mathematical and physical problems Moreover, the initial and boundary conditions defined by the Caputo fractional derivatives are specific in physical meaning, which can provide great convenience for computation and application in the field of physics and engineering. Therefore, this paper adopts Caputo fractional derivatives directly.
In summary, it has been known that the fractional derivative is widely used in the expanding research of the constitutive model to solve mathematical and physical problems, but the research of applying the fractional-order derivative to the Oldroyd-B model is still relatively rare, especially when the analytical and numerical solutions are given in specific arithmetic cases and compared. This research idea is executed in this paper, and four sections are mainly included here. The finite method is given in Section 2, including different schemes, the proof of the stability, and convergence, as well as the analysis of unique solvability. Taking the unsteady Poiseuille flow and oscillatory flow in a circular pipe as example, numerical simulations are developed in Section 3. Then, the final conclusion is given in Section 4. Moreover, the proof of the numerical simulation method and analytical solution is given in the appendix.

Finite Difference Method
2.1. Governing Equations with Caputo Derivative. The fundamental equations describing the unsteady motion of an incompressible fluid are given [21]: in which V = ðu, v, wÞ is the velocity vector and u, v, and w that are, respectively, along the x, y, and z directions in the Cartesian coordinate system. ρ is the uniform density, and d/dt denotes the material derivative with respect to time. T = −pI + σ is the Cauchy stress tensor, in which p represents the pressure, I is the unit tensor, and σ is the extra stress tensor. The constitutive equation of the fractional Oldroyd-B fluid [3] is where μ is the viscosity and λ 1 and λ 2 are the relaxation time and retardation time, respectively. α and β are the order of the fractional derivatives with the condition 0 ≤ α ≤ β ≤ 1. A 1 is the first Rivlin-Ericksen tensor given by 2 Advances in Mathematical Physics where δ α σ/δt α and δ β A 1 /δt β are defined as Here, the fractional differential operator δ α /δt α based on Caputo's definition [8] is expressed as Here, Γð⋅Þ is the Gamma function. Besides, the Laplace transform of the time-fractional Caputo derivative is given as and the time-fractional Caputo derivative is defined as The definition is stricter than that of the Riemann-Liouville fractional-order derivative, in which f ðtÞ is an integrable function of order n.
For the unidirectional flow, the velocity and stress are considered as follows: where i is the unit vector along the x direction. Equations (5) and (9) satisfy the initial condition σðy, 0Þ = 0; thus, σ xz = σ yz = σ yy = σ zz = 0 is established. Corresponding to Equation (1), the following expressions are obtained: In the absence of body forces, the momentum Equation (2) for the unidirectional flow of the fractional Oldroyd-B fluid is written as According to Equations (12) and (14), the motion equation of the fractional Oldroyd-B fluid is expressed as where A = ð1/ρÞð∂p/∂xÞ, ν = μ/ρ.

Difference Scheme.
In the space region ½0, d, let us set y i = ih, i = 0, 1, 2, ⋯, M, t n = nτ, n = 0, 1, ⋯, h = d/M is the space step, and τ is the time step. Here, u n i represents the approximate value of uðih, nτÞ, and it is assumed that u n = ðu n 1 , ⋯, u n M Þ is a M-dimensional vector. For convenience, the following marks are introduced: The Crank-Nicolson (C-N) scheme is used to discretize partial derivative terms in Equation (11). At node ðy i , t n−1/2 Þ, it is satisfied as t n−1/2 = ðn − 1/2Þτ. The time partial and space derivative with two orders are, respectively, expressed as In combination with Equation (9) and the above C-N scheme, the implicit difference scheme of Equation (13) at node ðy i , t n−1/2 Þ is expressed as The detailed derivation process is given in Appendix A.
In addition, the initial and boundary conditions are also discretized according to the above steps.

Advances in Mathematical Physics
The function u is defined in Theorem 2. The implicit difference scheme of Equation (18) is unconditionally stable.
Proof. Supposing u n i is the exact solution of Equation (18). Multiplying Equation (18) by hτδ t u n i , and summing up for i from 1 to M − 1; then, for n from 1 to N, it can be obtained as The following emphasis is focused on analyzing Equation (19). In combination with Equations (14) and (15), it can be concluded that The detailed proof of Equation (20) is given in Appendix B. According to Lemma 1 [22], we can obtain Therefore, the implicit difference scheme is unconditionally stable.

Theorem 3. Equation (18) is convergent when it is satisfied by h ≤ τ.
Proof. Suppose that uðy i , t n Þ is the exact solution of Equation (13) at node ðy i , t n Þ; u n i is the numerical solution, and the error between two solutions is defined as e n i = uðy i , t n Þ − u n i . Here, assume that there exists no error for the initial data, according to Equations (17) and (18); then, Here, 1 ≤ i ≤ M − 1, n ≥ 1. The initial and boundary values of the error are It follows from Equation (22) that We have Because ke n k ∞ ≤ ð ffiffiffi d p /2Þkδ y e n k 2 is established, it can be further obtained that where C 1 is a constant which has nothing to do with h and τ. Therefore, according to h ≤ τ, when it is satisfied that h ⟶ 0 and τ ⟶ 0, then ke n k ∞ ⟶ 0 is established, and the convergence is proven.

Advances in Mathematical Physics
Proof. Equation (16) can be written as Here, it is satisfied that Equation (27) contains M − 1 equations, and the part of the coefficient can be written as a matrix form: u is the velocity vector and b is the right-hand side Because the coefficient matrix is strictly diagonally dominant, the difference scheme Equation (18) is uniquely solvable.

Numerical Study of Poiseuille Flow with
Finite Difference Method Figure 1 shows the two-dimensional Poiseuille flow between two parallel plates. In order to verify the finite difference method, the analytical solution of the Poiseuille flow is firstly given: where The above analytical solution satisfies the following initial and boundary conditions: Because the form of the analytical solution is very complex, we only consider the approximate analytical solution under the condition that t is a small number. In order to execute the numerical simulation, it is supposed that ν = d = 1 and A = −1. In Figure 2, the velocity on the centerline of the two parallel plates is considered, and numerical comparison between the finite difference method and the approximate analytical solution is shown. Here, three groups of parameters are chosen: In the first two group parameters, the variation induced by α and β is compared. In the first and third   7 Advances in Mathematical Physics above four parameters. The velocity based on the parameter of group 1 is larger than other two group parameters. Meanwhile, Table 1 displays the relative error of the velocity on the centerline with time for three group coefficients. For example, Error 1 corresponds to group 1. It demonstrates that the relative error gradually turns to be within the scope of 2% as time increases, especially for group 2 and group 3. It powerfully states that the finite difference method is stable and reliable. The relative error of the velocity on the centerline with time for the above three group coefficients is given in Table 1, and the small relative error further reflects the accuracy of the numerical algorithm.
Next, the influence of four parameters on the velocity and fluid physical properties is considered with the finite difference method. Figure 3 displays the influence of the parameter α on the velocity at the center line with the other   Advances in Mathematical Physics three fixed parameters; here, it is supposed that β = 1:0, λ 1 = 0:1, and λ 2 = 0:01. When the parameter α = 0 or α = 1, the velocity represents the moving velocity of the secondorder fluid or integer-order Oldroyd-B fluid. It can be observed that the velocity distribution with α = 0 and α = 1 is in good agreement with that based on integral-order constitutive Equation (18). When it satisfies 0 < α < 1, the velocity of the fractional Oldroyd-B fluid is between that of the secondorder fluid and that of the integer-order Oldroyd-B fluid. It also shows in Figure 3(a) that the velocity gradually tends to be stable with time, and maximum values of the velocity increase with α. This phenomenon is called velocity overshoot, which appears as the increasing elasticity and a fluctuation in velocity. As displayed in Figure 3(b), the distribution of flow is filed symmetric along the y direction with different α. Meanwhile, the closer the parameter α is to 1, the more obvious is the velocity overshoot. The physical characteristics of the fluid is manifested as an increase in fluid elasticity and appearances of velocity fluctuations, and the flow field variation is also closer to that of the integer-order Oldroyd-B fluid. In   Figure 4 mainly studies the influence of the order parameter β on the velocity at the center line with fixed parameters; here, it is satisfied as λ 1 = 0:1 and λ 2 = 0:01. As shown in Figure 4(a), the velocity increases gradually at the steady state, but the variation of the amplitude reduces with β. It illustrates that the fluid viscosity reduces with β; this caused the amplitude of the velocity to reduce gradually. As shown in Figure 4(b), when α is large, the variation of the velocity is basically the same. The parameter β almost has no influence on the amplitude of velocity and the time to steady state.
There appears a velocity overshoot phenomenon which is identical with Figure 3. It can be concluded that variations of parameters α and β can lead to velocity overshoot; this is because an increase in α causes an increase in the elasticity of the fluid, while an increase in β causes a decrease in the viscosity of the fluid.

Numerical Simulations of the Oscillatory Flow in a Circular Pipe
The oscillatory flow in a circular pipe exists extensively in the blood, petroleum, and polymer solutions. This kind of flow not only has the viscosity of the fluid but also shows the elasticity of the solid. Therefore, it is difficult to use

10
Advances in Mathematical Physics one model to represent all the characteristics of viscoelastic fluids, and many viscoelastic models and constitutive equations have emerged. Fractional Oldroyd-B has made a breakthrough in the research of material constitutive theory. In general, only a simple analysis of the velocity and stress is performed in the literature, but the variation of periodic velocity and stress amplitude is rarely studied. In this paper, the velocity and shear stress of the oscillatory flow in a circular pipe is studied in detail. The laminar flow in a circular pipe with an infinite radius is considered, and it is supposed that (1) the flow is incompressible, (2) the velocity along the radius direction is equal to zero, (3) the flow is axisymmetric about the pipe, and (4) the axial velocity is only related with the radius.
In order to verify the correctness of the finite difference method, the analytical solution of the oscillatory flow is firstly derived in detail, and then, numerical simulations of the velocity and shear stress are executed and compared between two kind of methods. The cylindrical coordinate system with ðr, θ, zÞ is established here, in which the z axis is along the central axis of the circular pipe. Here, the processing of the difference format for shear stress is similar to that for the velocity.

Advances in Mathematical Physics
Starting from Equations (10) and (11), the variable separation method is used to derive the analytical solution.
Here, the dimensionless analytical amplitudes of the velocity and shear stress on the pipe wall are given: The detailed process is contained in Appendix C.
According to the expressions of A u * and A σ * , the amplitudes of velocity and stress are related with the pipe radius R * , oscillation frequency ω * , and specific value θ = λ 2 /λ 1 , as well as the order of the fractional Oldroyd-B constitutive model. In the numerical simulation, the influence of these parameters is studied.
4.1. The Influence of Parameters on the Velocity. Next, numerical simulations of the oscillatory flow in a circular pipe are executed. Firstly, the finite difference method is verified by comparing numerical results and the analytical solution. When the order parameters are chosen as α = β = 1, the fractional Oldroyd-B constitutive model degenerates to the classical Oldroyd-B constitutive model. Figure 5 shows the variation of the velocity amplitude with the oscillation frequency. It demonstrates that numerical results are in good agreement with the analytical solution for different pipe radii. The velocity amplitude basically shows a decreasing trend with the pipe radius. In Figure 5(a), the resonance phenomenon appears at ω * = 48, because the peak value of the velocity amplitude is much larger than at other frequencies.
In Figures 5(b)-5(d), the number of the formant is gradually increasing with the pipe radius. Meanwhile, the peak value substantially reduces especially for the velocity amplitude at R * = 5 . Compared with the case of R * = 0:05, 0:1, 0:5, 5, all formants concentrate at 0 ≤ ω * ≤ 10. The curve of the velocity amplitude gradually approximates to a monotonic decreasing curve with R * , which is very close to the curve of the Newtonian fluid.
Secondly, a different order parameter is chosen to study the variation of the velocity amplitude. Figure 6 shows the velocity amplitude with α = 0:8 and β = 0:95, and Figure 7 chooses the order parameter as α = 0:2 and β = 0:35. It can be observed that the formant also appears in the velocity amplitude curve of the fractional Oldroyd-B fluid similar to the classical Oldroyd-B fluid. On the one hand, the first peak appears in Figure 7(a) corresponding to R * = 0:05, and the number of the formant increases, but the peak value substantially reduces with R * . On the other hand, the location of the peak value moves left along the horizontal axis. Otherwise, when the parameter is chosen as α = 0:8 and β = 0:95 in Figure 7, the velocity amplitude also decreases quickly comparing with Figures 6(a) and 6(b). In brief, the velocity amplitude reduces gradually with R * and the parameters α and β. Comparison with Figure 5 shows that the physical variation laws of fractional-order Oldroyd-B fluids are similar to those of classical Oldroyd-B fluids. It further shows that the oscillation peak and oscillation frequency can be reduced by reasonably controlling the above parameters.

Advances in Mathematical Physics
In addition, we discuss the influence of α and β on the velocity amplitude. In Figure 8, the influence of α with fixed β and the influence of β with fixed α are all shown. It can be found that there are multiple formants in the curve with α = 0:9, but only one formant exists with α = 0:6 and α = 0:3 in Figure 8(a). However, there only exists one formant with different β in Figure 8(b). The same characteristic is that the first peak with different α and β corresponds to the same frequency. Meanwhile, the peaks are gradually decreasing with α and β. On the whole, the velocity amplitude increases with the increase of α and decreases with β at low frequencies.

The Influence of Parameters on the Shear Stress.
Zhu et al. [23] studied physical characteristics of the Maxwell flow oscillating in a pipe. Their research shows that the resonance phenomenon occurred on the amplitude of the wall shear stress at certain frequencies, and this discovery is helpful to improve the core oil displacement efficiency in the practical engineering field. Therefore, the variation of the shear stress with different parameters and frequency is further studied in this part.
At first, the parameters are chosen as α = 1:0 and β = 1:0; the influence of the frequency on the shear stress is given in Figure 9. It can be observed that numerical results are in good agreement with analytical solutions. Similar to the velocity image, the resonance phenomenon also exists in the amplitude of the wall shear stress, especially for R * = 0:5 as displayed in Figure 9(b). For different R * , the variation tendency of A σ * is identical with A u * , such as R * = 0:05 and R * = 0:5. The first peak value reaches the maximum value, and other peak value shows a decreasing trend on the whole. Meanwhile, the number of peaks gradually increases with R * . The first peak value decreases and moves to the left along the horizontal axis. It can be inferred that when R * is greater than a certain value, the cure will approximate to a monotone decreasing cure. Finally, after R is greater than a certain value, 14 Advances in Mathematical Physics the stress change curve will approach a monotonically decreasing curve, close to that of a Newtonian fluid. Secondly, the influence of parameters on the wall shear stress is considered. For Figures 10 and 11, α = 0:8, β = 0:95, and α = 0:2, β = 0:35, are chosen, respectively. Figure 10 shows the same variation trend as Figure 9, and the resonance peak increases with R * . The difference is that the amplitude of A σ * is smaller compared with that in Figure 9. It can be concluded that the amplitude of A σ * varies with α, β, and R * . Here, R * has the greatest impact compared to the other two parameters. The larger R * is, the faster A σ * decreases with little α and β. The curve corresponding to A σ * is basically monotonically decreasing with frequency for different R * .
As Figure 11 shows, when both α and β are small, the curves obtained for different R * are monotonically decreasing. And the larger the R * , the faster the magnitude decreases. This physical phenomenon fully illustrates that parameters α and β influence the appearance of the resonance phenomenon. The classical viscoelastic model can be considered as a special case of the fractional-order viscoelastic model, and this is theoretically consistent with Zhu's findings.
Next, the parameters are chosen as R * = 0:5 and θ = 0:001, and the influence of different order parameters on

Advances in Mathematical Physics
the wall shear stress is studied in detail. As shown in Figure 12(a), the shear stress with β = 0:95 results in a resonance phenomenon at ω * ≈ 5:6 Hz, and the peak value increases dramatically with α at this frequency. Multiple peak values appear especially for α = 0:9, and the amplitude decreases with α when the frequency satisfies ω * > 10 Hz. The similar phenomenon also exists in Figure 12

16
Advances in Mathematical Physics only one peak value at ω * ≈ 2 Hz, and this indicates that the parameter β almost has no influence on the resonance phenomenon. However, as β decreases, the amplitude of the wall shear stress decreases very slightly at low frequencies, but there is a significant increase in the amplitude at high frequencies.

Conclusions
In the present paper, the difference scheme based on the differential equation with fractional derivatives is obtained. Then, the stability, convergence, and unique solvability are proven by the energy method. Finally, the finite difference method is verified with the Poiseuille flow and oscillatory flow in a circular pipe. The following conclusions are obtained: (1) The fractional mixed partial derivative in the differential equation is solved reasonably by integration, successfully avoiding the problem of low-order accuracy caused by discretizing the mixed derivative directly (2) The consistency of the numerical and exact solutions indicates that the method is reliable, and it is effective and concise to deal with applicable engineering problems. Moreover, according to the results in this paper, the finite difference method can be applied to differential equations with fractional mixed partial derivatives similar to Equation (13) (3) The finite difference method accurately captures the velocity overshoot phenomenon in the Poiseuille flow and resonance phenomena that appeared in the oscillatory flow. Based on these studies, it can be inferred that the same phenomena that occur in the classical viscoelastic fluid flow also occur in the fractional viscoelastic fluid flow. Otherwise, the above phenomena have a significant dependence on the order of fractional derivatives. Meanwhile, the classical viscoelastic model can be regarded as a special case of the fractional viscoelastic model. For example, the simulation of the oscillatory flow can well describe the physical deformation of materials, which is very helpful for the study of the properties of several materials in terms of damping

Appendix A The Proof of Finite Difference Scheme
For Equation (13), the time partial derivative with the α + 1 order is given: ðA:1Þ Then, the following formula is obtained: ðA:3Þ Based on the following formulas, ðA:4Þ
The expression (A.10) is further arranged as ðA:11Þ

B The Proof of Equation (20)
For the second item on the LHS of Equation (A.11), it is satisfied: ðB:1Þ And the RHS of the above inequality is further simplified as ðB:2Þ

Advances in Mathematical Physics
For the formula (B.2), the following inequality is established: Therefore, formula (B.1) can be equally expressed as ðB:4Þ Meanwhile, the equality is found for the RHS of Equation (B.4): ðB:5Þ Then, the following result can be easily obtained: Because A is a negative, we can obtain Þδ y δ t u k i 2 " Þδ y δ t u n i ð Þ À Á 2 " # ≤ 0:

ðB:8Þ
Taking into account Equations (B.7) and (B.8), it can be concluded that ðB:9Þ Namely, where σ rz is the rz component of shear stress and u z is the velocity component along z direction. The influence of the volume force is ignored, and the momentum equations with the nonslip boundary condition are obtained