Exact, approximate and numerical solutions for a variant of Stokes׳ first problem for a new class of non-linear fluids

Author's Accepted Manuscript

Exact, Approximate and Numerical Solutions for a Variant of Stokes’ First Problem for a New class of Non-Linear Fluids K.V. Mohankumar, K. Kannan, K.R. Rajagopal

www.elsevier.com/locate/nlm

PII: DOI: Reference:

S0020-7462(15)00130-4 http://dx.doi.org/10.1016/j.ijnonlinmec.2015.07.004 NLM2517

To appear in:

International Journal of Non-Linear Mechanics

Received date: 7 July 2015 Revised date: 13 July 2015 Accepted date: 14 July 2015 Cite this article as: K.V. Mohankumar, K. Kannan, K.R. Rajagopal, Exact, Approximate and Numerical Solutions for a Variant of Stokes’ First Problem for a New class of Non-Linear Fluids, International Journal of Non-Linear Mechanics, http://dx.doi.org/10.1016/j.ijnonlinmec.2015.07.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Exact, Approximate and Numerical Solutions for a Variant of Stokes’ First Problem for a New class of Non-Linear Fluids K.V. Mohankumara , K. Kannana,∗, K.R. Rajagopalb a Department

of Mechanical Engineering, Indian Institute of Technology Madras, Chennai 600036 of Mechanical Engineering, Texas A&M University, College Station, Texas 77840

b Department

Abstract Stress power-law fluids are a special sub-class of fluids defined through implicit constitutive relations, wherein the symmetric part of the velocity gradient depends on a power-law of the stress (see Eq. (2.2)), was introduced recently to describe the non-Newtonian response of fluid bodies. Such fluids are counterparts to the classical power-law fluids wherein the stress is given in terms of a power-law for the symmetric part of the velocity gradient. Stress power-law fluids can describe phenomena that cannot be described by the classical power-law fluids (see Percalova and Prusa (2015) [1]). In this paper, first a new exact solution for a variant of Stokes’ first problem for stress power-law fluids, wherein the exponent n = 0 (Navier-Stokes fluid), is obtained. Such an exact solution for the stress is in terms of a convolution integral, for which we establish bounds. We then compute the convolution integral using Gauss-Kronrod quadrature by ensuring that its value always lies within the bounds. Using the validated quadrature, we can accurately evaluate the exact solution and we use it to validate the numerical scheme employed in solving the governing equations for the stress-power law fluids with arbitrary exponent n. Finally, for stress power-law fluids wherein the exponent n < 0 (stress-thickening fluids), we obtain an approximate solution for the stress that agrees well with the numerical solution. Keywords: Implicit constitutive theories, Stress power-law fluids, Variation of Stokes’ first problem, exact solution, approximate solution, parabolic cylinder function, modified theta function, hypergeometric function, convergence, bounds 1. Introduction Many real fluids exhibit response characteristics such as shear-thinning/shear thickening (which is reflected by the generalized viscosity depending on the Frobenius norm of the symmetric part of the velocity gradient), stress relaxation, non-linear creep, normal stress differences in simple shear flows, dependence of the material moduli on the invariants of the stress, yielding1 , thixotropy etc., that cannot be described by the classical Navier-Stokes fluid model. Complicated fluid models have been developed to describe the non-Newtonian behavior exhibited by many real fluids, and amongst them one of the most popular ∗ Corresponding

author Email address: [email protected] (K. Kannan ) 1 Yielding is supposedly the phenomenon wherein the fluid starts to flow once a certain threshold in stress is overcome. However, if one defines a fluid as a material that cannot resist shear, then a material that exhibits the phenomenon of yielding cannot be a fluid. In order to clearly define the phenomenon one needs to delineate how long the body can resist the applied stress. While it might seem to resist it for a very short duration it might flow on prolonged application of the stress. Thus, to clearly define what one means by a fluid, one needs to take into account both time and length scales (scale at which motion is measured). The concept of what is meant by a fluid is very far from being clear and definitely imprecise (see the extended discussion in Malek and Rajagopal (2005)[2]). Preprint submitted to International Journal of Non-Linear Mechanics

class of models is the class of classical power-law fluids2 (Huilgol (1975)[3], Schowalter (1978) [4]). The classical power law fluid (see equation (2.1) that follows) is one in which the generalized viscosity is a function of the symmetric part of the velocity gradient, the function being a power law function. This class of fluids can describe the shear thinning/thickening that is displayed by many fluids; it cannot however describe the other non-Newtonian characteristics such as stress relaxation, dependence of the material moduli on the invariants of the stress such as the mean normal stress (which in an incompressible fluid is the mechanical pressure), normal stress differences in simple shear flow, etc. The classical power law model cannot describe interesting data concerning the flows of colloids and other solutions wherein the generalized viscosity cannot be expressed as a function of the symmetric part of the velocity gradient (see Percalova and Prusa (2015) [1]). Recently, Rajagopal (2003) [5], (2006) [6], has introduced a very general class of implicit models3 for describing the non-linear response of fluids which include the clas2 The class of power-law fluids is a sub-class of the Generalized Stokesian fluids wherein the Cauchy stress is expressed as a nonlinear function of the symmetric part of the velocity gradient. 3 Prusa and Rajagopal (2012) [7] have considered a general class of implicit models wherein the history of the stress and the relative deformation gradient are related. The class of models considered by

July 14, 2015

sical fluid models such as the classical power-law fluids, but it also includes a large class of other new non-Newtonian fluid models. A sub-class of the same is models wherein one has an explicit expression for the symmetric part of the velocity gradient as a non-linear function of the stress, and a further sub-class is the class of stress power-law models wherein the symmetric part of the velocity gradient has the mathematical structure given by equation (2.2) that follows. This class of fluids exhibits very interesting response characteristics that cannot be captured by the classical power-law fluids described by equation (2.1), and the distinction in the characteristics of the flow of the two classes of power-law fluids undergoing simple shear flow is studied by Malek et al. (2010) [8]. Of particular interest with regard to the flow characteristics of a stress power-law fluid is the possibility of the existence of a limiting shear rate in simple shear flows for a range of values of the power-law index. We note that for the values of the power-law index n ∈ (−∞, − 12 ) in equation (2.2), the shear rate never exceeds a maximum value, however large the stress, that is the fluids are strain rate limited, a very interesting characteristics. Such a response cannot be captured by the classical stress power-law fluid model. As discussed by Percalova and Prusa (2015) [1], colloids and biological solutions containing DNA strands exhibit response wherein the strain rate is not a monotonic function of the stress and modification of the powerlaw model considered here, and that by Le Roux and Rajagopal (2013) [9] are good candidates to describe the flows of such systems. Malek et al. (2010) [8] studied the counterparts of the classical plane Couette and Poiseuille flow, HagenPoiseuille flow, and cylindrical Couette flow of stress powerlaw fluids. Even within the class of the above simple flows, Malek et al. (2010) [8] find palpable departures in the fluid response from that exhibited by the classical power-law fluids. The study by Malek et al. (2010) [8] has been followed by other analyses of the steady and unsteady flows of stress power-law fluids by Narayanan and Rajagopal (2013) [10] which essentially extended the seminal studies of unsteady flows due to an infinite accelerating plate and an infinite oscillating plate by Stokes (1851) [11] (Rayleigh (1911) [12] also studied the flow due to an infinite oscillating plate), and the present study is a natural extension of the same. Here we study a variant of the Stokes’ first problem for the stress power-law fluid, wherein the fluid is at constant stress until time t = 0 and is subjected to a constant time rate of shear stress at t > 0 at one of the boundaries while the other is undisturbed. The subsequent unsteady motion of the fluid is studied for various values of the power-law index and other material parameters. Unlike the classical power-law fluid wherein one substitutes the constitutive relation for the stress into the balance of linear momentum and obtains a partial differential equation for the ve-

