Stokes’ First Problem for Viscoelastic Fluids with a Fractional Maxwell Model

: Stokes’ ﬁrst problem for a class of viscoelastic ﬂuids with the generalized fractional Maxwell constitutive model is considered. The constitutive equation is obtained from the classical Maxwell stress–strain relation by substituting the ﬁrst-order derivatives of stress and strain by derivatives of non-integer orders in the interval ( 0, 1 ] . Explicit integral representation of the solution is derived and some of its characteristics are discussed: non-negativity and monotonicity, asymptotic behavior, analyticity, ﬁnite/inﬁnite propagation speed, and absence of wave front. To illustrate analytical ﬁndings, numerical results for different values of the parameters are presented.


Introduction
Viscoelastic fluid flows with the classical Maxwell constitutive model [1] have been the object of intense study for many years, for a short survey see for example [2,3]. In the case of a unidirectional flow, the Maxwell constitutive relation has the dimensionless form where σ and ε are shear stress and strain, respectively, a is the relaxation time, and b the dynamic viscosity. Fractional calculus, in view of its ability to model hereditary phenomena with long memory, has proved to be a valuable tool to handle viscoelastic aspects [4]. For instance, rheological constitutive equations with fractional derivatives play an important role in the description of properties of polymer solutions and melts [5]. A four-parameter generalization of the classical Maxwell model (1) has been proposed in [6], see also [7,8] and [9] (Chapter 7). It is obtained from (1) by substituting the first order time-derivatives of stress and strain with Riemann-Liouville fractional derivatives in time, which leads to the following generalized fractional Maxwell constitutive equation where 0 < α, β ≤ 1 and a, b > 0. Since for α > β, the corresponding relaxation function is increasing, which is not physically acceptable, the constitutive model (2) is considered only for fractional order parameters satisfying the constraints (see [6][7][8][9]) The fractional Maxwell constitutive Equation (2) with thermodynamic constraints (3) has been shown to be an excellent model for capturing the linear viscoelastic behavior of soft materials exhibiting one or more broad regions of power-law-like relaxation [8,10,11]. Stokes' first problem is one of the basic problems for simple flows. It is concerned with the shear flow of a fluid occupying a semi-infinite region bounded by a plate which undergoes a step increase of velocity from rest. One of the first studies on Stokes' first problem for viscoelastic fluids is [12]. Since then many works have been devoted to this problem for different types of fluids with linear constitutive equations. Stokes' first problem for viscoelastic fluids with general stress-strain relation is studied in [13], [14] (Section 5.4) and [4] (Chapter 4). For some recent results concerning unidirectional flows of viscoelastic fluids with fractional derivative models we refer to [15] (Chapter 7). Let us note that the fractional order governing equations are usually derived from classical models by substituting the integer order derivatives by fractional derivatives. The problem of correct fractionalization of the governing equations is discussed in [16].
The classical Maxwell constitutive Equation (1) and its fractional order generalization (2) are among the simplest viscoelastic models for fluid-like behavior and the corresponding Stokes' first problem is a good model of essential processes involved in wave propagation. For studies on Stokes' first problem for classical Maxwell fluids we refer to [2,3] and the references cited therein. Stokes' first problem for generalized fractional Maxwell fluids is studied in [17,18], where explicit solutions in the form of series expansions are obtained (in [18] only the particular case β = 1 is considered). Explicit solutions of other types of problems for unidirectional flows of fractional Maxwell fluids with constitutive Equation (2) can be found in [19][20][21].
Despite the abundance of works devoted to exact solutions of unidirectional flows of linear viscoelastic fluids, to the best of our knowledge, there is very little analytical work enlightening the basic properties of solution to Stokes' first problem for fractional Maxwell fluids, such as propagation speed, non-negativity, monotonicity, regularity, and asymptotic behavior. In addition, it was pointed out in [22,23] that in a number of recent works the obtained exact solutions to Stokes' first problem for non-Newtonian fluids contain mistakes. In the present work we provide a new explicit solution to Stokes' first problem for generalized Maxwell fluids with a constitutive model (2) and compare it numerically to the one given in [17]. Moreover, based on the theory of Bernstein functions, we study analytically the properties of the solution and support the analytical findings by numerical results. This paper is organized as follows. In Section 2, the material functions of the fractional Maxwell model are studied. In Section 3, the solution to Stokes' first problem is studied analytically based on its representation in Laplace domain. Section 4 is devoted to derivation of explicit integral representation of solution and numerical computation. Some facts concerning Bernstein functions and related classes of functions are summarized in an Appendix.

