A Graphical Technique for Solving the Couette-Poiseuille Problem for Generalized Newtonian Fluids

This paper addresses the mixed Couette-Poiseuille problem, that is the flow between two parallel plates, in the presence of simultaneous pressure gradient and wall motion. Instead of the wall-normal coordinate y, we use the local shear stress as our primary variable and rewrite the corresponding formulae for the velocity profile, flow rate, etc. This gives rise to a graphical technique for solving the problem in the case of arbitrary (possibly measured) generalized Newtonian fluid rheology. We demonstrate the use of the proposed technique on two problems: (a) Bingham fluid and (b) a non-Newtonian fluid with general, nonmonotonous viscosity function.


Introduction
The general solution of the Couette-Poiseuille problem (CPP), i.e. the flow between two parallel plates driven by pressure gradient and wall motion simultaneously is not only a fundamental problem of analytical fluid mechanics but it is also relevant in numerous practical applications, for example tape casting [1][2][3][4][5][6][7][8][9][10], damper modelling [11][12][13][14][15][16] or drilling [17][18][19][20][21][22][23].If the rheological properties of the fluid are known, the flow properties such as e.g.velocity profile or volume flow rate can be calculated by evaluating two integrals and determining the integration constant numerically from boundary conditions.Following this method, several expressions were derived and published for different rheological models, e.g.Newtonian [24], Bingham [7], power law [25] etc.These semi-analytical-methods are applicable in the case of simpler non-Newtonian rheology, however, up to the authors' best knowledge, there exists no general approach for coping with more complex fluids, e.g.non-monotonous shear-thickening fluids [26][27][28][29][30]. Therefore, our goal in this paper is to give a general approach to calculate the velocity profile, volume flow rate in the case of generalized Newtonian fluids.Moreover, the proposed graphical technique does not need elaborate mathematical tools.
Non-Newtonian fluids are present in a vast number of engineering areas.In ceramic tape casting engineers must calculate the casting thickness precisely to design the proper devices and process, see [31].To accomplish this aim, it is necessary to model the flow between the moving conveyor and the so-called doctor blade.These materials can exhibit different rheological behaviour according to [31] e.g.Newtonian, power law, Bingham, Herschel-Bulkley, Casson, Prandtl-Eyring or Yasuda rheology, but only the first 4 were studied, see [1][2][3][4][5][6][7][8][9][10].
In vibration isolation design an adequate damper should be developed, see [32] for details.The basic mechanisms of these devices can be explained simply: due to external forces pressure builds up in inner chambers and fluid starts to flow through gaps providing the required damping.In the past these fluids were mostly Newtonian (water, oil), however, with the increasing vibration isolation demand, magnetorheological (MR) fluids became more significant.These fluids have viscoplastic properties and are modelled by Bingham, Herschel-Bulkley or Casson fluids [11][12][13][14][15][16].However, some MR fluids exhibit non-monotonous behaviour, which can be described only by complex models as in [30].
The importance of the CPP is especially high in wellbore hydraulics when trying to determine the surge and swab pressures.The most common geometry is an annulus with moving inner wall [17][18][19][20][21][22][23].This type of flow is investigated by [23] for Robertson-Stiff, [17] for power law and in [45] for Bingham fluid.The parallel plate assumption is valid in case of small clearances (compared to the radius), see [46], which is then used as a submodule in more complicated models as in [18,[20][21][22].
For pure Poiseuille problem Rabinowitsch [42] derived a general expression for the velocity profile and volume flow rate for arbitrary rheology.Furthermore, an expression was also given for the inverse problem, that is to determine the rheological curve from volume flow rate-pressure gradient data.This method is often referred to as differentiation method [47] and is used to evaluate capillary viscometer measurement data [43].This idea is applied to plane Poiseuille flow by [44] and is used commercially in e.g.chip rheology measurement [48].
Despite the existence of a general solution for Poiseuille flow there is no general derivation of flow behaviour in the case of generalized Couette (mixed Poiseuille-Couette, CPP) flow.In the case of Newtonian flow the pure Poiseuille and pure Couette flow can be superposed, therefore it can be solved fully analytically.However, in the case of non-Newtonian fluid, the describing equations are non-linear, thus the solution cannot be superimposed.The full analytical solution exists only for specific cases but mostly even these are semi-analytic solutions, meaning that the differential equations are solved analytically, but the integration constants (and hence the final solution) could be determined only numerically from the boundary conditions.In more complex cases, even the integration itself has to be performed numerically.
To overcome this massive need of mathematical machinery in the case of general fluids, the present study introduces a practical method, which can be used in engineering practice without any limitation on the rheological model.Furthermore, if the resolution of viscosity data is sufficiently fine the velocity profile and volume flow rate can be reconstructed directly from viscosity data.
The rest of the paper is organized as follows.Section 2 describes the mathematical toolbox used in which the novelty is that we use the shear stress as our main variable and build all the equations on it.Section 3 gives an example for simplified, purely Poiseuille flow, while Section 4 gives the main result, which is a technique based on graphical approach.Section 5 gives some examples for general fluids while Section 6 concludes the work.

Generalized Couette flow for general fluid 2.1 Governing equations
We consider the flow of an incompressible liquid in a channel driven by a pressure gradient and a moving wall.The governing equations are the continuity and the momentum equation in Cartesian coordinate system (x, y, z).To obtain an analytical solution the following assumptions are made: 1) the fluid is incompressible 2) steady state solution, no time dependence, 3) 2D flow, i.e. all properties are constant in the z direction, 4) the velocity component in the x (streamwise) direction depends only on the y (spanwise) coordinate (u = u (y)), the other components are zero v = w = 0 and 4) the pressure gradient is constant, i.e. dp dx const