locity field and the indeterminate pressure, which is solved jointly with the constraint of incompressibility for the velocity and the mechanical pressure, the stress power-law fluids satisfy the constraint of incompressibility implicitly. Therefore, one can combine the stress power-law form of constitutive equation and the balance of linear momentum equation into a governing equation in terms of the stress. We find that the nature of the solution for the unsteady stress and the velocity changes dramatically with the power-law index. It is necessary to generate a body of literature for the stress power-law fluid so that we can evaluate the efficacy and usefulness of the model in describing natural phenomena and this paper is one such addition to the literature. In the next section we introduce stress power-law fluids and in Section 3, we derive the equations governing the problem under consideration. In order to derive the governing equations, we adopt a semi-inverse approach wherein we assume special forms for both the stress and the velocity fields. The problem reduces to solving two equations, one for the balance of linear momentum and the other for the constitutive relation, simultaneously. The equations can however be reduced to a single equation of higher order for the shear stress. In Section 4, we consider the special power-law exponent n = 0, and in this case the higher order equation for the stress reduces to that for the Navier-Stokes fluid. This equation can be solved exactly using Laplace Transforms and the solution is given in terms of a term involving the parabolic cylinder function and a convolution integral involving the modified theta function and the complimentary error function (see equation (4.15)). The integrand has a singularity but to show that the integral is finite, we establish upper and lower bounds for the solution and use these bounds to verify the validity of the numerical evaluation of the solution (4.15). We can do this as the bounds are very tight that they lie practically on one another for small values of time and thus act as a good check of the numerical calculations for small values of time. Under certain conditions, the higher order stress equation can be approximated by an equation with a simpler structure and we obtain a solution to the approximate equation for negative values of the exponent n in Section 5. Finally, we obtain numerical solutions for the variant of the Stokes’ problem that is considered in the paper. 2. Stress Power-Law Fluids Before introducing the stress power-law fluid, we shall document the classical power-law fluid so that we can clearly see the differences between the two models. In the classical power-law fluid, the Cauchy stress T is related to the symmetric part of the velocity gradient through 1

T = −pI + µ[1 + γ(tr D 2 ) 2 ]m D.

(2.1)

where the term −pI denotes the indeterminate part of the stress due to the constraint of incompressibility, µ, γ, m

Prusa and Rajagopal includes most of the popular models used in non-Newtonian fluid mechanics as special sub-classes.

2

where ρ denotes the density. The term 31 tr(T ) is the mechanical pressure in the fluid and we shall denote it by p, that is

are constants, and D is the symmetric part of the velocity gradient. The stress power-law fluid is defined through the constitutive relation 1

D = α[1 + β(tr T d 2 ) 2 ]n T d .

(2.2)

p :=

where α, β, n are constants and T d is the deviatoric part of the stress defined through 1 T d = T − (tr T )I. 3

(2.3)

ρ

√ ∂u = 2α[1 + 2β|S|]n S. ∂y

y¯ =

We are interested in studying unidirectional unsteady flows wherein the stress and velocity take the form:

(3.2)

y u ¯ S t , u ¯= , S¯ = , , t= h Uo (h/Uo ) ρUo 2 √ β¯ = 2βρUo 2 , and Re = ρUo hα

(3.10)

where Re is the Reynold’s number. On substituting (3.10) into the equations (3.7) and (3.8), we obtain the dimensionless forms of the balance of linear momentum and constitutive relation to be: ∂ S¯ ∂u ¯ = , ¯ ∂t ∂ y¯

(3.11)

