Lubricated flow of a viscous liquid in a pipe

This notebook has been written in  Mathematica 4.0  by 

Mark J. McCready
Professor and Chair of Chemical Engineering
University of Notre Dame
Notre Dame, IN 46556
USA

mjm@nd.edu
http://www.nd.edu/~mjm/


It is copyrighted to the extent allowed by which ever laws pertain to the World Wide Web and the Internet.

I would hope that as a professional courtesy, this notice remain visible to other users. 
There is no charge for copying and dissemination   

Version: 10/15/99
More recent versions of this notebook may be available at:
http :  // 
                www . nd . edu/~mjm/lubricatedflow . nb

This notebook uses the notation of and is intended as a companion to:
 An Introduction to Fluid Dynamics by Stanley Middleman, John Wiley, 1999.

Keywords

lubricated transport, two-layer flows, multifluid flow, two-phase flow, laminar flow

Motivation

We have found that pressure drop is proportional to viscosity and diameter to the fourth power for laminar flows.  Thus as the viscosity of a liquid gets much larger than water, it will require either very high pressure drop (high operating costs) or large pipes (high capital costs) to pump.  Obviously, if the liquid is valuable enough then it will be worth it and pipeline built.  However for a less valuable material, pumping as a single phase liquid may not be economically viable.  Such is the case for very heavy crudes.  Pumping costs can be so large that the oil could never be sold for profit.

However, what if your company or country has access to only this kind of heavy crude oil.  It is in your interests to solve this problem and produce a new technology that transport the crude efficiently enough so that it can be sold for profit.

The idea that could be used is "lubricated" transport where the viscous oil is confined to a "core" region in the center of the pipe and a much less viscous liquid (e.g. cheap also, water) flows in the annular region.  The pressure drop is reduced because while the water is sheared highly, the viscous oil core is not subject to much deformation thus does not contribute to pressure drop.

For more about the real problem you will need to read elsewhere.  Real core-annular flows have waves on the interface which usually act to stabilize the core flow, although could sometimes break up the core. 

A general reference for this topic is Fundamentals of two-fluid dynamics by D. D. Joseph and Y. Y. Renardy, Springer Verlag, 1993.  

There is also a CNN science/technology video on this topic "staring" Dan Joseph and his lab.

Preview of major points

This notebook solves the "lubricated" flow in a pipe problem and in doing so demonstrates several major fundamental points in Transport Phenomena.  

1.  For flows where more than one fluid is flowing together, there is a separate equation for the velocity profile for each fluid.  These equations differ just by the fluid properties.  

2.  To solve these two-fluid problems you need to recognize:
	•  The pressure drop is the same in each phase if gravity can be ignored.
	•  The velocity of the two fluids matches at the interface.
	•  The shear stress in each fluid matches at the interface.  
The last two of these are common boundary conditions.  

For this specific situation the major results are:

Lubricated transport can cause a flow rate enhancement of the viscous fluid (compared to this fluid flowing alone in the pipe) on the order of the viscosity ratio of two fluids and the maximum enhancement occurs when the less viscous fluid occupies the outer 30% of the radius.

Flow geometry

The problem we are considering is the flow of a very viscous liquid (Fluid I) confined to the center of the pipe, with a much less viscous liquid (Fluid II) flowing around Fluid I and contacting the walls.  The pipe diameter is D and the thickness of the lubricating layer is h.

Solution of the example

Basic equations

The momentum equations for steady, unidirectional flow

We have derived the equation for steady laminar flow in a tube.   When we did this, we did not make restrictions on the type of fluid or even if there is more than one fluid.  The equation is

eq1 = (1/r) D[r (Subscript[τ, zr])[r], r] -  
        D[P[z], z]  +  ρ g

The Newtonian stress constitutive equation