=
. Based on these assumptions the continuity and momentum equation in the y and z-directions are automatically fulfilled.The remaining simplified momentum equation is where τ (y) is the local shear stress.Theoretically, the above equation, together with the rheology curve τ = g (γ) can be used to obtain the velocity profile and calculate the volume flow rate for a prescribed pressure gradient.By integrating Eq. ( 1) from y = 0 (axis of the channel) the following expression can be obtained to local shear stress: where τ 0 = τ (0) is the local shear stress at the axis of the slot.
We wish to describe the flow of fluids with general rheology, hence we introduce the f (τ) shear rate-shear stress function, which is the inverse function of the above defined g  γ ( ) . With the help of this inverse rheology curve, computing the velocity profile is equivalent of solving the following boundary value problem (BVP): where h is the half width of the slot and u w is the velocity of the moving wall.The sign function is used to cope with negative shear stresses, as the rheology curve is defined only for positive shear stress values.The sign f τ τ ( ) ( ) function is the f odd (τ) odd extension of f (τ) inverse rheological function.

Analytical solution
The general solution for the velocity profile can be obtained easily from Eq. ( 3) by integrating with respect to y.It can be observed, that in case of  4) is Equation ( 5) is a well-known expression that can be found in handbooks (for example, [24]) and it is also clear that the velocity profile does not depend on the rheological properties of the fluid, only the required force to maintain u w wall velocity depends on these properties (which is the working principle of parallel plate viscometer).The volume flow rate per unit width can be obtained after the definite integral of the velocity profile from −h to h with respect to y where W is the width of the channel.If pressure gradient is applied on the slot, the general solution for velocity profile has the form of after integration of Eq. ( 3) respecting y and rewriting.2), where τ p is the wall shear for the pure Poiseuille case In this case Eq. ( 7) becomes where C CP is the integration constant for CPP problem and F odd is the integral of f odd , that is Typical f odd (τ) and F odd (τ) functions of a shear thinning power law fluid can be seen on Fig. 2. The particular solution considering boundary conditions ( 4) is Fig. 1 Qualitative sketch of the possible velocity profiles depending on velocity of moving wall and applied pressure gradient Equation (10) contains the yet unknown τ 0 parameter which is chosen to satisfy the boundary condition (moving wall), . ( Equation ( 12) is a nonlinear algebraic equation which can be solved for a given τ P , u w , and h for τ 0 by e.g.Newton-Raphson iteration.Depending on the actual value of τ 0 , there are several cases that might arise, see Fig. 1.
The y location at which τ (y) vanishes shows the local maximum of velocity profile.The case when the shear is positive or negative in the whole slot (−h ≤ y ≤ h) will be referred to as Couette dominated flow (Fig. 1 a) and Couette dominated backflow (Fig. 1 e).The sufficient condition of this case is τ 0 > τ P and τ 0 < −τ P , respectively.The τ 0 solution of Eq. ( 12) equals to zero only if the righthand side of the equation is zero as well i.e. u w = 0 , which is the pure Poiseuille type (Fig. 1 c).The remaining cases are −τ P < τ 0 < τ P , which means the velocity profile has local maximum in the slot.If u w positive and therefore τ 0 > 0, the velocity profile has no negative value (Fig. 1 d).Nevertheless, if u w is negative and therefore τ 0 < 0 the velocity profile has negative and positive values as well (Fig. 1 b).
The volume flow rate per a unit width (in the z direction) is the definite integral of the velocity profile from −h to h with respect to y.Similarly, to integrate Eq. ( 7), substitution integral can be used.The main benefit of this step is to obtain a relation between the volume flow rate and the pressure gradient, which is especially advantageous in the case of pure Poiseuille flow as Section 3 will describe.The flow rate per unit channel width is where )