∂u ¯ ¯ S|] ¯ n S¯ , = 2Re [1 + β| ∂ y¯

(3.12)

ρ

On substituting equations (3.1) and (3.2) into the balance of linear momentum, and on ignoring the body force, we obtain ∂( 1 tr(T )) ∂u ∂S ρ = + 3 . ∂t ∂y ∂x

(3.8)

When S¯ = 0, the absolute value function |S| is not differentiable and therefore study of cases where stress changes sign is an issue. In eliminating the velocity field from equations (3.7) and (3.8), we have increased the order of the equation for the stress and thus need some additional initial and boundary conditions for the stress. We now proceed to render the above equations non-dimensional by introducing the following dimensionless quantities:

3. Governing equations

v = u(y, t)ex .

(3.7)

On differentiating the balance of linear momentum equation (3.7) partially with respect to y and substituting the constitutive relation (3.8) in it, we obtain the differential equation in terms of stress as √ 1 ∂2S ∂((1 + 2β|S|)n S) = for S¯ 6= 0 . (3.9) ∂t 2ρα ∂y 2

where κ is the shear rate and τ is the shear stress. Usually one expresses the shear stress as the product of the generalized viscosity (which is a function of the shear stress) and the shear rate. Thus, we can view γ as the inverse of the generalized viscosity only when the original expression is invertible. Also as the generalized viscosity increases, γ decreases and vice-versa. Thus, we can think of γ as a measure of the fluidity of the body. We can then associate stress-thickening when dγ dτ < 0 as the fluidity is decreasing, while if dγ > 0 the fluidity increases and hence the fluid dτ stress-thins. In the case of a stress power-law fluid, the fluid stress-thickens when n < 0, and stress-thins when n > 0.

(3.1)

∂u ∂S = . ∂t ∂y

On substituting (3.1) and (3.2) into the constitutive relation (2.2), we obtain

(2.4)

T d = S(y, t)(ex ⊗ey + ey ⊗ex ).

(3.6)

Assuming that the mechanical pressure in the fluid (3.6) is a constant, the equations of balance of linear momentum (3.3) through (3.5) reduce to

Notice that the model defined through equations (2.2) and (2.3) automatically satisfies the constraint that the fluid can only undergo isochoric motions; in other words the fluid is incompressible. We do not need to introduce a Lagrange multiplier associated with the constraint within this framework. In the simple unidirectional flow that is being considered, one can express the constitutive equation in the form: κ = γ(τ )τ,

1 (tr(T )). 3

(3.3)

and the equation (3.9) takes the form (for S¯ = 6 0) 0=

∂( 13 tr(T )) . ∂y

(3.4)

¯ S|) ¯ n S) ¯ ∂((1 + β| 1 ∂ 2 S¯ = . ∂ t¯ 2Re ∂ y¯2

0=

∂( 13 tr(T )) . ∂z

(3.5)

Without loss of generality, we take β¯ = 1 in (3.12) and 1 (3.13) (which implies β = √2ρU in (3.9) and (3.10)), and 2 o

3

(3.13)

the dimensionless form of the constitutive relation reduces to

approximate the series, as a sum of its first two terms, i.e. if (3.18) and the condition

∂u ¯ ¯ n S¯ , = 2Re [1 + |S|] ∂ y¯

¯ − |S|| ¯ 1 >> |n|S|

(3.14)

are satisfied, one can show the terms in the series (3.17) are such that     ¯ − |S|| ¯ >> n − 1 |S| ¯ k >> n − 1 |S| ¯ k+1 1 >> |n|S| k k+1 (3.20)

and the dimensionless form of (3.13) (for S¯ = 6 0) reduces to ¯ n S) ¯ 1 ∂ 2 S¯ ∂((1 + |S|) . = ∂ t¯ 2Re ∂ y¯2

(3.15)

where n and Re are non-dimensional parameters. One can either solve the system (3.11) and (3.14) for the unknown dimensionless velocity and stress simultaneously or one can solve (3.15) for the stress and then solve for the velocity field.

where k ≥ 2 and an approximation of the series (3.17) is given through ¯ n−1 ≈ 1 + (n|S| ¯ − |S|) ¯ . (1 + |S|)

¯ + |n||S| ¯ ≥ |n|S| ¯ − |S|| ¯ , |S|

¯ + |n||S| ¯ , 1  |S|

(3.23)

it immediately implies ¯ , 1  |S|

(3.24)

¯ . 1  |n||S|

(3.25)

and

Further, on substituting (3.22) in (3.23), we obtain ¯ + |n||S| ¯ ≥ |n|S| ¯ − |S|| ¯ . 1  |S|

(3.26)

Therefore, for the special case where the dimensionless stress S¯ and n satisfy (3.23), which implies that (3.18) and (3.19) are satisfied through (3.26), one can substitute (3.21) in (3.16) and obtain an approximation of the PDE to be   2¯ 1 1 ∂ S ∂ S¯ ≈ . ¯ ¯ ¯ ¯ ¯ ∂t 2Re (1 + n|S| − |S|) (1 + (|S| + n|S|)) ∂ y¯2 (3.27)

¯ n−1 in the equaThe series expansion of the term, (1 + |S|) tion (3.16) is given through  ∞  X n−1 ¯k |S| (3.17) k

k=2



n−1 and n is arbi:= (n−1)(n−2)(n−3)...(n−k) k! k trary. In general, the largest term (absolute value) in the expanded series (3.17) depends on n, but one can show that the series (3.17) has absolute convergence if

where

¯ 1 > |S|

(3.22)

¯ + |n||S| ¯ to be of much Now, if we let the upperbound |S| smaller value than 1, i.e.

1 ∂ 2 S¯ ∂ S¯
2 . (3.16) = ¯ n−1 (1 + |S| ¯ + n|S|) ¯ ∂ y¯ ∂ t¯ 2Re (1 + |S|)



(3.21)

However, using the triangle inequality, one can show that an upperbound to the absolute term in (3.19) is given through

3.1. An Approximation of the governing PDE In this section we first establish sufficient conditions under which the partial differential equation (3.15) can be approximated by the partial differential equation (3.28) that follows. These conditions will not be met in general. These conditions pertain to the magnitude of the stress that is engendered due to applied boundary conditions and thus by appropriately controlling the boundary condition it should be possible to ensure that these conditions are met. Whether this is indeed so can only be determined a’posteriori and as we shall see, we can ensure the fulfilment of the condition and hence the validity of the approximation. If one considers that the dimensionless stress in (3.15) satisfy sufficient smoothness condition with respect to time and space without changing sign, the non-linear differential equation (3.15), for S¯ 6= 0, reduces to the form

¯ n−1 = 1 + n|S| ¯ − |S| ¯ + (1 + |S|)

(3.19)

Using (3.24) and (3.25) in (3.27), we neglect the higher ¯ 2 , n|S| ¯ 2 and n2 |S| ¯ 2 and further reduce the order terms |S| above PDE to the form   2¯ ∂ S¯ 1 1 ∂ S ≈ (3.28) ¯ ∂ y¯2 . ∂ t¯ 2Re 1 + 2n|S|

(3.18)

is satisfied. However, in addition to satisfying the above criterion, if we assume that, the first term of the series (3.17) has the highest absolute value such that it is much larger than the absolute value of the second term, it implies that the absolute value of k + 1th term is smaller than the kth term where k ≥ 1 and therefore, one can simply

where Re and n are arbitrary. We will show that a solution to the above approximate PDE, satisfying the criterion (3.23), is an approximate solution to the original differential equation given in (3.15) or its reduced form (3.16). 4

4. Variant of Stokes’ first problem

Laplace transform method. Applying the Laplace trans¯ y , t¯)) = s¯(¯ form L(S(¯ y , r) to equations (4.4), (4.1) through (4.3), reduces them to

We shall suppose that the fluid is considered to be bounded between two parallel planes and assume that the flow is unidirectional parallel to the planes. We solve the dimensionless governing differential equations (3.11) and (3.14) for the dimensionless stress and velocity, in a finite domain, satisfying specific initial and boundary conditions. Let us consider a variant of Stokes’ first problem, for the new stress power-law fluid. We call it as ‘a variant of Stokes’ first problem’ as the flow is confined to the finite region between the two parallel plates and appropriate initial and boundary conditions are assumed on the shear stress unlike the original Stokes’ problems where the domain is semi-infinite and the initial and boundary conditions are imposed on the kinematic variable, namely velocity. In the original Stokes’ first problem, at initial time, the fluid is assumed to be at rest everywhere and a sudden increase of velocity to a constant value is imposed on one of its boundaries while the other semi-infinite boundary is kept at rest. However, in this variant of Stokes’ first problem that we define here, we assume that at initial time (t = 0), the fluid is at a constant stress (i.e stress is uniform and therefore no stress gradient) and at time t > 0, a constant time rate of shear stress is imposed on one of its boundary while the other finite boundary is always maintained at its initial stress value. The appropriate dimensionless form of the initial and boundary conditions, are given through ¯ y , 0) = S¯o , for y S(¯ ¯ ∈ (0, 1),

d2 s¯ − (2Re)r¯ s + (2Re)¯ s=0, d¯ y2

(4.5)

s¯(¯ y , 0) = S¯o ,

(4.6)

s¯(0, r) =

¯ t¯) = S¯o for ¯t > 0. S(1,

(4.7)

S¯o r

(4.8)

and s¯(1, r) =

respectively. We find that the exact general solution to the equation (4.5), satisfying the initial condition (4.6), is of the form √

s¯(¯ y , r) = C3 e−¯y

2Rer

√

+ C2 ey¯

2Rer

+

S¯o , r

(4.9)

where the constants C2 , C3 are determined from the boundary conditions (4.7) and (4.8). The above solution is exact in a finite domain y¯ ∈ (0, 1) and the constants C2 , C3 have to satisfy the conditions C3 + C2 =

(4.1)

λo , r2

(4.10)

and √

¯ t¯) = λo t¯ + S¯o for ¯t > 0 and S(0,

S¯o λo + 2 , r r

C 3 e−

(4.2)

2Rer

√

+ C2 e

2Rer

=0,

(4.11)

where Re ≥ 0 and r > 0. The constants C2 , C3 obtained from (4.10) and (4.11) are given through

(4.3)

√

λo e− 2Rer √ C2 = − r2 2 sinh( 2Rer)

We carry out a parametric study of the non-linear PDE (3.15) for different values of the non-dimensional parameters. We exhibit an explicit exact solution when the exponent n = 0. Further, for n < 0 we find an approximate solution. We also obtain numerical solutions for n < 0, n = 0 and n > 0.

and √

λo λo e− 2Rer √ C3 = 2 + r r2 2 sinh( 2Rer)

4.1. Exact unsteady solution for n = 0 We recall that when n = 0, the governing equation (3.15) reduces to the special case of an incompressible Navier-Stokes fluid as 1 ∂ 2 S¯ ∂ S¯ = . ∂ t¯ 2Re ∂ y¯2

(4.12)

(4.13)

On applying (4.12) and (4.13) in (4.9), the Laplace transform reduces to ! √ e−¯y 2Rer s¯(¯ y , r) = λo r2 ! ! √ √ e− 2Rer sinh(¯ y 2Rer) S¯o √ . − λo + 3/2 r r r1/2 sinh( 2Rer) (4.14)

(4.4)

To our knowledge, an exact solution to the Navier-Stokes equation for the boundary conditions (4.1) through (4.3) and the initial condition for the stress has not been documented and hence we proceed to do so. Let us determine the exact solution to the above equation using the

¯ y , t¯) Applying inverse Laplace transform L−1 (¯ y , r) = S(¯ and using the convolution theorem to obtain inverse Laplace 5

transform of the second term on RHS in (4.14) (which is a product of two transformed functions), we obtain r !! 3/2 −y ¯2 Re Re 2 ¯ ¯ S(¯ y , t¯) = λo √ t¯ e 4t D−3 y¯ t¯ π Z t¯   y¯ τ  −1 √ − λo | θˆ4 2 2Re 2Re 0 s !! r √ t¯ − τ 2(−Re Re ¯ 2 dτ e t−τ ) − 2Re Erfc π 2(¯t − τ ) + S¯o

Table 1: Some specific values of θˆ4 (ˆ z |tˆ):

(4.15)

0

0

−∞

tˆ > 0

0

finite

finite

1 2

zˆ =

1 2

√

τ On substituting z = y¯ √Re in (4.16), zˆ = y2¯ and tˆ = 2Re in t¯ (4.17) and using them in (4.15), one can obtain the final ¯ y , t¯). solution S(¯ In the integral term of the solution (4.15), the integrand has singularity at τ = 0 (for y¯ = 1, we have
1 τ ˆ lim θ4 2 | 2Re = −∞) and the exact integral is not known τ →0

explicitly for y¯ 6= 1. However, to show that the integral is finite, we determine the bounds for the convolution integral using bounds for the integrand. We proceed, next, to obtain the bounds for the modified theta function, the convolution integral and the dimensionless stress. These bounds will be used later to accurately evaluate the solution (4.15). From (4.17), one can show that, for a given zˆ and tˆ, the modified theta function reduces to evaluation of a continuous function of n ˆ,   2n−1 ˆ 2n−1 ˆ 2 ˆ 2 ˆ 1 e−( 2 +ˆz) /t − e−( 2 −ˆz) /t , (4.20) f (ˆ n) = √ π tˆ

e . On applying the recurrence and special relations for ν = −3, the parabolic cylinder function reduces to r   z 2 + 1 π z2 z z −z2 4 D−3 (z) = e Erfc √ − e 4 . (4.16) 2 2 2 2 The modified theta function, θˆ4 (ˆ z |tˆ) is defined by ∞  X 2n−1 ˆ ˆ 2 ˆ 1  −( 2n−1 +ˆ z )2 /tˆ 2 √ − e−( 2 −ˆz) /t . e π tˆ n ˆ =1 (4.17)

Since equation (4.15) is the exact solution to the linear PDE (4.4) subject to the boundary and initial conditions (4.1) through (4.3), it follows that

at discrete values of n ˆ (i.e. at n ˆ = 1, 2, 3, ..∞) and summing them up. For n ˆ ≥ 1, 0 ≤ zˆ ≤ 21 and tˆ ≥ 0, one can show that  2n−1 ˆ 2 ˆ 1  −( 2nˆ −1 +ˆz)2 /tˆ 2 0 ≥ f (ˆ n) = √ e − e−( 2 −ˆz) /t (4.21) π tˆ

Z t¯ 

=

lim

0 < zˆ <

From (4.17), one can show that the modified theta function is a monotonically decreasing function of zˆ (i.e. ∂∂zˆ (θˆ4 (ˆ z |tˆ) ≤ 0) for 0 ≤ zˆ ≤ 12 with its maximum value being zero, which implies that the modified theta function is non-positive valued such that 1 1 0 ≥ θˆ4 (ˆ z |tˆ) ≥ θˆ4 ( |tˆ) for ˆt ≥ 0 and 0 ≤ ˆz ≤ . (4.19) 2 2

−z 2 4

0

zˆ = 0

tˆ→0

where Dν (z) is the parabolic cylinder function, θˆ4 (ˆ z |tˆ ) is the modified theta function, Erfc(¯z) denotes the complimentary error function and S¯o is obtained from the initial condition (4.6). The parabolic cylinder function, Dν (z) is defined as a 00 solution to the Weber differential equation g (z)+(ν + 12 − 2 z 4 ) g(z) = 0. This function satisfies the recurrence relation Dν+1 (z)−zDν (z)+νDν−1 (z) = 0andreduces to spez2 p cial relations D−1 (z) = e 4 π2 Erfc √z2 and D0 (z) =

θˆ4 (ˆ z |tˆ) =

θˆ4 (ˆ z |tˆ)

  −1 y¯ = 1 τ √ θˆ4 | 2 2Re 2Re s !! r √ t¯ − τ 2(−Re Re ¯ 2 e t−τ ) − 2Re Erfc dτ π 2(¯t − τ ) r !! 23/2 ¯ −Re Re √ t e 4t¯ D−3 . (4.18) t¯ π

and therefore the modified theta function (4.17) is a nonpositive valued function. Since f (ˆ n) is a continuous function, one can show that the bounds for the summation of discrete values of f (ˆ n) exist in terms of definite integral (whose integrand is f (ˆ n)) such that Z ∞ 0≤ f (ˆ n) dˆ n ≤ θˆ4 (ˆ z , tˆ)

In other words, the exact integral in (4.15), which we refer as the convolution term, can be determined at y¯ = 1. Some specific values (for zˆ = 0 or tˆ → 0) obtained for the modified theta function using (4.17) are given in Table 1 and we note that it has a singularity at zˆ = 12 , i.e θˆ4 ( 12 |tˆ → 0) = −∞. Also, from the work that follows, we reproduce the Table 1 using the bounds and show that the convolution integral (4.18) is indeed finite in spite of the integrand being singular when τ → 0.

1

=

∞ X n ˆ =1

Z f (ˆ n) ≤

∞

f (ˆ n) dˆ n + f (1).

(4.22)

1

The definite integral in the above equation is given through      Z ∞ 1 − 2ˆz 1 + 2ˆz 1 √ √ Erf − Erf . f (ˆ n) dˆ n= 2 1 2 ˆt 2 ˆt (4.23) 6

non-negative Re, y¯ and t¯ we immediately conclude that the convolution integral in (4.15) is finite. We note that, for two specific cases, 1) when t¯ = 0 and 2) when y¯ = 0 in (4.26), the exact value of the convolution integral and its upper-bound are equal to zero, i.e.

On substituting (4.21) and (4.23) in (4.22) we obtain the bounds for the modified theta function (4.17) to be      1 1 − 2ˆz 1 + 2ˆz √ √ 0≤ Erf − Erf ≤ θˆ4 (ˆ z |tˆ) 2 ˆt ˆt 2 2      1 + 2ˆz 1 1 − 2ˆz √ √ − Erf ≤ Erf 2 2 ˆt 2 ˆt  2 ˆ 1 1  −( 1 +ˆz)2 /tˆ +√ − e−( 2 −ˆz) /t . (4.24) e 2 π tˆ

Cupper (¯ y , 0) = 0 for 0 ≤ y¯ ≤ 1 and Re ≥ 0, and Cupper (0, t¯) = 0 for t¯ ≥ 0 and Re ≥ 0. (4.27) y , t¯)) ≥ 0 and Also, one can show that ∂∂y¯ (Cupper (¯ ¯ ¯ Cupper (0, t) ≤ Cupper (¯ y , t) ≤ Cupper (1, t¯) which implies that the maximum value of Cupper (¯ y , t¯) occurs at y¯ = 1. ¯ Further, when t << Re, the maximum value Cupper (1, t¯), is very close to zero, i.e.

and the bounds are finite for 0 ≤ zˆ < 12 and tˆ ≥ 0. At zˆ = 1 ˆ 2 and t → 0, the lower bound value is −∞ (as obtained in Table 1). Therefore, the modified theta function θˆ4 (ˆ z , tˆ) is absolutely convergent for 0 ≤ zˆ ≤ 12 and tˆ > 0.
τ Referring to the equation (4.15), √−1 θˆ4 y2¯ | 2Re ≥ 0 it 2Re q q  √ −Re t¯−τ 2(t¯−τ ) Re − 2ReErfc is easy to show that 2 π e 2(¯ t−τ ) is non-negative and a monotonically decreasing function of τ . It immediately follows that, the upper and lower bounds of the convolution term can be shown to be Z t¯  y¯ τ  −1 √ 0≤ θˆ4 | 2 2Re 2Re 0 s !! r √ t¯ − τ 2(−Re Re ¯−τ ) t e dτ 2 − 2Re Erfc π 2(¯t − τ ) r r !! ¯ −Re √ Re t e 2t¯ − 2Re Erfc ≤ Cupper (¯ y , t¯) = 2 π 2¯t Z t¯ −1 ˆ  y¯ τ  √ | θ4 dτ . (4.25) 2 2Re 2Re 0

0 ≤ Cupper (¯ y , t¯) ≤ Cupper (1, t¯) → 0, ¯ when t << Re and 0 ≤ y¯ ≤ 1 .

(4.28)

From (4.19), one can show that the integrand of the convolution term is a monotonically increasing function of y¯ with its minimum and maximum value occurring at y¯ = 0 and y¯ = 1 respectively. Using (4.19), (4.18) and (4.16), the values for the bound are obtained to be Z t¯ −1 ˆ  y¯ τ  √ 0≤ | θ4 2 2Re 2Re 0 s !! r √ t¯ − τ 2(−Re Re ¯−τ ) t e dτ 2 − 2Re Erfc π 2(¯t − τ )   Z t¯ −1 ˆ 1 τ √ ≤ Cmax (t¯) = | θ4 2 2Re 2Re 0 s !! r √ Re t¯ − τ 2(−Re ¯ 2 e t−τ ) − 2Re Erfc dτ π 2(¯t − τ ) ! r ! r  −Re Re Re 2Re . = t¯ Erfc + 1 − e 2¯t ¯t 2¯t π¯t

On applying the lower-bound of the modified theta function (4.24) into (4.25), we obtain the upper bound of the convolution term, Cupper (¯ y , t¯) to be Z t¯  y¯ τ  −1 √ | θˆ4 2 2Re 2Re 0 s !! r √ t¯ − τ 2(−Re Re ¯ 2 e t−τ ) − 2Re Erfc dτ π 2(¯t − τ ) r r !! ¯ −Re √ t Re ≤ Cupper (¯ y , t¯) = 2 e 2t¯ − 2Re Erfc π 2¯t √ √ Re(1−y) ¯ 2 Re(1+y) ¯ 2 ¯ − 2t¯ (3 + y¯) t¯ e− 2t¯ (3 − y¯)  te √ √ − 2 π 2 π ! √ √ t¯ Re(1 + y¯) Re(1 − y¯)2 √ √ √ − Erfc + 2 2 2 2Re 2t¯ ! √  √ Re(1 − y¯) t¯ √ √ − 2Re + Erfc 2 2Re 2t¯ !! √ 2 √ Re(1 − y¯) √ + − 2Re (1 − y¯) (4.26) 2 2

(4.29) In eq. (4.25) and (4.29), we have obtained two upper bounds and therefore, we use the lower of the two as the upper bound for the convolution term. This implies that the bounds of the convolution term are given through t¯

−1 ˆ  y¯ τ  √ θ4 | 2 2Re 2Re 0 s !! r √ t¯ − τ 2(−Re Re ¯−τ ) t 2 e dτ − 2Re Erfc π 2(¯t − τ )  ≤ min Cupper (¯ y , t¯), Cmax (t¯) (4.30) Z

0≤

where 0 ≥ y¯ ≥ 1, t¯ ≥ 0 and Re ≥ 0. On using (4.30), (4.26), (4.29) in (4.15), we get the bounds

for 0 ≥ y¯ ≥ 1, Re ≥ 0 and t¯ ≥ 0. In view of the upperbound for the integral, Cupper (¯ y , t¯), which is finite for any 7

for the dimensionless stress, for λo ≥ 0, to be ! r ! r  −¯ y2 Re Re y ¯2 Re 2Re ¯ ¯ 2 t λo t Erfc y ¯ +1 −y ¯e + S¯o ¯t 2¯t π¯t r !  Re y ¯2 Re ¯ y , t¯) ≥ λo t¯ Erfc y ≥ S(¯ ¯ +1 ¯t 2¯t ! r  −y ¯2 Re 2Re − λo min Cupper (¯ y , t¯), Cmax (t¯) + S¯o − y¯e 2t¯ ¯ πt

n =0, n =0, n =0, n =0, n =0, n =0,

Dimensionless stress, S¯

1.2

(4.31)

Re =1, Re =1, Re =1, Re =5, Re =5, Re =5,

t¯=0.3 t¯=0.3 t¯=0.3 t¯=1.4 t¯=1.4 t¯=1.4

(Upperbound) (Eq. (4.15)) (Lowerbound) (Upperbound) (Eq. (4.15)) (Lowerbound)

1

0.8

0.6

0.4

where 0 ≥ y¯ ≥ 1, Re ≥ 0 and t¯ ≥ 0. We note that on applying (4.27) in (4.31), the solution satisfies the initial and the boundary conditions at y¯ = 0 exactly (the bounds are equal). And, the boundary condition at y¯ = 1 lies within the bounds. However, referring to the equation (4.28), we note that, for a special case where t << Re in (4.31), the difference between the bounds of the dimensionless stress ¯ y , t¯) diminishes (as the convolution term Cupper → 0) S(¯ ¯ y , t¯) to good accuracy using its and one can compute S(¯ upper-bound relation, i.e. for t << Re, r !  Re y ¯2 Re ¯ y , t¯) ≈λo t¯ Erfc y +1 S(¯ ¯ ¯t 2¯t !! r −y ¯2 Re 2Re − y¯ e 2t¯ + S¯o . (4.32) π t¯

0.2

0

0.2

0.4 0.6 Dimensionless distance, y¯

0.8

Figure 1: Plots of the bounds for the dimensionless shear stress and its numerical evaluation for n = 0, λo = 1 when time t¯ is small compared to Re

numerical integration procedure. We note that when t¯ is comparable to Re, the interval between the bounds gets wider and it cannot serve as evaluator for numerical integration. In Fig. 1, we observe that the interval of bounds t¯ widens when Re increases. In Fig. 2, we plot (for n = 0), the numerical evaluations of dimensionless stress equation (4.15) and their numerical counterparts from COMSOLr where in Fig. 2a, the stress is plotted for Re = 1 and three different t¯, and in Fig. 2b we display the stress for t¯ = 1.4 and two different Re. We note that the numerical solutions from COMSOLr matches extremely well with those obtained from numerical evaluation of solution (4.15). However, when t¯ is small compared to Re, the convolution integral in solution (4.15) is of extremely small value which results in equal bounds and the solution (4.15) reduces to t¯ (4.32). Therefore for small Re , the solution (4.32) is plotted in Fig. 2a and the corresponding numerical results match very well with it.

Also, when t¯ << Re in (4.31), the boundary condition at y¯ = 1 is satisfied exactly as both the bounds approach equal value of So . 5. Numerical solutions We obtain numerical solutions for the dimensionless stress and velocity by solving the two PDEs (3.11) and (3.14) simultaneously that satisfy the conditions (4.1) through (4.3) using COMSOLr , for various n (for n < 0, n = 0 and n > 0) and Re so as to study the influence of these parameters on the unsteady solutions for the variant of Stokes’ first problem. Also, the numerical evaluation of the close-form solution (4.15) for the dimensionless stress for n = 0 and the approximate solution (5.8) for n < 0 are obtained for comparison with their numerical counterparts wherever applicable. We note that any exact or approximate solution for the dimensionless velocity involves integration (numerical) of the stress solutions and are not included in the present work. In order to validate the numerical integration scheme in evaluating (4.15), we compared the numerical calculations against the bounds (4.31) that we established. We find that the numerical evaluation lies within the bound and in Figure 1, we plot them for small time t¯ for different Re. What is important is that at early times, the lower and upper bounds are very tight and lie so close to one another that they serve as very good evaluators for the

5.1. Numerical solutions corresponding to PDEs (3.11) and (3.14) For the variant of Stokes’ first problem, the numerical unsteady solutions for dimensionless shear stress and dimension velocity are obtained (by solving PDEs (3.11) and (3.14) satisfying conditions (4.1) through (4.3) using COMSOLr ) for three different values of n: n < 0, n = 0 and n > 0 and for different values of Re. 5.1.1. Solutions for varying n and same Re In Fig. 3a, the variation of the numerical solutions for stress obtained for the three different values of n at two different instants of time for a fixed Re are plotted. Also, for n = 0, the plots of the numerical approximation for the t¯ stress obtained from Eq. (4.32) for small Re and from Eq. (4.15) for all other cases are included in the Fig. 3a. And, 8

1 0.9 0.8

0.4

0.7 0.6 0.5 0.4 0.3

0.35 0.3 0.25 0.2 0.15

0.2

0.1

0.1

0.05

0

0.2

0.4 0.6 Dimensionless distance, y¯

n =-1, t¯=0.2 (Numerical) n =0, t¯=0.2 (Eq. (4.32)) n =0, t¯=0.2 (Numerical) n =1, t¯=0.2 (Numerical) n =-1, t¯=0.5 (Numerical) n =0, t¯=0.5 (Eq. (4.15)) n =0, t¯=0.5 (Numerical) n =1, t¯=0.5 (Numerical)

0.45

Dimensionless stress, S¯

Dimensionless stress, S¯

0.5

t¯=0.2 (Numerical) t¯=0.2 (Eq. (4.32)) t¯=0.5 (Numerical) t¯=0.5 (Eq. (4.15)) t¯=1 (Numerical) t¯=1 (Eq. (4.15))

0.8

0

0.2

(a)

1.4

0.8

(a)

¯ Re=1 ¯ Re=1 ¯ Re=5 ¯ Re=5

1.2

0.4 0.6 Dimensionless distance, y¯

(Numerical) (Eq. (4.15)) (Numerical) (Eq. (4.15))

−0.05

Dimensionless velocity, S¯

Dimensionless stress, S¯

−0.1 1

0.8

0.6

0.4

−0.15 −0.2 −0.25 −0.3 n =-1, t¯=0.2 n =0, t¯=0.2 n =1, t¯=0.2 n =-1, t¯=0.5 n =0, t¯=0.5 n =1, t¯=0.5

−0.35 0.2 −0.4 0

0.2

0.4 0.6 Dimensionless distance, y¯

0.8

0

0.2

0.4 0.6 Dimensionless distance, y¯

0.8

(b)

(b)

Figure 2: Plots of numerical solutions for the dimensionless shear stress for n = 0, λo = 1 when a) Re = 1 and t¯ = [0.2, 0.5, 1] and b) Re = [1, 5] and t¯ = 1.4

