Start-up plane Poiseuille flow of a Bingham fluid

Journal of Non-Newtonian Fluid Mechanics 265 (2019) 133–139 Contents lists available at ScienceDirect Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm Short Communication Start-up plane Poiseuille flow of a Bingham fluid Raja R. Huilgol a,∗, Andreas N. Alexandrou b,†, Georgios C. Georgiou c a College of Science and Engineering, Flinders University, GPO Box 2100, Adelaide SA 5001, Australia Department of Mechanical and Manufacturing Engineering, University of Cyprus, PO Box 20537, Nicosia, 1678 Cyprus c Department of Mathematics, University of Cyprus, PO Box 20537, Nicosia, 1678, Cyprus b a r t i c l e i n f o a b s t r a c t Keywords: Bingham plastic plane Poiseuille flow start-up flow analytical solution The start-up flow of a Bingham plastic in a channel is considered and Safronchik’s solution [1] for the initial evolution of the yield surface and the core velocity is revisited. Stricter time bounds for the validity of the above solution are derived and the solution is extended to include the velocity profile in the evolving yielded zone. Comparisons are made with another approximate solution derived under the assumption that the velocity in the yielded zone is parabolic adjusting with the evolving yield surface. This approximation performs well for small values of the yield stress, or, equivalently, for large values of the imposed pressure gradient. 1. Introduction solution by Safronchik [1] holds over a period of time interval less than that found by him. Next in Section 2, using the Safronchik approximation, the velocity field 𝑢 = 𝑢(𝑦, 𝑡) in the yielded region can be found; see Eq. (2.26) below. This result complements the solution found in [1] and [5], where only the velocity in the core was given. It is worth noting that obtaining the analytical solution for 𝑢 = 𝑢(𝑦, 𝑡) takes considerable effort. In this connection, one notes that three integrals, 𝐼𝑗 (𝑦, 𝑡), 𝑗 = 1, 2, 3, appear in the solution; see Eq. (2.24). It turns out that the last two are respectively one and two orders of magnitude less than that of the first one and can be ignored as long as the yield surface is close to the boundary of the channel. See Fig. 3, where these integrals are depicted when the Bingham number Bn = 1 and 𝐺 = 1.2. In Section 3, a comparison of the velocity in the core uc (t) given by Eq. (2.18) with that obtained from u(𝛿(t), t) in Eq. (2.26) is made when the Bingham number Bn = 1 and the pressure drop G ∈ {5, 10, 100}. It is found that these two values agree with one another as G increases. Since the derivation of the velocity field in the yielded region is complicated, one may assume that u(y, t) is parabolic and varies with time as the size of the core decreases; see Eq. (3.2). In Figs. 5 and 6, we have compared the evolution of the velocity profile in Eq. (2.26) with that in Eq. (3.2) over the time interval [0, t3 ]. Once again, the Bingham number Bn = 1, while G ∈ {10, 100}. It is found that the two profiles are similar for small t, when 𝐺 = 10, and diverge as t → t3 . The differences are greater when 𝐺 = 100. It is clear that the comparisons in Figs. 5 and 6 are not accurate, for the velocity field u(y, t) in Eq. (2.26) is based on the assumption that 𝛿(𝑡) = 1. Therefore, in Figs. 7–9, we compare the velocity profile In non-Newtonian fluid mechanics, analytical solutions of initialboundary value problems are very rare. When one considers the flows of a Bingham fluid, only four are known. Three of them were found by Safronchik [1–3], and the other by Sekimoto [4]. The problem of interest here is the solution to the start-up flow in a channel, or the plane Poiseuille flow [1]. It is assumed that at 𝑡 = 0+ , a constant nondimensional pressure drop per unit length G > 0 is suddenly applied and the Bingham fluid is set in motion. The yield surface moves into the fluid from the upper and lower planes of the channel and a core forms in the interior. This core region continues to shrink until the flow becomes steady. The solution to the problem requires both the velocity field as well as the location of the yield surface 𝛿(t) be found at any given time t. Safronchik’s method [1] delivers these results for a finite period of time only due to the way the solution has been constructed; for a detailed description of the method and its limitations, see Huilgol [5]. Since the location of the yield surface at 𝑦 = 𝛿(𝑡) has to be found from the solution of a nonlinear integral equation, an approximate solution valid for a short period of time only can be derived. In Section 2, Safronchik’s method [1] is revisited and the result for 𝛿(t) is given by Eq. (2.15), and the velocity uc (t) of the core appears in Eq. (2.18) and its time of validity t1 is given by Eq. (2.19). If one assumes that the velocity uc (t) in the core is an increasing function of time till the flows becomes steady, a different upper bound t2 for the period of approximation arises; see Eq. (2.20). Finally, a third bound t3 can be derived if the size of the core is assumed to shrink and approach its final value when the flow is steady; see Eq. (2.23). It is found that t3 < t2 < t1 , which means that the ∗ † Corresponding author. E-mail address: raj.huilgol@flinders.edu.au (R.R. Huilgol). Professor Alexandrou has passed away. https://doi.org/10.1016/j.jnnfm.2018.10.009 Received 26 July 2018; Received in revised form 15 October 2018; Accepted 26 October 2018 Available online 27 October 2018 0377-0257/© 2018 Elsevier B.V. All rights reserved. R.R. Huilgol, A.N. Alexandrou and G.C. Georgiou Journal of Non-Newtonian Fluid Mechanics 265 (2019) 133–139 predicted by Eq. (2.26) with the parabolic approximation of the velocity field in Eq. (3.2) over the time interval [0, t09 ], where t09 is the time interval when 𝛿(𝑡09 ) = 0.9. It is seen that the parabolic velocity approximation is adequate for small Bn/G ratios. The Appendix describes the three integrals which appear in the velocity field in the yielded region; see Eq. (2.24). It is found that one of them can be evaluated analytically, while the other two cannot be. 2. Formulation and solution of the problem 2.1. The problem We assume that an incompressible Bingham fluid occupies a channel of width 2H in the 𝑦−direction. The flow occurs in the x-direction and is symmetric about the centreline with a velocity field 𝑢 = 𝑢(𝑦, 𝑡). When the fluid has yielded, the constitutive equation for the shear stress 𝜏 in the Bingham fluid is given by 𝜏 = 𝜏𝑦 + 𝜂 𝜕𝑢 , 𝜕𝑦 (2.1) Fig. 1. The parameter 𝛼 as a function of the ratio Bn/G. where 𝜏 y is the yield stress and 𝜂 is the plastic viscosity of the fluid. In the core where the fluid moves as a rigid body, the shear stress 𝜏 is not defined by a constitutive equation; rather it is bounded by the yield stress. That is |𝜏| < 𝜏 y . We consider the start-up flow of the Bingham fluid in the channel initiated by a constant pressure gradient G > 𝜏 y /H. That is, 𝜕 𝑝∕𝜕 𝑥 = −𝐺. The flow is divided into a yielded region (h(t), H] and an unyielded one [0, h(t)]. In the latter, the central core moves with an unknown velocity uc (t). Balance of linear momentum leads to the following equation: 0 ≤ 𝑦 ≤ 𝐻, 𝑡 > 0, 𝜕𝑢 𝜕𝜏 𝜌 = + 𝐺, 𝜕𝑡 𝜕𝑦 Let the Bingham number be defined through Bn = 𝜏 y H/𝜂U. Hence, the right side becomes ∫0 𝑡 𝑢𝑐 (𝑡) = 𝐺𝑡 − Bn 𝜕𝑢 𝜕2 𝑢 = − 𝐺, 𝜕𝑡 𝜕𝑦2 𝑢(𝑦, 0) = 0, (2.5) 𝜂 𝑡, 𝜌𝐻 2 𝛿(𝑡̃) = 𝐻 𝜏̃ = 𝜏, 𝜂𝑈 ℎ(𝑡) , 𝐻 𝐺𝑡 1 𝜂𝑈 ̃ 𝐺⋅ = ⋅ 𝜌 𝜌 𝐻2 𝜂 𝑡 > 0. (2.13) (2.14) Here the constant 𝛼 is the solution of the equation (2.6) 𝐻2 𝐺̃ = 𝐺. 𝜂𝑈 𝑡̃ = 𝑈 𝐺̃ 𝑡̃. − 2𝛼 ∫𝛼 ∞ 2 𝑒−𝑠 𝑑𝑠 = Bn , 𝐺 (2.16) (2.17) where erf (⋅) is the error function. In Fig. 1, we have depicted the way 𝛼 varies with the ratio Bn/G. The velocity in the core is given by [1] [ √] 1√ 1 ln 1 − 2𝛼 𝑡 . 𝑢𝑐 (𝑡) = 𝐺𝑡 + Bn 𝑡 + Bn (2.18) 𝛼 4𝛼 2 (2.8) ̃ leading to Next in analogy with 𝑡 = (𝜌𝐻 2 ∕𝜂)𝑡̃, we let 𝜉 = (𝜌𝐻 2 ∕𝜂)𝜉, 𝑡 𝜏 𝑡̃ 𝑡̃ 𝜏 𝐻 𝜏𝑦 𝜌𝐻 2 ̃ 𝑦 𝑦 1 1 ̃ 𝑑𝜉 = ⋅ 𝑑 𝜉. 𝑑𝜉 = ̃ ̃ ∫0 𝜂𝛿(𝜉) 𝜌 ∫0 ℎ(𝜉) 𝜌 ∫0 𝐻𝛿(𝜉) 𝜂 −𝛼 2 which can also be written as follows: [ ] √ 2 Bn 𝑒−𝛼 − 𝛼 𝜋 1 − erf (𝛼) = , 𝐺 (2.7) Thus, 𝜌𝐻 2 𝑡 ≥ 0, In his important work, Safronchik [1] found that the location of the yield surface in the upper-half of the channel is given in a nondimensional form by √ 1 𝛿(𝑡) = 1 − 2𝛼 𝑡, 0 ≤ 𝑡 ≤ . (2.15) 4𝛼 2 𝑒 𝑢 𝑢̃ = , 𝑈 𝑢(1, 𝑡) = 0, 2.3. Safronchik’s solution for the core and time of validity We shall now introduce the following non-dimensional variables, us√ ing H as the characteristic length scale and 𝑈 = 𝜏𝑦 ∕𝜌 as the characteristic velocity scale: 𝑡̃ = (2.12) The required equations have now been assembled. These are Eqs. (2.12)–(2.14); their solution provides the velocity in the core (2.11), valid in 0 ≤ y ≤ 𝛿(t). 2.2. Non-Dimensionlisation 𝑦 , 𝐻 𝛿(𝑡) < 𝑦 < 1, 𝑡 > 0, 0 ≤ 𝑦 ≤ 1; 𝜕𝑢 (𝛿(𝑡), 𝑡) = 0, 𝜕𝑦 Of course, 𝑢(𝑦, 𝑡) = 𝑢𝑐 (𝑡), 0 ≤ 𝑦 ≤ ℎ(𝑡). 𝑦̃ = (2.11) where the tildes have again been dropped. The following conditions apply: In this core, the velocity 𝑢 = 𝑢𝑐 (𝑡) only, whence the equation of motion (2.2) reduces to [ ] 𝜏𝑦 𝑑𝑢𝑐 1 = 𝐺− . (2.4) 𝑑𝑡 𝜌 ℎ(𝑡) 𝑥 , 𝐻 1 𝑑𝜉. ∫0 𝛿(𝜉) In the yielded region, the partial differential Eq. (2.2) has the form: where 𝜌 is the density of the fluid. In the core, while the shear stress does not satisfy a constitutive relation, its distribution is linear. That is 𝜏𝑦 𝑦, 0 ≤ 𝑦 ≤ ℎ(𝑡). (2.3) 𝜏=− ℎ(𝑡) 𝑥̃ = (2.10) Since 𝑢 = 𝑈 𝑢̃ , dropping the tildes for simplicity, the velocity in the core given by Eq. (2.5) takes on the following form: (2.2) On integration, one finds that ] [ 𝑡 𝜏 𝑦 1 𝑑𝜉 . 𝐺𝑡 − 𝑢𝑐 (𝑡) = ∫0 ℎ(𝜉) 𝜌 𝑡̃ 𝜏𝑦 𝐻 1 ̃ 𝑑 𝜉̃ = 𝑈 Bn 𝑑 𝜉. ̃ ̃ ∫0 𝛿(𝜉) 𝜂𝛿(𝜉) 𝑡̃ (2.9) 134 R.R. Huilgol, A.N. Alexandrou and G.C. Georgiou Journal of Non-Newtonian Fluid Mechanics 265 (2019) 133–139 Fig. 3. The plots of 𝐼𝑖 (𝛿(𝑡), 𝑡), 𝑖 = 1, 2, 3 for Bn = 1 and 𝐺 = 1.2. For this choice of parameters, 𝛼 = 0.09962, 𝑡3 = 0.6997, and 𝛿(𝑡3 ) = 0.8333. Fig. 2. The three critical times t1 , t2 and t3 versus the ratio Bn/G. We see that the solutions for 𝛿(t) and the velocity field 𝑢 = 𝑢(𝑦, 𝑡) are valid provided t < t1 , where In Fig. 3, noting that the integrals are all negative, the absolute values of 𝐼𝑗 (𝛿(𝑡), 𝑡), 𝑗 = 1, 2, 3, for Bn = 1, 𝐺 = 1.2 have been plotted versus time t up to 𝑡3 = 0.6997. For this choice of Bn and G, it is clear that −𝐼2 and −𝐼3 are respectively one and two orders of magnitude less than −𝐼1 . Hence, the velocity field in the yielded region is given by the following approximation: 1 . (2.19) 4𝛼 2 Next, the velocity in the core, uc (t), must be an increasing function of time till the flow becomes steady. In Eq. (2.18), by putting 𝑑 𝑢𝑐 (𝑡)∕𝑑 𝑡 = 0, we find that this velocity increases provided t ≤ t2 , where )2 (√ 1 + Bn2 ∕𝐺2 + 1 − Bn∕𝐺 . (2.20) 𝑡2 = 16𝛼 2 From Eq. (2.15), we find that ] [ √ 1 Bn (2.21) 𝛿(𝑡2 ) = − 1 + Bn2 ∕𝐺2 . 1+ 2 𝐺 𝑡1 = 𝑢(𝑦, 𝑡) = 𝐺𝑡 + 𝐼1 (𝑦, 𝑡), 𝑢(𝑦, 𝑡) = 𝐺 Bn . (2.22) 𝐺 Since 𝛿(t) in Eq. (2.15) decreases with time, it will attain the value when 𝑡 = 𝑡3 given by 𝛿∞ = (2.25) ∫0 𝑡 erf(𝛽(𝑦, 𝜎)) 𝑑𝜎, 𝛿(𝑡) ≤ 𝑦 ≤ 1, 𝑡 > 0, (2.26) From the above, it follows that the core velocity u(𝛿(t), t) is given by 2 2𝛼 𝑢(𝛿(𝑡), 𝑡) = 𝐺𝑡[erf(𝛼) + √ 𝑒−𝛼 − 2𝛼 2 erfc(𝛼)], 𝜋 (1 − Bn∕𝐺)2 . (2.23) 𝑡3 = 4𝛼 2 In Fig. 2, we have compared the three critical times t1 , t2 and t3 as they vary with ratio Bn/G. It is apparent that t3 < t2 < t1 . Thus, the velocity in the core holds over a time interval less than t1 found by Safronchik [1]. (2.27) where erfc( · ) is the complementary error function. That is, the velocity in the core is proportional to Gt and does not agree with that given by uc (t) in Eq. (2.18). The reason lies in the various approximations made to arrive at these two values. In Fig. 4, we have compared the two when the Bingham number Bn = 1, and the pressure drop G ∈ {5, 10, 100}. √ When G is small, the profile of uc (t) given by Eq. (2.18) depends on 𝑡 and is curved, and at large values of G, the linear term is dominant and the profile of uc (t) is almost linear. In this situation, the core velocities in Eq. (2.18) and (2.42) are almost identical. 2.4. Velocity field in the yielded region We now turn to the determination of the velocity field in the yielded region. This is given by Eqs. (6.1.19) and (6.1.34) in [5], and takes the following form: 𝑢(𝑦, 𝑡) = 𝐺𝑡 + 𝐼1 (𝑦, 𝑡) + 𝐼2 (𝑦, 𝑡) + 𝐼3 (𝑦, 𝑡). 𝑡 > 0. ∞ √ 2𝑎 √ −𝑎2 ∕𝑡 4𝑎 3 −𝑎2 𝑧2 𝑑𝑧] = 𝐺[𝑡 erf(𝑎∕ 𝑡) + √ 𝑡𝑒 −√ √ 𝑒 𝜋 𝜋 ∫1∕ 𝑡 √ (1 − 𝑦) √ [−(1−𝑦)2 ∕4𝑡] = 𝐺[𝑡 erf[(1 − 𝑦)∕2 𝑡] + √ 𝑡𝑒 ] 𝜋 √ (1 − 𝑦)2 −𝐺 (1 − erf[(1 − 𝑦)∕2 𝑡)]). 2 Next, the width of the rigid core when the flow is fully developed is given by 𝛿∞ 𝛿(𝑡) ≤ 𝑦 ≤ 1, Consequently, from Eqs. (A.1) and (A.16) in the Appendix, we obtain 3. Numerical comparison (2.24) Since the velocity field in Eq. (2.26) is difficult to find, one may assume that the velocity field in the yielded region is parabolic, with the size of the core, 𝛿(t), and the velocity in the core, uc (t), given by Eq. (2.15) and (2.18) respectively. That is, the parabolic approximation has the form Here, the 𝐼𝑗 , 𝑗 = 1, 2, 3, are integrals with complicated expressions. In the Appendix, we have listed these integrals from Eqs. (6.1.35)–(6.1.37) in [5]. At small t, the yield surface is close to 𝑦 = 1. Since I2 (y, t) and I3 (y, t) approach zero as 𝑦 → 1− , we can discard them [1]; see also Eqs. (6.1.44) and (6.1.52) in [5] in this connection. The above assumption has also been tested by evaluating 𝐼𝑗 (𝑦, 𝑡), 𝑗 = 1, 2, 3, for various values of Bn and G numerically. These tests have verified that this assumption is reasonable as long as 𝛿(t) is close to unity. 𝑢(𝑦, 𝑡) = 𝑎(𝑡) + 𝑏(𝑡)𝑦 + 𝑐(𝑡)𝑦2 , 𝛿(𝑡) ≤ 𝑦 ≤ 1, 𝑡 ≥ 0, (3.1) where the functions a(t), b(t) and c(t) have to be determined. Here, one appeals to the following boundary conditions: 135 R.R. Huilgol, A.N. Alexandrou and G.C. Georgiou Journal of Non-Newtonian Fluid Mechanics 265 (2019) 133–139 Fig. 5. Evolution of the velocity in the interval [0, t3 ] for Bn = 1 and G = 10 using the parabolic approximation (3.2) (solid lines) and approximation (2.26) for the velocity (dashed lines). Fig. 6. Evolution of the velocity in the interval [0, t3 ] for Bn = 1 and G = 100 using the parabolic approximation (3.2) (solid lines) and approximation (2.26) for the velocity (dashed lines). 1. The velocity at the boundary is zero. Hence, 𝑢(1, 𝑡) = 0 ⇒ 𝑎(𝑡) + 𝑏(𝑡) + 𝑐(𝑡) = 0. 2. The velocity at the yield surface is that of the core. Thus, 𝑢(𝛿(𝑡), 𝑡) = 𝑢𝑐 (𝑡) ⇒ 𝑎(𝑡) + 𝑏(𝑡)𝛿(𝑡) + 𝑐(𝑡)𝛿(𝑡)2 = 𝑢𝑐 (𝑡). 3. The shear rate at the yield surface is zero. Or, 𝜕 𝑢(𝑦, 𝑡)∕𝜕 𝑦 = 0 at 𝑦 = 𝛿(𝑡). Thus 𝑏(𝑡) + 2𝑐(𝑡)𝛿(𝑡) = 0. Using the above conditions, we obtain the velocity field in the fluid as follows: ⎧ ]𝑢𝑐 (𝑡), ⎪[ 𝑢(𝑦, 𝑡) = ⎨ 2 ⎪ 𝑢𝑐 (𝑡)∕(1 − 𝛿(𝑡)) (1 − 𝑦)[1 + 𝑦 − 2𝛿(𝑡)], ⎩ 0 ≤ 𝑦 ≤ 𝛿(𝑡), 𝑡 > 0, 𝛿(𝑡) ≤ 𝑦 ≤ 1, 𝑡 > 0. (3.2) Fig. 4. Comparisons of the core velocity predicted by Eq. (2.18) (solid line) and Eq. (2.27) (dashed) line for Bn = 1: (a) G = 5 (𝛼 = 0.7231); (b) G = 10 (𝛼 = 0.9627); G = 100 (𝛼 = 1.6056). We shall now compare the velocity profile in the yielded region given by Eq. (2.26) with the parabolic approximation in Eq. (3.2). From Eqs. (2.22) and (2.23), we note that 𝛿(𝑡3 ) = 𝛿∞ = Bn∕𝐺. Hence, in Fig. 5, where Bn = 1 and 𝐺 = 10, one has 𝛿(𝑡3 ) = 𝛿∞ = 0.1, with 𝛼 = 0.9627 and 𝑡3 = 0.2185. In Fig. 6, where Bn = 1 and 𝐺 = 100, we see that 136 R.R. Huilgol, A.N. Alexandrou and G.C. Georgiou Journal of Non-Newtonian Fluid Mechanics 265 (2019) 133–139 Fig. 7. Evolution of the velocity in the interval [0, t09 ] for Bn = 1 and G = 10 using the parabolic approximation (3.2) (solid lines) and approximation (2.26) for the velocity (dashed lines). Fig. 9. Evolution of the velocity in the interval [0, t09 ] for Bn = 1 and G = 1000 using the parabolic approximation (3.2) (solid lines) and approximation (2.26) for the velocity (dashed lines). to a constant pressure gradient is now complete. The earlier work of Safronchik [1] had delivered the location of the yield surface and the velocity in the core; for a summary, see [5]. Here, we have found the velocity in the yielded zone. This exact solution approximates quite nicely the velocity in the core for high values of the pressure gradient. A comparison of this exact solution with that obtained by a parabolic approximation shows that the two differ from one another at low values of the pressure drop G and converge as G increases. Finally, solutions to initial boundary value problems in Bingham fluids using the Laplace Transform have appeared in the literature. Using this method, Daprà and Scarpi [7] examined the start-up of channel flow, which is the same as that studied here. Subsequently, they applied the same technique to the start-up flow in a pipe of circular cross-section [8]. More recently, Wu and Liu [9] employed the Laplace transform technique to the start-up flow of a Bingham fluid between coaxial cylinders under a constant wall shear stress. These solutions are incorrect because one cannot use the Laplace transform method for these initial boundary value problems. For a detailed explanation, see Huilgol [10]. Fig. 8. Evolution of the velocity in the interval [0, t09 ] for Bn = 1 and G = 100 using the parabolic approximation (3.2) (solid lines) and approximation (2.26) for the velocity (dashed lines). Acknowledgements This paper is dedicated to the memory of our friend and colleague, Andreas Alexandrou. We would also like to thank the two reviewers for their helpful suggestions which led to several clarifications and improvements of our manuscript. 𝛿(𝑡3 ) = 𝛿∞ = 0.01, with 𝛼 = 1∕6056 and 𝑡3 = 0.09505. In Figs. 5 and 6, we have plotted 10 velocity profiles corresponding to 𝑡𝑖 = 𝑡3 ∕10, 𝑖 = 1, 2, ⋯ , 9, and 𝑡10 = 𝑡3 . It is clear that the comparisons in Figs. 5 and 6 are not accurate, for the velocity approximation in Eq. (2.26) is based on the assumption that 𝛿(t) ≈ 1. Hence, in Figs. 7–9, we have compared the velocity profile in Eq. (2.26) with that given by the parabolic approximation in Eq. (3.2) in the time interval [0, t09 ], where t09 is the time at which 𝛿(𝑡09 ) = 0.9. It turns out that 𝑡09 = 0.0025∕𝛼 2 . In these figures, we have let Bn = 1, and G ∈ {10, 100, 1000}. Again, ten profiles are shown as t increases from 𝑡1 = 0.1𝑡09 to 𝑡10 = 𝑡09 in equal instalments. The corresponding values of 𝛼 are 0.0002698, 0.000009698 and 0.00005601 respectively. It is clear that the parabolic approximation for the velocity field performs well for small Bn/G ratios. Appendix In this Appendix, we shall list the three integrals I1 (y, t), I2 (y, t) and I3 (y, t), which appear in Eq. (2.24). We begin with [5] 2𝐺 𝐼1 (𝑦, 𝑡) = √ 𝜋 ∫0 = −𝐺 ∫0 𝑡 ( ∫∞ ) 2 𝑒−𝛼 𝑑 𝛼 𝑑 𝜎 𝑡 [1 − erf (𝛽(𝑦, 𝜎))] 𝑑𝜎 = −𝐺𝑡 + 𝐺 4. Concluding remarks 𝛽 ∫0 𝑡 erf (𝛽(𝑦, 𝜎))] 𝑑𝜎, where For a time interval of short duration, the determination of the velocity profile in the start-up flow of a Bingham fluid in a channel due 1−𝑦 𝛽 = 𝛽(𝑦, 𝜎) = √ , 2 𝑡−𝜎 137 𝑦 < 1, 𝑡 ≥ 0. (A.1) (A.2) R.R. Huilgol, A.N. Alexandrou and G.C. Georgiou Journal of Non-Newtonian Fluid Mechanics 265 (2019) 133–139 We now turn to the other two integrals in Eq. (2.24). First of all, I2 is given by (Huilgol [5]) ) ( 𝑧1 (𝑦,𝜎) 𝑡 2 1 (A.16) 𝜙′ (𝜎) 𝑒−𝛽 𝑑 𝛽 𝑑 𝜎, 𝐼2 (𝑦, 𝑡) = − √ ∫𝑧1 (𝑦,𝑡) 𝜋 ∫0 We can now evaluate the integral on the right side in Eq. (A.1) as follows. First of all, the error function erf(x) is defined through 𝑥 2 2 𝑒−𝑠 𝑑𝑠. erf(𝑥) = √ 𝜋 ∫0 (A.3) Hence, using integration by parts, we obtain 2 2 1 2 𝑥 𝑒−𝑥 𝑑𝑥 = 𝑥 erf(𝑥) + √ 𝑒−𝑥 . erf(𝑥) 𝑑𝑥 = 𝑥 erf(𝑥) − √ ∫ 𝜋∫ 𝜋 where 𝑦 − 𝛿(𝜎) 𝑧1 (𝑦, 𝜎) = √ , 2 𝑡−𝜎 (A.4) Next, and the lower limits are given by 1 erf(𝑥) 𝑑𝑥, erf(𝑎𝑧) 𝑑𝑧 = ∫ 𝑎∫ ⎧ ∞, ⎪ 𝑧1 (𝑦, 𝑡) = ⎨ 0, ⎪−∞, ⎩ (A.5) for any constant a. In Eq. (A.5), let 1−𝑦 , 2 𝑎= 𝑧= √ 1 (A.6) . 𝑑𝑧 1 = 𝑧3 . 𝑑𝜎 2 Hence, ∫ erf(𝑎𝑧)𝑧−3 𝑑𝑧. √ √ 𝛼( 𝑡 − 𝜎) ≤ 0. 𝑧1 (𝛿(𝑡), 𝜎) = − √ 𝑡−𝜎 (A.8) ∞ 𝑡 erf(𝑎𝑧)𝑧−3 𝑑𝑧 ∫1∕√𝑡 [ ] ∞ √ erf(𝑎𝑧) 2𝑎 1 −𝑎2 𝑧2 = − lim 𝑑𝑧. 𝑒 +𝑡 erf(𝑎∕ 𝑡)+ √ √ 2 𝑧→∞ 𝑧 𝜋 ∫1∕ 𝑡 𝑧2 erf(𝛽(𝑦, 𝜎)) 𝑑𝜎 = 2 ∫0 Next, it is easy to see that ] [ 2 2 2 2 𝑑 1 −𝑎2 𝑧2 1 𝑒 = − 𝑒−𝑎 𝑧 − 2𝑎2 𝑒−𝑎 𝑧 . 𝑑𝑧 𝑧 𝑧2 𝜙′ (𝜎) = − ∞ 1 −𝑎2 𝑧2 2𝑎 −𝑎2 𝑧2 𝑑𝑧 𝑒 𝑑𝑧 = − √ √ √ √ 𝑒 𝜋 ∫1∕ 𝑡 𝑧2 𝜋 ∫1∕ 𝑡 [ −𝑎2 𝑧2 ] 𝑒 2𝑎 √ −𝑎2 ∕𝑡 2𝑎 +√ 𝑡𝑒 . − √ lim 𝑧 →∞ 𝑧 𝜋 𝜋 Since 𝑒−𝑎 2 𝑧2 ∞ ⎧ ∞, ⎪ 𝑧2 (𝑦, 𝑡) = ⎨ 0, ⎪−∞, ⎩ (A.12) Now, (A.13) 𝑧2 (𝛿(𝑡), 𝜎) = Next, 𝐺 ∫0 (A.24) 𝑦 > 𝛿(𝑡), 𝑦 = 𝛿(𝑡), 𝑦 < 𝛿(𝑡). √ √ 𝛼( 𝑡 + 𝜎) ≥ 1. √ 𝑡−𝜎 (A.25) (A.26) Thus, one can derive ( 𝑧2 (𝛿(𝑡),𝜎) ) 𝑡 2 1 𝜙′ (𝜎) 𝑒−𝛽 𝑑 𝛽 𝑑 𝜎, 𝐼3 (𝛿(𝑡), 𝑡) = √ ∫0 𝜋 ∫0 (A.14) Consequently, 𝑡 2 − 𝑦 − 𝛿(𝜎) , √ 2 𝑡−𝜎 ∞ [ ∞ ∞ √ ] 2 2 2 2 1 2 𝑒−𝑎 𝑧 𝑑 𝑧 = √ 𝑒−𝑥 𝑑 𝑥 = 1 − erf(𝑎∕ 𝑡) . √ √ √ 𝑎 𝜋 ∫1∕ 𝑡 𝑎 𝜋 ∫𝑎 ∕ 𝑡 (A.22) and → 0 as z → ∞, we obtain 2𝑎 2𝑎 √ −𝑎2 ∕𝑡 4𝑎 3 1 −𝑎2 𝑧2 −𝑎2 𝑧2 𝑡𝑒 −√ 𝑒 𝑑𝑧 = √ 𝑑 𝑧. √ √ √ 𝑒 2 𝜋 ∫1∕ 𝑡 𝑧 𝜋 𝜋 ∫1∕ 𝑡 Bn Bn =− √ . 𝛿(𝜎) 1 − 2𝛼 𝜎 where (A.11) 𝑧2 (𝑦, 𝜎) = 4𝑎3 (A.21) Unlike I1 (y, t), we have been unable to evaluate I2 (y, t) analytically. Next, ) ( 𝑧2 (𝑦,𝜎) 𝑡 2 1 𝐼3 (𝑦, 𝑡) = √ 𝑒−𝛽 𝑑 𝛽 𝑑 𝜎, (A.23) 𝜙′ (𝜎) ∫𝑧2 (𝑦,𝑡) 𝜋 ∫0 (A.10) Thus, ∞ 1 𝜙′ (𝜎) erf (−𝑧1 (𝛿(𝑡), 𝜎)) 𝑑𝜎, 2 ∫0 where [5] Since erf(az) is bounded as z → ∞, the first term goes to zero in the limit. Hence, ∞ √ 1 −𝑎2 𝑧2 2𝑎 𝑒 𝑑 𝑧. erf(𝛽(𝑦, 𝜎)) 𝑑 𝜎 = 𝑡 erf(𝑎∕ 𝑡) + √ √ 𝜋 ∫1∕ 𝑡 𝑧2 (A.20) 𝑡 = (A.9) 𝑡 (A.19) Hence, noting that when x < 0, one has erf (𝑥) = − erf (−𝑥), we obtain ) ( 𝑧1 (𝛿(𝑡),𝜎) 𝑡 2 1 𝐼2 (𝛿(𝑡), 𝑡) = − √ 𝑒−𝛽 𝑑 𝛽 𝑑 𝜎, 𝜙′ (𝜎) ∫0 𝜋 ∫0 Here, we appeal to the Tables of Integrals due to Ng and Geller [6]. Using Eq. (14) in Section 4.1, we obtain ∫0 (A.18) Thus, (A.7) erf(𝛽(𝑦, 𝜎)) 𝑑𝜎 = 2 𝑦 > 𝛿(𝑡), 𝑦 = 𝛿(𝑡), 𝑦 < 𝛿(𝑡). We recall from Eq. (2.15) that √ 𝛿(𝑡) = 1 − 2𝛼 𝑡. 𝑡−𝜎 Thus, ∫ (A.17) ] ∞ √ 2𝑎 √ −𝑎2 ∕𝑡 4𝑎 3 −𝑎 2 𝑧 2 𝑒 𝑑𝑧 𝑡𝑒 −√ erf(𝛽(𝑦, 𝜎)) 𝑑𝜎 = 𝐺 𝑡 erf(𝑎∕ 𝑡)+ √ √ 𝜋 𝜋 ∫1∕ 𝑡 [ ] √ (1 − 𝑦) √ [−(1−𝑦)2 ∕4𝑡] = 𝐺 𝑡 erf[(1 − 𝑦)∕2 𝑡] + √ 𝑡𝑒 𝜋 2 √ (1−𝑦) −𝐺 (1−erf[(1−𝑦)∕2 𝑡]), 𝛿(𝑡)≤𝑦≤1, 𝑡>0. 2 (A.15) [ 𝑡 = 1 𝜙′ (𝜎) erf (𝑧2 (𝛿(𝑡), 𝜎)) 𝑑𝜎. 2 ∫0 (A.27) Once again, we have been unable to evaluate I3 (y, t) analytically. Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.jnnfm.2018.10.009. 138 R.R. Huilgol, A.N. Alexandrou and G.C. Georgiou Journal of Non-Newtonian Fluid Mechanics 265 (2019) 133–139 References [6] E.W. Ng, M. Geller, A table of integrals of the error function., J. Res. Natl. Bur. Standards B. Math. Sci. 73B (1969) 1–20. [7] I. Daprà, G. Scarpi, Start-up of channel-flow of a bingham fluid initially at rest., Rend. Mat. Acc. Lincei, Ser. 9 (15) (2004) 125–135. [8] I. Daprà, G. Scarpi, Start-up flow of a bingham fluid in a pipe., Meccanica 40 (2005) 49–63. [9] Y.-H. Wu, K.-F. Liu, Start-up flow of a bingham fluid between coaxial cylinders under a constant wall shear stress., J. Non-Newt. Fluid Mech. 223 (2015) 116–121. [10] R.R. Huilgol, Comments on the paper by Y-H. WU and K-F Liu : start-up flow of a Bingham fluid between coaxial cylinders under a constant wall shear stress., J. Non-Newt. Fluid Mech. 223 (2015) 116–121. J. Non-Newt. Fluid Mech. 244 (2017) 116–121. [1] A.I. Safronchik, Non-steady flow of a visco-plastic material between parallel walls, J. Appl. Math. Mech. (PMM) 23 (1959) 1314–1327. [2] A.I. Safronchik, Rotation of a cylinder with a variable angular velocity in a visco– plastic medium., J. Appl. Math. Mech. (PMM) 23 (1959) 1504–1511. [3] A.I. Safronchik, Unsteady flow of visco-plastic material in a circular tube., J. Appl. Math. Mech. (PMM) 24 (1960) 200–207. [4] K. Sekimoto, An exact non-stationary solution of simple shear flow in a bingham fluid, J. Non-Newt. Fluid Mech. 39 (1991) 107–113. [5] R.R. Huilgol, Fluid mechanics of viscoplasticity, Springer Verlag, Berlin Heidelberg, 2015. 139