Vertical Laminar Flow of Bingham Liquid Film 
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :
Since in this laminar flow problem we have split boundary conditions, its resolution is based on the Shooting Method developed in Mathematica by Professor Brian Higgins (Chemical Engineering Department, University of California, Davis, http://www.higgins.ucdavis.edu/chemmath.php  ).

Vertical Flow of Bingham liquid film

g = 9.81; ρ = 950; (Subscript[τ, 0]) = 5; (Subscript[μ, 0]) = 0.15; δ = 
      0.005;

system[Ω_] = {D[(Subscript[τ, xz])[x], {x, 
                1}] ==  ρ g, 
          D[(Subscript[v, z])[x], {x, 1}] ==  
            If[Abs[(Subscript[τ, xz])[x]] ≤  (Subscript[τ, 0]), 0, 
              If[(Subscript[τ, xz])[
                    x] > (Subscript[τ, 0]), ((Subscript[τ, 0]) - (Subscript[τ, xz])[
                        x])/(Subscript[μ, 0]), -((Subscript[τ, 0]) + (Subscript[τ, xz])[
                          x])/(Subscript[μ, 0])]], (Subscript[τ, xz])[0] ==  0, 
          (Subscript[v, z])[0] == Ω};
myODEsoln[Ω_] := 
      NDSolve[system[Ω], {(Subscript[v, z]), (Subscript[τ, xz])}, {x, 
          0, δ}]
yend[Ω_?NumericQ] := 
      Flatten[((Subscript[v, z])[r] /. myODEsoln[Ω]) /. 
          r -> δ]
bc = FindRoot[
            yend[Ω] == 0, {Ω, 0, 
              0.5}][[1, 2]];

sol1 = myODEsoln[bc]

plt1 = 
      Plot[(Subscript[v, z])[x] /. sol1, {x, 0, δ}, 
        PlotStyle -> RGBColor[0, 0, 1]]

plt2 = 
      Plot[(Subscript[τ, xz])[x] /. sol1, {x, 0, δ}, 
        PlotStyle -> RGBColor[0, 0, 1]]

Comparison with Newtonian fluid velocity profile obtained using the analytical solution. Velocity profile for the Bingham fluid is flat near the surface of the liquid film. In fact, only when (Subscript[τ, xz])[x] > (Subscript[τ, 0]), we have a non zero velocity gradient. This behavior is typical of Bingham fluids.

(Subscript[v, analytical])[
        x_] := ρ g δ^2/(2 (Subscript[μ, 0])) (1 - (x/δ)^2)

plt3 = 
      Plot[(Subscript[v, analytical])[x], {x, 0, δ}, 
        PlotStyle -> RGBColor[1, 0, 0]]

Show[plt3, plt1]