Figure 3: Plots of numerical solutions (corresponding to PDEs (3.11) and (3.14)) for a) dimensionless stress and b)dimensionless velocity obtained for n = [1, 0, −1], λo = 1, Re = 1 and t¯ = [0.2, 0.5]

in Fig 3b, we plot the numerical solutions for the velocity distribution corresponding to the numerical solutions for stress plotted in Fig 3a. From Fig 3a, we note that the disturbance set in by the stress boundary condition at the boundary y¯ = 0 propagates but has not reached the boundary y¯ = 1 for short times. The undisturbed region, where the stress and its gradient are zero, shrinks with increasing time. At a given time, the stress at the boundary y¯ = 0, for different fluids of varying n, is the same but their corresponding velocity gradients, observed from Fig 3b, varies reflecting their stress thinning/thickening behavior. We observe that for the same magnitude of stress at the boundary y¯ = 0, the velocity gradient for the stress thinning fluid (n > 0) is higher than that for linear vis-

cous fluid (n = 0) while that for the stress thickening fluid (n < 0) is lower when compared to that for linear viscous fluid. This disparity between the classical linearly viscous fluid and fluids whose exponent n is not zero, becomes more significant with increase in time. 5.1.2. Effect of Re on solutions for varying n at a fixed t¯ The effect of Re on the solutions at a fixed time are portrayed in Figures 4a and 4b for the dimensionless shear stress and velocity respectively, for three different n and two different Re values. The dimensionless stress in the fluid, irrespective of n, at a given distance y¯ (where 0 < y¯ < 1) and time t¯ > 0 is such that the magnitude of the dimensionless stress increases as Re decreases. At a given time, the disturbance set in by the stress boundary 9