Pure Poiseuille flow 3.1 Determination of volume flow rate from rheology
In this section we demonstrate the use of the above derived equations by reproducing the results available in the literature.The velocity profile for pure Couette flow is already presented in Eq. (5).A more interesting case is the pure Poiseuille flow when the velocity at the walls is zero, which means τ 0 = 0 according to Eq. (12).To obtain velocity profile and volume flow rate for a prescribed pressure gradient the appropriate F odd (τ) functions should be known.For Newtonian fluid, we have and for power law fluid and finally, for Bingham fluid As the profile is symmetric, we also have Substituting the F odd (τ) functions corresponding to Newtonian ( 16), power law (17) and Bingham (18) model into Eq.( 11) and Eq. ( 14) the velocity profiles and volume flow rate values are for Newtonian for power law and for Bingham fluid which are identical to the ones presented in handbooks, e.g.[49].

Determination of rheology from volume flow rate measurement data
There are cases, when the flow rate is known and the rheology is to be found.Substituting Eq. ( 19) into Eq.( 14) we obtain a first-order ordinary differential equation for S (τ).This type of differential equation can be written in the form of which has an analytic solution where C is an arbitrary constant and In our case the meanings of notations are ) and the solution is By differentiating twice the right-hand side of Eq. ( 28), the rheological behaviour f (τ P ) is This expression contains only the volume flow rate-pressure gradient relationship which can be measured by e.g.parallel plate viscometer (rheometer-on-a-chip [48]).The shear rate-volume flow rate relationship was first derived for pipe in [42] and then in parallel plate slot.The above presented relationship is exactly the same as published one [44].

Graphical representation of generalized Couette flow
Let us assume in this section, that the f (τ) inverse rheological curve of an arbitrary fluid is given by the measured points that are sufficiently close to each other.If one wanted to determine the volume flow rate for a prescribed u w wall velocity and pressure gradient directly, the following steps are needed: 1. choose and fit a rheological model to the measurement data, 2. substitute it back to (1) and solve the corresponding boundary value problem (BVP).However, this can be challenging in the case of complex fluids, and can be performed only numerically.Moreover, this method requires a BVP solver code, which (although available in most scientific programming packages) might be cumbersome to use.
However, the velocity profile and volume flow rate values can be obtained without any programming if the above presented steps are made graphically.Based on Eq. ( 11) and Eq. ( 14) it is clear, that only F odd (τ) and S odd (τ) functions are needed.If the measured data is sufficiently fine, F odd (τ) function can be prepared by numerical integration of measured data, e.g. by trapezoid method.Then S odd (τ) function can be also plotted by another numerical integration.To obtain τ 0 value we need the root of which is also easy to locate graphically: we are searching for a "window" of width 2τ P and height u h w P τ as shown in Fig. 2 b, where τ P is defined by (8).Let F A = F odd (τ 0 − τ P ) and where The velocity profile can be plotted using where F y corresponds to shear stress value at y If the wall velocity is negligible, then there is no need to read root of Eq. ( 30) from graph, thus the volume flow rate-pressure gradient connection can be depicted directly.The presented analysis can be easily performed with a spreadsheet software as demonstrated in the next section.

Application examples 5.1 Bingham fluid
To verify the validity of the graphical method, the tape casting thickness prediction presented in [8] will be repeated here by the analytic and graphical approach.The authors of [8] determined the theoretical and experimental values for 24 parameter combinations.In our study the 12 th row of Table 1 in [8] is chosen.These parameters are summarized in Table 2.

Analytical solution
To obtain the exact value of δ tp tape thickness [8] (where K = 0.7464 is a constant related to side-flow effect), the τ 0 should be calculated from Eq. ( 12) with the help of Eq. ( 18).Because of the absolute values in Eq. ( 18) it is important to identify which case of Fig. 1 will occur.Let us assume the case "e" of Fig. 1 applies.Then Eq. ( 12) has the form of which has the solution satisfying τ 0 > τ P , (thus the previous assumption of "e" was correct).The volume flow rate per unit width from Eq. ( 14) is from which the tape thickness is δ tp = 1.3536 [mm], see also [8].

