Obtaining the Velocity Profile for Pipe Flow of a Non-Newtonian Fluid

> restart;

Area in r-direction for a pipe

> A1:=r->2*Pi*r*L;

Area in z-direction for a pipe

> A2:=r->2*Pi*r*dr;

Left-hand side of momentum balance due to shear stress

> LHS:=A1(r)*tau[rz](r)-A1(r+dr)*tau[rz](r+dr);

Right-hand side of momentum balance due to pressure forces

> RHS:=A2(r)*pL-A2(r)*p0;

Take limit as dr->0 to obtain differential equation

> eq1:=limit(LHS/dr,dr=0)=RHS/dr:

> eq1:=simplify(-eq1/2/Pi/L);

Solve the differential equation to obtain an expression for shear stress as a function of r

> s1:=dsolve(eq1,tau[rz](r)); assign(s1); tau[rz]:=unapply(tau[rz](r),r);

Take the constant of integration to be zero, to make sure that the shear stress does not blow up at r=0

> _C1:=0;

Ellis Model for non-Newtonian fluids

> eq2:=-D(vz)(r)=(phi[0]+phi[1]*abs(tau[rz](r))^(alpha-1))*tau[rz](r);

For given conventions, the shear stress is positive

> assume(tau[rz](r),positive);

> eq2:=simplify(eq2):

Solve differential equation for velocity profile as function of radius, with the boundary condition of zero velocity at pipe wall

> s2:=dsolve({eq2,vz(R)=0},vz(r)): assign(s2); vz:=unapply(vz(r),r);

Velocity at the wall, check expression for velocity in z-direction

> simplify(vz(R));

Solve for maximum velocity at centerline

> vmax:=simplify(vz(0));