solution to (3.28), i.e. 1.2

n =-1, Re=1 (Numerical) n =0, Re=1 (Eq. (4.15)) n =0, Re=1 (Numerical) n =1, Re=1 (Numerical) n =-1, Re=5 (Numerical) n =0, Re=5 (Eq. (4.15)) n =0, Re=5 (Numerical) n =1, Re=5 (Numerical)

Dimensionless stress, S¯

1

0.8

1 ∂ Sˆ = ∂ t¯ 2Re

!

1 ˆ 1 + 2n|S|

∂ 2 Sˆ . ∂ y¯2

(5.1)

ˆ ˆ and further satisfy the criterion (3.23), i.e. 1  |S|+|n|| S|. ˆ ¯ We refer to S(¯ y , t) as an approximate solution to (3.16). In (5.1), we let

0.6

ˆ = q¯(¯ 1 + 2n|S| y , t¯) .

(5.2)

0.4

and transform the PDE (5.1) to 0.2

0

∂ q¯ 1 ∂ 2 q¯ . = ∂ t¯ 2Re q¯ ∂ y¯2 0.2

0.4 0.6 Dimensionless distance, y¯

0.8

(5.3)

An exact, multiplicative separable solution for the nonlinear PDE (5.3) exists and is of the form

(a)

ˆ fˆ(¯ q¯(¯ y , t¯) = (λt¯ + β) y)