Graphical solution
The main advantage of graphical method is that the volume flow rate values can be determined from read values from diagrams, which are drawn from measured data, thus model fitting and solving boundary value problem is not needed.Let us assume the f (τ) measurement data as shown in Fig. 3.a.The first step is to plot F (τ) (Fig. 3 b) and S (τ) (Fig. 3 c).The τ 0 value can be read if the 2τ

Fluid with nonmonotonous viscosity
Shear thickening effect of non-Newtonian fluids are advantageous in vibration isolation or impact protection application.These fluids are usually suspensions i.e. solid particles are solved in carrier fluids.However, these types of fluids can exhibit non-monotonous viscosity behaviour, i.e. in low range of shear rate slightly shear thinning can be observed, but after critical shear rate the viscosity increases rapidly, as Fig. 4 shows.An easily available example is the cornstarch suspension [27], however this phenomenon also appears in magnetorheological fluids [30].
There is no general model to investigate nonmonotonous viscosity fluids (NMF).David and Filip had published few papers about possible empirical modelling facilities.In [50], a 6-parameter model is presented which can be used in situations, only when one local extreme appears.Galingo-Rosales [26] divides the viscosity curve into 3 regions and describes the fluid by models with different parameters in these regions, so an 11-parameter model is employed to describe the rheology.In [29] the authors decrease the number of parameters to 10, but their model is cumbersome to handle analytically.
Because of the fact that typical application of NMF is vibration isolation and shock absorbing, it can be necessary Fig. 3 The f odd (τ), F odd (τ), S odd (τ) functions corresponding to Bingham fluid [8], and the applied values in to cope with CPP problem with NMF.As we see, adequate modelling of NMF is complicated and these complex rheologies make the analytical solution presented in Section 2.2 cumbersome to derive.Another way is the previously introduced graphical approach, where only the viscosity measurement data is needed and no model fitting is required.Let us determine volume flow rate and velocity profile in case of viscosity measurement data of [28] (Fig. 4) for Then the appropriate f odd (τ) directly, F odd (τ) S odd (τ) after numerical integration can be plotted (Fig. 5).
From the graph the following data can be read From the points of Fig. 5, the velocity profile can also be sketched with the help of Eq. ( 34) and (35) to illustrate the flow.This profile can be seen on Fig. 6, where we also depicted (solid line) the numerical solution of the boundary value problem defined by Eq. ( 3) and (4).

Conclusion
A graphical technique was presented to solve the mixed Couette-Poiseuille problem, that is the flow between two parallel plates in the simultaneous presence of a pressure gradient and wall motion.We have reformulated the  governing equation in terms of the shear stress (instead of the traditional wall-normal coordinate y) that allowed us rewrite the flow rate, velocity profile etc. purely using the rheological curve.After two numerical integrations (e.g. by trapezoidal rule) of this (possibly measured) rheological curve, one can read the quantities required to compute e.g. the flow rate.The main advantage of the proposed technique is that there is no need for elaborate analytical or programming skills.The proposed technique was demonstrated on two problems (Bingham fluid and on a fluid with nonmonotonous viscosity) and its viability was proven.
dp dx = 0 i.e. pure Couette flow the shear stress is constant across the slot, which means that d dy τ = 0 and the general solution for velocity profile is a linear function u y y C C ( ) = + τ 0 for any rheology, where C C is the integration constant for pure Couette flow.The particular solution of pure Couette case considering boundary conditions (

Fig. 2
Fig. 2 Method of determining volume flow rate in case of a prescribed pressure gradient and wall velocity.If the measured inverse rheological points (a) are closed, the F odd (τ) (b) and F d odd τ τ ( ) ∫ and vice versa), hence only one solution exists, despite the symmetry of F odd (τ) function.Then one should read the values of F A , F B , S A , S B and calculate the volume flow rate using notations of Fig.

×
window is adjusted on plotted F (τ) curve (Fig. 3 b), giving Then the F B , S A and S B values turn out to be tape thickness is δ tp = 1.3544 [mm].As we can see the difference between analytical and graphical solution comes only from the inaccuracy of figure reading.

Table 1
Summary of studies on the generalized Couette flow and its application in engineering problems.