> restart;

This problem deals with steady-state laminar flow in which the fluid is moving in a circular pattern. Therefore, with cylindrical coordinates, the velocity components vr and vz are zero. In addition, there is no pressure gradient in the theta direction.

All the terms in the equation of continuity: D(rho)(t)+1/r*D(rho*r*vr)(r)+1/r*D(rho*vtheta)(theta)+D(rho*vz)(z)=0

are equal to 0

The equations of motion found in Table 3.4-4 reduce to:

> eqr:=-rho*vtheta(r)^2/r=-diff(p(r),r);

> eqtheta:=0=diff(1/r*diff(r*vtheta(r),r),r);

> eqz:=0=-diff(p(z),z)+rho*gz;

> s:=dsolve({eqtheta,vtheta(kappa*R)=0, vtheta(R)=Omega*R},vtheta(r));

> assign(s); vtheta:=unapply(vtheta(r),r);

> vtheta(r):=simplify(vtheta(r));

> vthetabook:=Omega*R*(r/(kappa*R)-kappa*R/r)/(1/kappa-kappa);

> eq1:=vtheta(r)-vthetabook;

> simplify(%);

This results in the velocity distribution of Eq 10.5-4

The energy equation found as Eq. B in Table 10.2-3 reduces to:

> en:=0=k/r*diff(r*diff(T(r),r),r)+mu*(r*diff(vtheta(r)/r,r))^2;

This is the differential equation for the temperature distribution.

> s1:=dsolve({en,T(0)=T[o],T(b)=T[b]},T(r));

On trying to solve equation en, our results were too complex, thus we are switching to dimensionless units as shown below:

> with(PDEtools,dchange);

> transf:={r=R*xi, T=theta*(T1-Tk)+Tk};

> Newde:=dchange(transf,R^2*en/(k*(T1-Tk)),[theta(xi),xi]);

> simplify (Newde);

Introducing N into our equation:

> mu:=N*k*(T1-Tk)*(1-kappa^2)^2/(Omega^2*R^2*kappa^4);

> en1:=simplify(Newde);

> s1:=dsolve({en1,theta(kappa)=0,theta(1)=1},theta(xi));

> assign(s1); theta:=unapply(theta(xi),xi);

> books1:=((N+1)-N/xi^2)-((N+1)-N/kappa^2)*log(xi)/log(kappa);

> simplify(theta(xi)-books1);

This shows that theta(xi) is the same as shown in Equation 10.5-14