(5.4)

where λ,βˆ are arbitrary constants and fˆ(¯ y ) is a solution to the ordinary differential equation

Dimensionless velocity, S¯

−0.5

∂ 2 fˆ = A¯ y2 ∂ y¯2

−1

where A = 2λRe. An implicit solution for (5.5) involves an improper integral and is given through Z 1 y¯ = ± q dfˆ + C2 (5.6) 2Afˆ3 + C 1 3

−1.5

−2 n =-1, Re=1 n =0, Re=1 n =1, Re=1 n =-1, Re=5 n =0, Re=5 n =1, Re=5

−2.5

−3 0

0.2

0.4 0.6 Dimensionless distance, y¯

(5.5)

where C1 and C2 are integration constants to be obtained from the boundary conditions. The improper integral in the above solution is evaluated using the hypergeometric function, and the implicit solution (5.6), valid for a finite domain, reduces to q 3 fˆ3 1 1 4 −2Afˆ fˆ 2A 3C1 + 1 F ([ 3 , 2 ], 3 , 3C1 ) q y¯ = C2 ± (5.7) 2Afˆ3 + C 1 3

0.8

(b) Figure 4: Plots of numerical solutions (corresponding to PDEs (3.11) and (3.14)) for a) dimensionless stress and b)dimensionless velocity obtained for n = [1, 0, −1], λo = 1, t¯ = 1.2 and Re = [1, 5]

