ABSTRACT (original article): We study four non-Newtonian fluid mechanics problems using Mathematica[R]. Constitutive equations describing the behavior of power-law, Bingham and Carreau models are recalled. The velocity profile is obtained for the horizontal flow of power-law fluids in pipes and annuli. For the vertical laminar film flow of a Bingham fluid we determine the velocity profile. Both problems involve the use of the shooting techniques because they have split boundary conditions. Since Mathematica permits symbolic computations, we determine analytical expressions of volumetric flow rates for pipe flow of the Bingham and power-law fluids. We use the built-in optimization command of Mathematica, FindMinimum, in order to find the non-Newtonian fluid model from representative data of flow rates measured under different applied pressure gradients in a horizontal pipe. These pedagogic problems are used to introduce the field of non-Newtonian fluid mechanics to students at the National Institute of Applied Sciences in Tunis. The Mathematica[R] notebooks are available from the corresponding author upon request or at Wolfram Research. (Contains 1 table and 6 figures). CITATION (original article): Binous, Housam (2007). “Introducing Non-Newtonian Fluid Mechanics Computations with Mathematica in the Undergraduate Curriculum.” Chemical Engineering Education 41(1):59–64. ISSN-0009-2479. https://www.researchgate.net/publication/237109031
We propose to find the velocity profiles for the laminar flow of Bingham fluids in a pipe (as depicted in below). Two numerical methods are proposed: the shooting method and the orthogonal collocation technique. We also compare the numerical solution with the analytical solution (see Appendix G of Ref. 2). We use the following values for the pressure difference ΔP, the pipe length L, and radius R: ΔP =100; L = 50; and R = 0.02
Parameters for Bingham fluid are taken as follows:
τ
0
=0.005;
μ
0
=0.5;
In[]:=
Clear[r,Ds,d,y,V,v,tau,T,τ,μ]

Flow of Bingham fluids in a horizontal pipe.

Bingham Fluids: solution based on the shooting technique
The blue dashed vertical line gives the loci where the velocity profile is flat (i.e., the velocity gradient is nil).
In[]:=
τ
0
=0.005;
μ
0
=0.5;ΔP=100;L=50;R=0.02;
In[]:=
system[Ω_]={D[r
τ
rz
[r],{r,1}]ΔP/Lr,D[
v
z
[r],{r,1}]If[Abs[
τ
rz
[r]]≤
τ
0
,0,If[
τ
rz
[r]>
τ
0
,(
τ
0
-
τ
rz
[r])/
μ
0
,-(
τ
0
+
τ
rz
[r])/
μ
0
]],
τ
rz
[10^-5]0,
v
z
[10^-5]Ω};​​myODEsoln[Ω_]:=NDSolve[system[Ω],{
v
z
,
τ
rz
},{r,10^-5,R}]​​yend[Ω_?NumericQ]:=Flatten[(
v
z
[r]/.myODEsoln[Ω])/.rR]​​bc=FindRoot[yend[Ω]0,{Ω,0.005}][[1,2]];
In[]:=
sol1=myODEsoln[bc]
Out[]=

v
z
InterpolatingFunction
Domain: {{0.00001,0.02}}
Output: scalar
,
τ
rz
InterpolatingFunction
Domain: {{0.00001,0.02}}
Output: scalar