From the previous example in Middleman (3.2.1, page 68-74)  you may recall two points.  First, that P[z] was a linear function and thus P'[z] = const.  We will call this constant, dpdz.  Second, we need a constitutive equation for the stress component, (Subscript[τ, zr]), in terms of the velocity gradient.  Since we cannot easily derive this, we just stated that it is: 
		\(TraditionalForm\`τ\_zr = \(-μ\)\ ∂v\_z\/∂r\).
		
Now we need to recognize that there are two fluids and thus two values of viscosity for the fluids in different regions.  Thus we will have two versions of eq1, one for fluid I and one for fluid II.  Since we know the constitutive equations and we know that P'[z] is a constant, we might as well substitute and start to solve it.  We also make the restriction that gravity is not important and thus ignore its possible effects by making g=0. 

The z momentum equations are:

Back to preview

Back to conclusions

fluid I equation

eq11 = 
      eq1 /. {Derivative[i_][(Subscript[τ, zr])][r] :>  
            D[-(Subscript[μ, 1]) D[(Subscript[v, z]^I)[r], r], {r, 
                i}], 
(Subscript[τ, zr])[r] -> -(Subscript[μ, 1]) 
              D[(Subscript[v, z]^I)[r], r], P'[z] -> dpdz1, 
          g -> 0}

fluid II equation

eq12 = 
      eq1 /. {Derivative[i_][(Subscript[τ, zr])][r] :>  
            D[-(Subscript[μ, 2]) D[(Subscript[v, z]^II)[r], r], {r, 
                i}], 
(Subscript[τ, zr])[r] -> -(Subscript[μ, 2]) 
              D[(Subscript[v, z]^II)[r], r], P'[z] -> dpdz2, 
          g -> 0}

A note on Mathematica

The"/." command is a substitution command, which is used frequently when doing computer algebra.  On the left side of /. is the term that is receiving the substitution.  On the right side of the /. is the designation for what we are replacing, either any derivative of (Subscript[τ, zr])or just (Subscript[τ, zr])(we can't do the 0 derivative substitution with this command), with the velocity gradient substituted in.

Solution of the governing equations

General solutions to these differential equations

Let's solve what we can.

Start with the middle of the fluid.  We can first get a general solution to this equation.

Equation 1 solution

eq1soln = DSolve[{eq11 == 0}, (Subscript[v, z]^I)[r], r]

Back to conclusions

It is convenient to examine this solution to see if an obvious simplification arises.  It does.  At r=0, we don't anticipate any reason for the velocity to be infinite.  Since there is no other term to cancel log(r), we must make (Subscript[c, 1])=0.  We can't do anything with the other constant at this moment.  Recall however that in the single phase example, (Subscript[c, 2]) was adjusted to make the velocity = 0 at the fixed wall.   

This real ugly command sets (Subscript[c, 1])=0, and changes the name of (Subscript[c, 2])to f1 (to sort of resemble Middleman's notation.).  It also picks out just the equation without the braces.

eq1a  =  
      (Subscript[v, z]^I)[
          r] /. (eq1soln /. {C[1] -> 0, 
                C[2] -> f1})[[1]]

Equation 2 solution

We can solve the second equation the same way.  This time the domain of r does not go to 0 and there is no need (or reason) to exclude the log(r) term -- so it stays.

eq2soln = DSolve[{eq12 == 0}, (Subscript[v, z]^II)[r], r]

Change the name of the constants to e2 and f2 and extract the equation as "eq2a"

eq2a  =  
      (Subscript[v, z]^II)[r] /. (eq2soln /. {C[1] -> e2, 
                C[2] -> f2})[[1]]

We can take stock at this point.  We don't know the three constants, f1, e2, f2 and we expect to need three boundary conditions.  We also worry about and will need to resolve if dpdx1 = dpdx2.  (They do, but we should discuss this -- see page 80).

Boundary conditions

Back to conclusions

Velocity match at the interface

Back to preview

Just as fluids stick to walls, they also stick to each other so we require that at the interface, r =  D/2 - h, the velocities in the two fluids be the same.

vmatcheq = (eq1a /. 
              r -> (D/2 - 
                    h))  -  
                  (eq2a /. r -> (D/2 - h)) == 0

Back to conclusions

No slip at the wall

The outer fluid, phase II, is touching the pipe wall and thus the velocity there is 0.

noslipeq = (eq2a /. r -> (D/2)) == 0

Back to conclusions

Stress match at the interface

Back to preview

We have used up conditions on the velocity at locations that we can think of (center of pipe, wall, interface) thus we think of possible boundary conditions on the first derivative.  This, of course, relates to the shear stress in the fluid and we have condition that arises from the equation that was derived originally.   The original equation (eq1), tells us that the stress profile is linear across the pipe with a value of 0 in the middle and a largest (in magnitude) value at the wall.  This result is exclusive of the values of viscosity of the fluid or if there are different fluids in the pipe.  

Thus a condition on the first derivative is that the shear stress is continuous across the pipe.  This means that the shear stress in each fluid is equal at the interface.

stresseq = ((Subscript[μ, 1])  D[eq1a, r] /. 
              r -> (D/2 - 
                    h))  -       
                     ((Subscript[μ, 2]) D[eq2a, r] /. 
              r -> (D/2 - h)) == 0

Back to conclusions

Solve for the constants using the three boundary conditions.

We can solve for the three constants using the 3 boundary condition equations.  The answer is in the form of a substitution rule telling what, f1, f2 and e2 should be replaced by.

bceqs1 =  {vmatcheq, noslipeq, stresseq}

bceqs2 =  Simplify[bceqs1]

bcans = Solve[bceqs2, {e2, f2,  f1}]

Pressure drop is same in both phases

Back to preview

We can make an argument that the pressure across a cross-section must be constant, (or just increase linearly downward because of the gravitational hydrostatic force) otherwise there will be flow in the radial direction to relax the pressure.  Since we are requiring only one dimensional flow, there can be no pressure change in the r direction.  Now that the pressure is constant across a cross section, we realize that the pressure change in the axial (z) direction must be the same in each phase.

Velocity profiles

Let's get the answers for the two velocity profiles by substituting in the answers from the boundary conditions.  Note that "[[1]]" takes the expression that is inside one set of braces, { $$$ }.

v1ans = 
      FullSimplify[(eq1a /. bcans[[1]]) /. {dpdz2 -> dpdz, 
            dpdz1 -> dpdz}]

v2ans = 
      FullSimplify[(eq2a /. bcans[[1]]) /. {dpdz2 -> dpdz, 
            dpdz1 -> dpdz}]

Dimensionless profiles

Following Middleman, we can simplify the equations and get a better feel for how they behave if we nondimensionalize these.  Note how we do substitutions to get these forms.

v1profile = 
      FullSimplify[(v1ans /. {D -> 2 R,  
                r -> s R,  (Subscript[μ, 2]) -> (Subscript[μ, 1])/M, 
                h -> η R})/(dpdz R^2/
                  4/(Subscript[μ, 1]))]

v2profile = 
      FullSimplify[(v2ans /. {D -> 2 R,  
                r -> 
                  s R,  (Subscript[μ, 2]) -> (Subscript[μ, 1])/M})/(dpdz R^2/
                  4/(Subscript[μ, 1]))]

These match equations 3.2.53 and 3.2.54, although you have to re arrange the first one a little.

Investigation of the solution

Plots of the velocity profiles

Let's make some plots to show how it behaves.  First the inside fluid is twice as viscous as the outside fluid. M=2

vp22 = Plot[v2profile /. M -> 2, {s, 1,  .8}];

vp21 = 
        Plot[v1profile /. {M -> 2, η ->  .2}, {s,  .8, 
            0}];

Show[vp22, vp21];

Now M=10

vp102 = Plot[v2profile /. M -> 10, {s, 1,  .8}];

vp101 = 
        Plot[v1profile /. {M -> 10, η ->  .2}, {s,  .8, 
            0}];

Show[vp102, vp101];

Now M=100

vp1002 = Plot[v2profile /. M -> 100, {s, 1,  .9}];

vp1001 = 
        Plot[v1profile /. {M -> 100, η ->  .1}, {s,  .9, 
            0}];

With this one we see very clearly that the center of the pipe is flowing with much less shear that the outside fluid

Show[vp1002, vp1001];

If you were worried, we can check to see if M=1, if the profile looks correct, it does

vp12 = Plot[v2profile /. M -> 1, {s, 1,  .8}];

vp11 = 
        Plot[v1profile /. {M -> 1, η ->  .2}, {s,  .8, 
            0}];

Show[vp11, vp12];

Finally we could envision the reverse case, the outside fluid is more viscous.  We don't want to use this pump our fluid, but we might be interested to see what it looks like. M=.1

vp012 = Plot[v2profile /. M ->  .1, {s, 1,  .8}];

vp011 = 
        Plot[v1profile /. {M ->  .1, η ->  .2}, {s,  .8, 
            0}];

Show[vp011, vp012];

We can make this last case extreme to where the outside fluid is almost a solid pipe.  We would not consider a real flow like this but this kind of limit is often useful to check answers or help our conceptual thinking. M=.001

vp00012 = 
        Plot[v2profile /. M ->  .001, {s, 1,  .8}];

vp00011 = 
        Plot[v1profile /. {M ->  .001, η ->  .2}, {s,  .8, 
            0}];

It looks like a laminar flow in a smaller pipe!

Show[vp00011, vp00012];

Flowrate enhancement by lubrication

Again following Middleman, we would like the enhancement in the inner fluid flowrate for a fixed pipe size and pressure gradient.  We realize that while adding more outside fluid reduces shear, and thus increases flow, it also reduces the area available for flow and thus could ultimately reduce the flow of the inner flow.

Thus a fair measure of the enhancement is the ratio of the volumetric flow of the inner fluid for the lubricated case to its volumetric flow when η=0, that is when there is no outside fluid.

Expression for the flow rate enhancement

To get an expression for the volumetric flow rate, we integrate the velocity of the inner fluid in the region that there flow of the inner flow.

q1  =  
      Simplify[Integrate[
          2 π s v1profile, {s, 0, 1 - η}]]

If we see η =0 we can get the dimensionless flow rate with no outside phase.

q1η0 =  q1 /. η -> 0

Now the enhancement is this ratio:

Ψ  =  FullSimplify[q1 /q1η0]

Which matches equation 3.2.58 and gives the relative enhancement due to lubrication.

Plots of the enhancement ratio.

We would like to make some log-log plots and we need to load the Graphics package.

Needs["Graphics`Graphics`"]

Note for Mathematica 3 users: the Log-Log plotting package may give trouble and not show correct axis numbering.  Sorry, I don't know how to fix this.

qp10 = 
        LogLogPlot[Ψ /. M -> 10, {η,  .01,  .9}, 
          PlotStyle -> 
            Dashing[{ .02,  .02}], 
AxesLabel -> {"\
η", "Ψ  "    Enhancement}];

qp100 = 
        LogLogPlot[Ψ /. 
            M -> 100, {η,  .01,  .9}, 
PlotStyle -> 
            Dashing[{ .04,  .04}]];

qp1000 = 
        LogLogPlot[Ψ /. 
            M -> 
              1000, {η,  .01,  .9}, 
PlotStyle -> 
            Dashing[{ .08,  .08}]];

qp1 = 
        LogLogPlot[Ψ /. 
            M -> 1, {η,  .01,  .9}];

qpx = 
        LogLogPlot[1, {η,  .01,  .9}, 
          PlotStyle -> Dashing[{ .01,  .01}], 
          AxesLabel -> {"η", "Ψ  \
"    Enhancement}];

Back to preview

Back to conclusions

Show[qpx, qp1000, qp100, qp10, qp1, 
        PlotRange -> All];

We see that the enhancement is highest when the outside fluid occupies about 30% of the pipe.  (How much of the area is this??)  Below this value, there is an advantage to increasing the less viscous flow to reduce shear.  Above this point, there is not enough flow area available for the desired fluid and the enhancement drops.

Velocity profiles for different values of η  at constant M to show how the enhancement occurs.

It is now worth taking a look at the velocity profiles along one of these curves. 

We make a series of plots but don't show them until the end

vp2hp1 = 
        Plot[v2profile /. M -> 10, {s, 1,  .9},  
          DisplayFunction -> Identity];

vp2hp2 = 
        Plot[v2profile /. M -> 10, {s, 1,  .8}, 
          PlotStyle -> Dashing[{ .01,  .01}],  
          DisplayFunction -> Identity];

vp2hp3 = 
        Plot[v2profile /. M -> 10, {s, 1,  .7}, 
          PlotStyle -> Dashing[{ .02,  .02}],  
          DisplayFunction -> Identity];

vp2hp4 = 
        Plot[v2profile /. M -> 10, {s, 1,  .6}, 
          PlotStyle -> Dashing[{ .03,  .03}],  
          DisplayFunction -> Identity];

vp2hp5 = 
        Plot[v2profile /. M -> 10, {s, 1,  .5}, 
          PlotStyle -> Dashing[{ .04,  .04}],  
          DisplayFunction -> Identity];

vp2hp6 = 
        Plot[v2profile /. M -> 10, {s, 1,  .4}, 
          PlotStyle -> Dashing[{ .05,  .05}],  
          DisplayFunction -> Identity];

vp1hp1 = 
        Plot[v1profile /. {M -> 10, η ->  .1}, {s,  .9, 0},  
          DisplayFunction -> Identity];

vp1hp2 = 
        Plot[v1profile /. {M -> 10, η ->  .2}, {s,  .8, 0}, 
          PlotStyle -> Dashing[{ .01,  .01}],  
          DisplayFunction -> Identity];

vp1hp3 = 
        Plot[v1profile /. {M -> 10, η ->  .3}, {s,  .7, 0}, 
          PlotStyle -> Dashing[{ .02,  .02}],  
          DisplayFunction -> Identity];

vp1hp4 = 
        Plot[v1profile /. {M -> 10, η ->  .4}, {s,  .6, 0}, 
          PlotStyle -> Dashing[{ .03,  .03}],  
          DisplayFunction -> Identity];

vp1hp5 = 
        Plot[v1profile /. {M -> 10, η ->  .5}, {s,  .5, 0}, 
          PlotStyle -> Dashing[{ .04,  .04}],  
          DisplayFunction -> Identity];

vp1hp6 = 
        Plot[v1profile /. {M -> 10, η ->  .6}, {s,  .4, 0}, 
          PlotStyle -> Dashing[{ .04,  .04}],  
          DisplayFunction -> Identity];

Show[vp1hp1, vp1hp2, vp1hp3, vp1hp4, vp1hp5, 
        vp1hp6, 
     vp2hp1, vp2hp2, vp2hp3, vp2hp4, 
        vp2hp5, vp2hp6, 
PlotRange -> All,  
        DisplayFunction -> $DisplayFunction];

The series of velocity profile plots is revealing.  We see that as the fraction of less viscous fluid is increased, the maximum velocity of the core is increased.  This allows a higher flow rate until finally the area available for flow is reduced to too small a value.

Reverse ratio showing reduction of flow rate.

Maybe it is interesting to see the behavior for the ratio being reversed.  We see it is always less than one.

LogLogPlot[Ψ /. 
          M ->  .01, {η,  .01,  .9}];

Scaling of the enhancement.  How much can you expect?

Finally, it is tempting to try to scale the enhancement.  It looks almost linear in M.  If we divide out M, we find that for large M , the there is no M variation.

Limit[Ψ/M,  M -> Infinity]

Plots of the scaled enhancment

qp10m = 
        LogLogPlot[(Ψ/M) /. 
            M -> 10, {η,  .01,  .9}, 
          PlotStyle -> 
            Dashing[{ .02,  .02}], 
AxesLabel -> {"\
η", "Ψ/M  "    scaling}];

qp100m = 
        LogLogPlot[(Ψ/M) /. 
            M -> 100, {η,  .01,  .9}, 
PlotStyle -> 
            Dashing[{ .04,  .04}]];

qp1000m = 
        LogLogPlot[(Ψ/M) /. 
            M -> 
              1000, {η,  .01,  .9}, 
PlotStyle -> 
            Dashing[{ .08,  .08}]];

qpxm = 
        LogLogPlot[1, {η,  .01,  .9}, 
          PlotStyle -> Dashing[{ .01,  .01}], 
          AxesLabel -> {"η", "Ψ/M  \
"    scaling}];

Back to preview

Back to conclusions

Show[qpxm, qp1000m, qp100m, qp10m, 
        PlotRange -> All];

We see that the scaling is correct and that the maximum enhancement is about 50% of the viscosity ratio.  Now you know enough to tell if this process for transporting oil will be worthwhile!

Where is the flowrate maximum located?

I can't resist one more calculation.  We can find the maximum enhancement readily.  First take the derivative.

max1 = FullSimplify[D[Ψ, η]]

Set it to 0 and find the value of η.  Isn't Mathematica great!!.

max2 = Solve[max1 == 0, η]

We have observed that it does not change much with M.  Lets see why by looking at the solution for large M.

max3 = Limit[max2, M -> Infinity]

The value of the maximum is about 2.9.  You now know even more about this problem.

N[max3]

Does it make sense that increasing one of the flow rates will decrease the pressure drop?

There is one more interesting point to make about this calculation.  We can examine how the pressure drop changes as the flowrate of the less viscous fluid is increased.  If we increase the flowrate of a fluid in a single phase flow the pressure drop will increase.  This will be linearly for a laminar flow.  However, for a two-fluid flow will this always be the case?  We will find that it is not the case.  For a range of flow rates, increasing the less viscous fluid flowrate will decrease the pressure drop.  

Calculate the dimensional flowrate of the inside fluid starting with the dimensional velocity profile.

v1prof = ((4 dpdz (D - h) h (Subscript[μ, 1]) + dpdz ((((D - 2 h))^2) - 4 (r^2)) (Subscript[μ, 2]))/(16 (Subscript[μ, 1]) (Subscript[μ, 2])))

The flowrate is,

volflow1 = 
      Simplify[Integrate[2 π r v1prof, {r, 0, D/2 - h}]]

Now do the same for the outside fluid.

v2prof = ((dpdz ((D^2) - 4 (r^2)))/(16 (Subscript[μ, 2])))

volflow2 = 
      Simplify[Integrate[2 π r v2prof, {r, D/2 - h, D/2}]]

To make the comparison we will have to solve for the pressure drop, dpdz and h.  We will then plot how pressure drop changes with changes in the volumetric flow rate of the outside fluid.  We will use numbers to show this.

Here is what is solved to get the plot.  Note the order of h.

testeqs = {volflow1 - qq1, volflow2 - qq2}

First find the pressure drop for a single phase flow.  We can do this if h=0.

ansx1 = Solve[(testeqs /. h -> 0) == 0, dpdz]

This looks correct and if we pick some numbers in cgs units for the variables we might find,

N[ansx1 /. {D -> 10, qq1 -> 1000, 
          qq2 -> que2, (Subscript[μ, 2]) ->  .01, (Subscript[μ, 1]) -> 
            10}]

So this is the pressure drop with no lubricating flow.  Adding outside fluid will decrease dpdz from this value.  

Now find the pressure drop as a function of qq2.  We need to use a nonlinear solution technique to get the answer since there are at least 6 solutions for different values of h.  You can explore this more if you wish, or just trust me.  Note that h must be positive and less than D.  

(To make this work in Mathematica 4, I needed to pick out root 5 "[[5]]" of the 6 roots.  I also chose the range of the entire calculation "i,1,20" and the scaling of que2, "  3+(i-1)*Exp[i/3.3]" to show the particular region that I wanted.  This type of behavior will exist for almost any "qq1" but you might imagine that it is difficult to show because it involves a six-order equation for h.)

(To make this work in Mathematica 3, you need to pick the first root, [[1]]. )

Here is what the set of two equations look like for one set of numbers:

testeqs /. {D -> 10, qq1 -> 1000, 
        qq2 -> que2, (Subscript[μ, 2]) ->  .01, (Subscript[μ, 1]) -> 
          10}

ExpandAll[%]

Now get the answers for the pressure drop decrease

Note for Mathematica 3 users: You need to change the [[5]] to [[1]] to pick out the first root.

Do[
que2 = 
        3 + (i - 1)*Exp[i/3.3]; 
ansq1[i] = {que2, 
          dpdz /.  NSolve[(testeqs == 0) /. {D -> 10, 
                    qq1 -> 1000, 
                    qq2 -> que2, (Subscript[μ, 2]) ->  .01, (Subscript[μ, 1]) -> 
                      10}, 
{h,  
                  dpdz}] [[5]]}, 
{i, 1, 
        20}]

plottinglist1 = Table[N[ansq1[i]], {i, 1, 20}]

Plot of pressure drop decrease

Back to conclusions

Note for Mathematica 3 users: The log-log plot labeling problem mare arise here.

LogLogListPlot[plottinglist1, 
        AxesLabel -> {"q2, flowrate", "pressure drop"}];

We see that increasing q2, the outside fluid flow rate decreases the pressure drop for quite a range.  This is not what is expected for a single fluid flow.

Just check h so that we are sure it is OK

Note for Mathematica 3 users: This one does not work very well. You need to change the [[5]] to [[1]] to pick out the first root, but some of the values are the wrong root.

Do[
que2 = 
        3 + (i - 1)*Exp[i/3.3]; 
ansq2[i] = {que2, 
          h /. NSolve[(testeqs == 0) /. {D -> 10, 
                    qq1 -> 1000, 
                    qq2 -> que2, (Subscript[μ, 2]) ->  .01, (Subscript[μ, 1]) -> 
                      10}, 
{h,  
                  dpdz}][[5]]}, 
{i, 1, 
        20}]

plottinglist2 = Table[N[ansq2[i]], {i, 1, 20}]

LogLogListPlot[plottinglist2, 
        AxesLabel -> {"q2", "h"}];

The minimum pressure drop occurs when h ≃ 1.5.  This makes sense because h is supposed to be about 30% of D/2 (= 1.5 for this calculation) for the most efficient flow.

Conclusions

1.  There is a separate equation for the velocity profile of each fluid, but only one equation for the stress profile.  

2.  The form of the differential operator (the differential equation) for a circular geometry leads to a logarithmic term in the solution.  This term does not give a physically realizable profile for the middle of the pipe and must be neglected.  (It is kept for fluids confined away from the middle.

3.  These are readily solved and the key is using the boundary conditions to match the constants.
	For a two-layer flow you expect:
	• Shear stress match at the interface
	• Velocity match at the interface
	• No slip at the solid surface

4.  The enhancement is highest when the outside fluid occupies about 30% of the pipe.  Below this value, there is an advantage to increasing the less viscous flow to reduce shear.  Above this point, there is not enough flow area available for the desired fluid and the enhancement drops. 

5.  The maximum enhancement scales as the value of the viscosity ratio and is approximately 50% of this value.  

6.  We can see explicitly that in contrast to single phase flow (and in contrast to intuition) there is a range where increasing the flowrate of the lubricating fluid will decrease the pressure drop.