where A = 2λRe and F (.) denotes the Gaussian hypergefˆ3 ometric function which is always convergent for | −2A 3C1 | ≤ 1. On substituting (5.4) in (5.2) and letting βˆ = 0 in it, we obtain an approximate solution of the form

condition (at y¯ = 0) in a fluid propagates through the distance such that the region of disturbance shrinks as Re increases. From Fig 4a and 4b, we observe that, for the same magnitude of dimensionless stress at the boundary y¯ = 0, the resulting velocity gradients increase as Re increases irrespective of the value of n. The influence of Re on the stress thinning and thickening behavior of the fluids is more significant for higher values of Re.

ˆ y , t¯)| = |S(¯

y) ¯ −1 λfˆ(¯ + t 2n 2n

(5.8)

where λ is an arbitrary constant and fˆ(¯ y ) satisfies (5.7). The initial condition satisfied by the approximate solution (5.8) is given by

5.2. Approximate solution for n < 0 We now derive an approximate solution for the variant of Stokes’ first problem for the stress power-law fluid. In addition to determining the nature of approximate solution for n < 0, we use it to check the numerical counterˆ y , t¯) be a parts obtained for different values of n. Let S(¯

ˆ y , 0)| = |S(¯

−1 . 2n

(5.9)

which implies that n < 0. Therefore to obtain a constant stress condition at the initial time, we must restrict the 10

approximate solution (5.8) to n < 0. For a finite domain where y¯ ∈ (0, 1), to obtain a constant stress condition at y¯ = 1, we let fˆ(1) = 0 in (5.8) which implies that the boundary condition at y¯ = 1 is given by

0.09 0.08 ¯ Dimensionless stress, |S|

ˆ t¯)| = −1 |S(1, 2n

0.1

(5.10)

where t > 0. Using fˆ(1) = 0 in (5.7) implies that C2 = 1 and y¯ = 1 ±

3 4 −2Afˆ 3 , 3C1 )