In[]:=
pos=x/.FindRoot[Abs[
τ
rz
[x]/.sol1]==
τ
0
,{x,0.0005}]
Out[]=
0.00500002
In[]:=
plt1=Plot[
v
z
[r]/.sol1,{r,10^-5,R},PlotStyle{Thick,Dashed,Darker@Red},Frame->True,AspectRatio->1,FrameLabel->{Style["position",14],Style["velocity",14]},GridLines->Automatic,ImageSize->350{1,1},Epilog->{Darker@Blue,Dashed,Line[{{pos,0},{pos,0.005}}]}]
Out[]=
In[]:=
rY=2
τ
0
L/ΔP(*sameasposdefinedabove*)
Out[]=
0.005
In[]:=
z=(r-rY)/(R-rY);
In[]:=
v
analytic
[r_]:=If[z>0,Vmax(1-z^2),Vmax](*analyticalsolutionfromAppendixGofRef.2*)
In[]:=
(*maximalvelocityobtainedusingthetwomethods*)
In[]:=
Vmax=ΔP/4/L/
μ
0
(R-rY)^2(*VmaxfromAppendixGofRef.2*)
Out[]=
0.000225
In[]:=
Quiet@FindMaximum[{First[
v
z
[r]/.sol1]/.r->a,a>0},{a,0.01}](*Vmaxfromshootingtechniquesolution*)
Out[]=
{0.000224999,{a0.00502434}}
In[]:=
ShowPlot
v
analytic
[r],{r,0,R},PlotStyle->Blue,Frame->True,AspectRatio->1,FrameLabel->{Style["position",14],Style["velocity",14]},GridLines->Automatic,ImageSize->350{1,1},Epilog->{Darker@Blue,Dashed,Line[{{pos,0},{pos,0.005}}]},plt1
Out[]=
Comparing the solution obtained using the shooting technique (dashed red curve) and the solution (black dots) using Chebyshev Orthogonal Collocation for Bingham fluids
In[]:=
Np=40;​​Ds=Table[d[i,j],{i,0,Np},{j,0,Np}];​​d[0,0]=(2Np^2+1)/6//N;​​d[Np,Np]=-(2Np^2+1)/6//N;​​d[i_,i_]:=-y[i]/(2(1-y[i]^2));​​y[i_]:=Cos[iPi/Np]//N;​​d[j_,k_]:=If[j≠k,c[j](-1)^(j+k)/(c[k](y[j]-y[k]))];​​c[i_]:=If[i0||iNp,2,1];
In[]:=
V=Table[v[i],{i,1,Np+1}]//Flatten;
In[]:=
T=Table[tau[i],{i,1,Np+1}]//Flatten;
In[]:=
r=Table[R/2(1+y[i]),{i,Np,0,-1}]
Out[]=
{0.,0.0000308267,0.000123117,0.000276301,0.000489435,0.000761205,0.00108993,0.0014736,0.00190983,0.00239594,0.00292893,0.00350552,0.00412215,0.00477501,0.0054601,0.00617317,0.00690983,0.00766555,0.00843566,0.00921541,0.01,0.0107846,0.0115643,0.0123345,0.0130902,0.0138268,0.0145399,0.015225,0.0158779,0.0164945,0.0170711,0.0176041,0.0180902,0.0185264,0.0189101,0.0192388,0.0195106,0.0197237,0.0198769,0.0199692,0.02}
In[]:=
Eq1[i_]:=2/R(Ds.V)[[i]]If[Abs[tau[i]]≤
τ
0
,0,If[tau[i]>
τ
0
,(
τ
0
-tau[i])/
μ
0
,-(
τ
0
+tau[i])/
μ
0
]]
In[]:=
Eq2[i_]:=r[[Np+2-i]]2/R(Ds.T)[[i]]+tau[i]==r[[Np+2-i]]ΔP/L;
In[]:=
BC1=v[1]==0;BC2=tau[Np+1]==0;
In[]:=
unknowns=Flatten[Join[{Table[{V[[i]],0.5},{i,1,Np+1}],Table[{T[[i]],0.5},{i,1,Np+1}]}],1];
In[]:=
SYS=Flatten[Join[{BC1,BC2,Table[Eq1[i],{i,2,Np+1}],Table[Eq2[i],{i,1,Np}]}],1];
In[]:=
sol=FindRoot[SYS,unknowns]//Chop
Out[]=
{v[1]0,v[2]9.21099×
-7
10
,v[3]3.67834×
-6
10
,v[4]8.20993×
-6
10
,v[5]0.0000144435,v[6]0.0000222539,v[7]0.0000315101,v[8]0.0000420337,v[9]0.0000536475,v[10]0.0000661349,v[11]0.0000792894,v[12]0.0000928741,v[13]0.000106672,v[14]0.000120447,v[15]0.00013399,v[16]0.000147084,v[17]0.000159549,v[18]0.000171203,v[19]0.00018191,v[20]0.000191535,v[21]0.0002,v[22]0.000207227,v[23]0.000213197,v[24]0.000217891,v[25]0.000221353,v[26]0.000223621,v[27]0.000224786,v[28]0.000225044,v[29]0.000225018,v[30]0.000225028,v[31]0.000225023,v[32]0.000225026,v[33]0.000225024,v[34]0.000225026,v[35]0.000225025,v[36]0.000225025,v[37]0.000225025,v[38]0.000225025,v[39]0.000225025,v[40]0.000225025,v[41]0.000225025,tau[1]0.02,tau[2]0.0199692,tau[3]0.0198769,tau[4]0.0197237,tau[5]0.0195106,tau[6]0.0192388,tau[7]0.0189101,tau[8]0.0185264,tau[9]0.0180902,tau[10]0.0176041,tau[11]0.0170711,tau[12]0.0164945,tau[13]0.0158779,tau[14]0.015225,tau[15]0.0145399,tau[16]0.0138268,tau[17]0.0130902,tau[18]0.0123345,tau[19]0.0115643,tau[20]0.0107846,tau[21]0.01,tau[22]0.00921541,tau[23]0.00843566,tau[24]0.00766555,tau[25]0.00690983,tau[26]0.00617317,tau[27]0.0054601,tau[28]0.00477501,tau[29]0.00412215,tau[30]0.00350552,tau[31]0.00292893,tau[32]0.00239594,tau[33]0.00190983,tau[34]0.0014736,tau[35]0.00108993,tau[36]0.000761205,tau[37]0.000489435,tau[38]0.000276301,tau[39]0.000123117,tau[40]0.0000308267,tau[41]0}

Acknowledgement

Since in this laminar flow problem, we have split boundary conditions. Problem’s resolution, based on the Shooting Method, is derived based on coding developed in Mathematica by Professor Brian Higgins (University of California, Davis).

References

◼
  • J. O. Wilkes, Fluid Mechanics for Chemical Engineers, Upper Saddle River, NJ: Prentice Hall, 1999.
    • ◼
  • John C. Dixon, The Shock Absorber Handbook, Appendix G: Bingham Flow, NY: Wiley, 2007.

    • CITE THIS NOTEBOOK

    Velocity profile for laminar flow of Bingham fluids using orthogonal collocation​
    by André Von-Held Soares, Lizandro de Sousa Santos, Housam Binous and Ahmed Bellagi
    Wolfram Community, STAFF PICKS, April 19, 2024
    ​https://community.wolfram.com/groups/-/m/t/3160813