Fractional Maxwell Model-Material Functions
Consider a unidirectional viscoelastic flow and suppose it is quiescent for all times prior to some starting time that we assume as t = 0. Since we work only with causal functions ( f (t) = 0 for t < 0), if there is no danger of confusion for the sake of brevity, we still denote by f (t) the function H(t) f (t), where H(t) is the Heaviside unit step function. We are concerned with the fractional Maxwell model, given by the linear constitutive Equation (2) with parameters satisfying (3). The material functions-relaxation function G(t) and the creep compliance J(t)-in a one-dimensional linear viscoelastic model are defined by the identities [4,24] σ( where the over-dot denotes the first derivative in time. Applying (A4), we derive from the constitutive Equation (2) expressions for material functions in the Laplace domain: Taking the inverse Laplace transform of the above representations of G(s) and J(s), we get, by using the properties (A5) and (A8), the following representations of the material functions (see also [6][7][8][9]) , t > 0, (6) in terms of the Mittag-Leffler function (A6) and power-law type functions, respectively. The second law of thermodynamics, stating that the total entropy can only increase over time for an isolated system, implies the following restrictions on the material functions: G(t) should be non-increasing and J(t) non-decreasing for t > 0. This behavior is related to the physical phenomena of stress relaxation and strain creep [4]. It has been proven that the fractional Maxwell model (2) is consistent with the second law of thermodynamics if and only if 0 < α ≤ β ≤ 1 [6][7][8][9]. For completeness, here we formulate and prove a slightly stronger statement. In fact, it appears that the monotonicity of the material functions of the fractional Maxwell model imply that the relaxation function G(t) is completely monotone and the creep compliance J(t) is a complete Bernstein function. The definitions and properties of completely monotone functions and related classes of functions are summarized in the Appendix.
is decreasing near 0, and thus (c) is violated. Therefore, any of conditions (b) and (c) implies (a). It remains to prove that (a) implies (d) and (e), i.e., G(t) ∈ CMF and J(t) ∈ CBF . Indeed, if 0 < α ≤ β ≤ 1 then t α−β ∈ CMF and E α,α−β+1 − 1 a t α ∈ CMF as a composition of the completely monotone Mittag-Leffler function of negative argument and the Bernstein function t α . Therefore, G(t) ∈ CMF is a product of two completely monotone functions. The fact that J(t) is a complete Bernstein function follows from the property of the power function Representations (6) together with the asymptotic expansion (A7) gives G(+∞) = 0 and J(+∞) = +∞. Therefore, the fractional Maxwell constitutive equation indeed models fluid-like behavior (for a discussion on the general conditions see e.g., [24], Section 2). More precisely, (6) and This means that only for β = 1 the relaxation function G(t) is integrable at infinity and the integral over (0, ∞) is finite, which is a stronger condition required for fluid behavior. In this case (β = 1) we obtain from (6) and (A9) the area under the relaxation curve The asymptotic expansions of the material functions for Therefore, if α < β then G(0 + ) = +∞ and J(0 + ) = 0, while if α = β different behavior is observed: This means that in the limiting case α = β the material possesses instantaneous elasticity and, therefore, finite wave speed of a disturbance is expected [4] (Chapter 4). This property will be discussed further in the next section.

Stokes' First Problem
Consider incompressible viscoelastic fluid which fills a half-space x > 0 and is quiescent for all times prior to t = 0. Rheological properties of the fluid are described by the fractional Maxwell constitutive Equation (2) with parameters satisfying the thermodynamic constraints (3). Taking into (2) can be rewritten in the form: The fluid is set into motion by a sudden movement of the bounding plane x = 0 tangentially to itself with constant speed 1. Denote by u(x, t) the induced velocity field. Assuming non-slip boundary conditions, the equation for the rate of strainε = ∂u/∂x, Cauchy's first law ∂u/∂t = ∂σ/∂x and the constitutive Equation (7), imply the following IBVP for the velocity field where H(t) is the Heaviside unit step function. Problem (8)-(10) is the Stokes first problem for a fractional Maxwell fluid, c.f. [17]. Let us note that by applying the Caputo fractional derivative operator C D 1−β t to both sides of Equation (8) and taking into account the zero initial conditions (10), and operator identities (A3), Equation (8) can be rewritten as the following two-term time-fractional diffusion-wave equation with Caputo fractional derivatives of orders 2 + α − β, 2 − β ∈ [1,2]. In a number of works various problems for the two-term time-fractional diffusion-wave equation are studied, e.g., [25], [26] (Chapter 6) and [27][28][29][30][31][32]. In this work, we study the Stokes' first problem following the technique developed in [33] for the multi-term time-fractional diffusion-wave equation.
Problem (8)-(10) is conveniently treated using Laplace transform with respect to the temporal variable. Denote by u(x, s) the Laplace transform of u(x, t) with respect to t. By applying Laplace transform to Equation (8) and taking into account the boundary conditions (9), initial conditions (10) and identities (A4), the following ODE for u(x, s) is obtained where Solving problem (11), with s considered as a parameter, it follows, c.f. [17] (Equation (22)) where g(s) is defined in (12).
Before taking inverse Laplace transform to obtain explicit representation of the solution, we first deduce information about its behavior directly from its Laplace transform (13).
First we prove that g(s) is a complete Bernstein function. Here this is deduced from the complete monotonicity of the relaxation function G(t) proved in Proposition 1. Therefore, the proposed method of proof can be used also for more general models, provided the corresponding relaxation function is completely monotone.
Proof. For the proof we use properties (D), (E) and (G) in the Appendix. Since G(t) ∈ CMF then G(s) ∈ SF by definition. Therefore, property (D) implies 1/ G(s) ∈ CBF . Since also s ∈ CBF (by (G)), g(s) = s/ G(s) is a product of two complete Bernstein functions, and (E) implies g(s) ∈ CBF . Theorem 1. The solution of Stokes' first problem u(x, t) is monotonically non-increasing function in x and monotonically non-decreasing function in t, such that 0 ≤ u(x, t) ≤ 1 and u(x, 0 + ) = 0, u(x, +∞) = 1. More precisely, for any x > 0, Proof. From Proposition 2 and properties (A), (B) and (C) in the Appendix we deduce (in the same way as in [33] (Theorem 2.2)): Therefore, (13) and Bernstein's theorem imply and, in this way, the non-negativity and monotonicity properties of the solution are proven. Taking into account that g(s) → 0 as s → 0 and g(s) → ∞ as s → ∞, it follows from (13) and, therefore, 0 ≤ u(x, t) ≤ 1. Since exp −x g(s) ∼ 1 − x g(s) and g(s) ∼ s −β /b for s → 0, then (13) implies which, by applying (A5) and Tauberian theorem, gives the asymptotic expansion (14) of u(x, t).
Next, propagation speed of a disturbance is discussed. From general theory, see e.g., [4] (Chapter 4) and [14] (Chapter 5), the velocity of propagation of the head of the disturbance is c = κ −1 , where From here, we easily obtain the propagation speed for the considered fractional Maxwell model Therefore, if 0 < α < β ≤ 1, the disturbance propagates with infinite speed. On the other hand, if the orders of the two derivatives in the model are equal 0 < α = β ≤ 1 (and thus the model possesses instantaneous elasticity as discussed earlier), the propagation speed c is finite. This means that u(x, t) ≡ 0 for x > ct. Let us consider the case of finite propagation speed in more detail. It is known that in the limiting case of classical Maxwell model (α = β = 1) the velocity field u(x, t) has a jump discontinuity at the planar surface x = ct, see e.g., [2]. What happens when the orders α and β are equal, but non-integer? According to [14] the amplitude of such a jump is exp(−ωx), where ω = lim s→∞ g(s) − κs . In our specific case (12) with α = β gives g(s) = as 2 b + s 2−β b = κ 2 s 2 1 + 1 as β and using the expansion Therefore, ω < ∞ only in the case β = 1, which means that only in the case of classical Maxwell fluid the wave u(x, t) exhibits a wave front, i.e., a discontinuity at x = ct. For equal non-integer values of the parameters α and β, 0 < α = β < 1, (18) implies ω = ∞, and thus there is no wave front at x = ct. Therefore, in this case, an interesting phenomenon is observed: coexistence of finite propagation speed and absence of wave front. This is a memory effect, which is not present for linear integer order models.
At the end of this section we consider in more detail the case 0 < α < β ≤ 1, for which we obtained c = ∞. We will prove next that in this case the solution u(x, t) is an analytic function in t.
The analyticity of the solution established in the above theorem implies that for any x > 0 the set of zeros of u(x, t) on t > 0 can be only discrete. This, together with the monotonicity of u(x, .) proved in Theorem 1, implies that u(x, t) > 0 for all x, t > 0. This observation confirms again that in the case 0 < α < β ≤ 1 a disturbance spreads infinitely fast.

Explicit Representation of Solution and Numerical Results
To find explicit representation of the solution we apply the Bromwich integral inversion formula to (13) and obtain the following integral in the complex plane where the Bromwich path has been transformed to the contour D ∪ D 0 , where In the same way as in [33] (Theorem 2.5) we deduce from (20) the following real integral representation of the solution. where Let us check that the integral in (21) is convergent. Indeed, since K − (r) ∼ r 1−β/2 as r → 0 then sin(rt − xK − (r))/r ∼ r −β/2 as r → 0. Therefore, the function under the integral sign in (21) has an integrable singularity at r → 0. At r → ∞ integrability is ensured by the term exp(−xK + (r)), since K + (r) > 0 and K + (r) ∼ r 1−(β−α)/2 → +∞ as r → +∞.
Next, the explicit integral representation (21) is used for numerical computation and visualization of the solution to the Stokes' first problem for different values of the parameters. For the numerical computation of the improper integral in (21) the MATLAB subroutine "integral" is used. The aim of our numerical computations is two-fold: to support the presented analytical findings and to compare our results to those in [17] (see their Figures 1 and 2), where plots of the solution are given, obtained from a different representation: a series expansion in terms of Fox H-functions. In order to compare our numerical solutions to the ones given in [17], we perform computations for the same parameters. In contrast to [17], where only plots of solution as a function of x are given, we plot it also as a function of t.
In Figure 1 the solution u(x, t) is plotted for β = 0.8 and four different values of α: 0.2, 0.4, 0.6, and 0.8; a = 10 α , b = 10 β−1 . The case of instantaneous elasticity α = β = 0.8 is added to those plotted in Figure 1 of [17], in order to show the different behavior. In this case, the disturbance at x = 0 propagates with finite speed c = √ b/a = 10 −1/2 and u(x, t) ≡ 0 for x > ct. There is no wave front at x = ct. In the other three cases with α < β the propagation speed is infinite and u(x, t) should be positive for all x, t > 0. However, far from the wall, it becomes negligible. As expected, the solution u(x, t) exhibits monotonic behavior: it is non-increasing with respect to x, see Figure 1a, and non-decreasing with respect to t, see Figure 1b. The plots presented in Figure 1a for α = 0.2, 0.4, 0.6 are identical with those given in Figure 1 of [17].
In Figure 2, the solution is plotted for α = 0.5 and three different values of β: 0.5, 0.7, 0.9, a = 10 α , b = 10 β−1 . The case of instantaneous elasticity now is α = β = 0.5 with finite propagation speed c having the same value as above. Again, as expected, there is no wave front at x = ct. Monotonicity of u(x, t) is clearly seen. The plots presented in Figure 2a for fixed time t = 4 are identical with those given in [17], Figure 2. We refer also to [33] (Figures 1-3) for more illustrations of the three different types of behavior, related to Maxwell fluids, namely: finite wave speed and presence of wave front (α = β = 1), finite wave speed and absence of wave front (0 < α = β < 1), and infinite propagation speed (0 < α < β ≤ 1). Let us note that the paper [18] also presents plots of the solution of Stokes' first problem for fractional Maxwell fluids with β = 1. However, the plots in Figures 5a, 7a and 9a of [18] do not seem to be monotonically decreasing with respect to the spatial variable, which is a contradiction with the expected monotonic behavior of the solution.

Conclusions
The solution of Stokes' first problem for a viscoelastic fluid with the generalized Maxwell model with fractional derivatives of stress and strain is studied based on its Laplace transform with respect to the temporal variable. The thermodynamic constraints on the fractional parameters imply physically reasonable monotonic behavior of solution with finite propagation speed (without wave front) when the two fractional orders coincide and infinite propagation speed otherwise. An explicit integral representation of the solution is derived and used for its numerical computation. The obtained numerical results are in agreement with those reported in [17].
The explicit integral representation of the solution can be further used to derive explicit or asymptotic expressions for other characteristics of the flow, such as shear stress at the plate, thickness of the boundary layer, etc.
The presented technique can be applied in the study of different generalizations of the considered model, e.g., the viscoelastic model based on Bessel functions in [34] and the fractional order weighted distributed parameter Maxwell model in [35].
The characterization of the class CMF is given by the Bernstein's theorem: a function is completely monotone if and only if it can be represented as the Laplace transform of a non-negative measure (non-negative function or generalized function).
The class of Stieltjes functions (S F ) consists of all functions defined on (0, ∞) which can be written as Laplace transform of a completely monotone function. Therefore, SF ⊂ CMF .
A non-negative function ϕ on (0, ∞) is said to be a Bernstein function (ϕ ∈ BF ) if ϕ (x) ∈ CMF ; ϕ(x) is said to be a complete Bernstein functions (CBF ) if and only if ϕ(x)/x ∈ SF . We have the inclusion CBF ⊂ BF .
A selection of properties of the classes of functions defined above is given next.
For proofs and more details on these special classes of functions we refer to [38].