fˆ F ([ 31 , 12 ], √ C1

0.07 0.06 0.05 0.04 0.03 0.02

(5.11)

t¯=0.5 (Eq.(4.40)) t¯=0.5 (Numerical) t¯=1 (Eq.(4.40)) t¯=1 (Numerical)

0.01 0 0

where C1 is an arbitrary constant to be obtained from another boundary condition. We note that C1 is fixed for a given fˆ(0), λ and Re, and we choose to use the lower sign form of the above equation in the present work. The boundary condition satisfied by the approximate solution (5.8) at y¯ = 0, for n < 0 is given by ! ˆ(0) λ f −1 ˆ t¯)| = + t¯ . (5.12) |S(0, 2n 2n

0.2

0.4 0.6 0.8 Dimensionless distance, y¯

1

(a)

0.1

t¯=0.5 (Eq.(4.40)) t¯=0.5 (Numerical) t¯=1 (Eq.(4.40)) t¯=1 (Numerical)

0.09

¯ Dimensionless stress, |S|

0.08

The above initial and boundary conditions (5.9), (5.10) and (5.12) are similar to those obtained for the variant of Stokes’ first problem in (4.1) through (4.3) except for their dependence on parameter n and the requirement that n < 0. For a fixed n < 0, the initial condition and the boundary condition at y¯ = 1 imply that the stress is a non-zero constant unlike the original Stokes’ first problem where the stress is zero at initial time and zero at the farthest boundary at all time.

0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0

0.2

0.4 0.6 0.8 Dimensionless distance, y¯

1

(b)

5.3. Comparison of the approximate solution (5.8) and numerical solution to PDE (3.15) for n < 0 We now compare the numerical solution for dimension¯ y , t¯) (corresponding to PDE (3.16) that satless stress S(¯ isfies the initial and boundary conditions (4.1) through ˆ y , t¯) in (5.8) for the (4.3)) with the approximate solution S(¯ ¯ ¯ particular case when 1  |S| + |n||S|, n < 0 and S¯ > 0. The initial and boundary conditions (4.1) through (4.3) are equivalent to (5.9), (5.10) and (5.12) when S¯o = −1 2n

0.1

¯ Dimensionless stress, |S|

0.08

ˆ

and λo = λf2n(0) . Therefore, for n < 0, the condition λo < 0 ensures that the initial dimensionless stress S¯ reduces with dimensionless time but is positive up to some limited dimensionless time. In Figure 5a and 5b, we plot the variation of dimensionless stress with distance, obtained by the numerical and approximate methods, at two different time instants when Re is fixed. As expected, the fluid with an initial positive shear stress (S¯o = −1 2n ), when subjected to a negative, con¯

Re=1 Re=1 Re=5 Re=5

0.09

(Eq.(4.40)) (Numerical) (Eq.(4.40)) (Numerical)

0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0

0.2

0.4 0.6 0.8 Dimensionless distance, y¯

1

(c) Figure 5: Plots of approximate and the numerical solutions of the dimensionless shear stress for 2nλo =0.93361 when a) n = −5, Re = 1 at t¯ = [0.5, 1] b) n = −10, Re = 1 at t¯ = [0.5, 1] and c) n = −10 at t¯ = 1 for Re = [1, 5]

ˆ

stant time rate of shear stress ( ∂∂St¯ |y¯=0 = λ0 = λf2n(0) < 0 ), shows a decrease in the magnitude of stress with increasing time. Of course, this reduction in magnitude occurs till the stress reaches zero, beyond which, the shear stress can flip sign and its absolute value can increase. It can 11

be observed that for a given n, at a given distance y¯ < 1, the stress increases linearly with time and the time rate of increase of dimensionless shear stress depends on the dimensionless distance (see equation (5.8)). We observe that the approximate solution is in good agreement with the corresponding numerical solution and the difference diminishes for smaller magnitude of stress. We also study the nature of the solutions at some fixed λo < 0 (we let λ = 10 and 2nλo =0.93361) and t¯ > 0 but for different values of Re. Figure 5c shows the influence of the parameter Re on the dimensionless stress with distance at a fixed time. When two different fluids with same n and different Re are subjected to the same initial and stress boundary conditions, it is observed that at a fixed distance and time, the difference in the shear stress from its initial value is higher for lower values of Re. We note that the initial stress is dependent on n and therefore to compare the behavior of stress power-law fluids for various n, one ¯ can define a scaled dimensionless stress, nS.

References [1] Perlacova, T., and Prusa, V., Tensorial implicit constitutive relations in mechanics of incompressible non-Newtonian fluids, J. Non-Newton. Fluid Mech., 216:13-21, (2015). doi: 10.1016/j.jnnfm.2014.12.006 [2] J. Malek and K. R. Rajagopal, Mathematical issues concerning the Navier-Stokes equations and some of their generalizations, in: Handbook of Evolutionary Equations, Vol II (eds. C. Dafermos and E. Feireisl), Elsevier (2005). [3] R. R. Huilgol, Mechanics of viscoelastic fluids, Hindustan Publishing Corporation, (1975). [4] W. R. Schowalter, Mechanics of non-Newtonian Fluids, Pergamon Press, Oxford (1978). [5] K. R. Rajagopal, On Implicit constitutive Theories, Application of Mathematics, 28, no. 4, 279-319 (2003). [6] K. R. Rajagopal, On implicit constitutive theories for fluids, Journal of Fluid Mechanics, 550, 243-249, (2006). [7] V. Prusa and K. R. Rajagopal, On implicit constitutive relations for materials with fading memory, Journal of Non-Newtonian Fluid Mechanics, 181-182, 22-29 (2012). [8] J. Malek, V. Prusa and K. R. Rajagopal, Generalizations of the Navier-Stokes fluid from a new perspective, International Journal of Engineering Science, 48, 1907-1924 (2010). [9] Le Roux,C., Rajagopal, K. R., Shear flows of a new class of power-law fluids, Appl. Math. 58 (2) (2013) 153177, doi: 10.100 7/s10492-013-0008-4. [10] S. P. A. Narayanan and K. R. Rajagopal, Unsteady flows of a class of novel generalizations of the Navier-Stokes fluid, Applied Mathematics and Computation, 219, 9935-9946 (2013). [11] G.G. Stokes, On the effect of the internal friction of fluids on the motion of pendulums, Transactions of the Cambridge Philosophical Society 9 (1851) 8-106. [12] Lord Rayleigh, On the motion of solid bodies through viscous liquid, Philosophical Magazine 6 (21) (1911) 697711. [13] Fritz Oberhettinger and Larry Badii, Tables of Laplace Transforms, Springer-Verlag, (1973).

6. Conclusions In (4.18), we have an exact integral wherein the integrand is singular when τ → 0 and is equal to zero when τ = t¯. The variation of the integrand is such that it is very sharp near τ = 0 as it varies from ∞ at τ = 0 to very small finite values adjacent to τ = 0. One can make use of the exact integral available in (4.18) to validate/enhance the performance of quadrature techniques (in terms of convergence and accuracy) in the numerical integration of this integral. Similarly, in (4.30), when time is small, we have an extremely small interval between the bounds for the integral, and one can use this to validate/enhance the numerit¯ is small, the cal integration techniques. In (4.31), when Re bounds for the solution are extremely close to each other and therefore one can use this result to test the numerical schemes in solving the PDE (3.15) for n = 0 satisfying the conditions (4.1) through (4.3). And in (5.8), we have an approximate solution to the non-linear PDE (3.15) for n < 1 which can be used as a guideline to test the numerical scheme used in solving the PDE as we have shown that the numerical solution is closer to the approximate one whenever the condition (3.23) is satisfied (see Fig. 5).

12