Section 10.8 Forced Convection
>  restart;
>  vz:=r->vzmax*(1-(r/R)^2); eq. 10.8-1
>  C:=r->2*Pi*r; The circumference of a circle of radius r
>  Qr:=(r,z)->C(r)*qr(r,z)*dz; Heat flow in radial direction through the surface of a ring at position (r,z)
>  Qz:=(r,z)->C(r)*qz(r,z)*dr; Heat flow in axial direction through the surface of a ring at position (r,z)
>  Qconv:=(r,z)->rho*Cp*vz(r)*(T(r,z)-To)*C(r)*dr; Heat flow in z direction by convection through the surface of a ring
>  eq:=Qr(r,z)-Qr(r+dr,z)+Qz(r,z)-Qz(r,z+dz)+Qconv(r,z)-Qconv(r,z+dz); Energy balance
>  de:=simplify(eq/(C(r)*dz*dr));
>  de1:=limit(de,dr=0); Taking the limit as dr goes to zero. Note that D[1] indicates differentiation wrt the first argument: r.
>  de2:=-limit(de1,dz=0); Taking the limit as dz goes to zero and changing the sign of the equation. We now get D[2] which implies differentiation wrt the second variable: z.
>  de2a:=(diff(r*qr(r,z),r)+r*diff(qz(r,z),z)+r*Cp*rho*vz(r)*diff(T(r,z),z))/r; the eq. in ed 1 of BS&L
>  simplify(de2-de2a);
>  qr:=(r,z)->-k*D[1](T)(r,z); qz:=0; Fourier's Law for the radial component of the heat flux, but zero for the axial term since we will neglect it.
>  de2; This is the dimensional form of the PDE. It should be equivalent to eq. 10.8-12
>  de2b:=rho*Cp*vz(r)*diff(T(r,z),z)-k/r*diff(r*diff(T(r,z),r),r);
>  simplify(de2-de2b);
>  trans:={r=R*xi,z=zeta*rho*Cp*vzmax*R^2/k,T=T1+Theta*q0*R/k}; All three transformations in eqs. 10.8-16,17,18
>  with(PDEtools,dchange);

>  newde:=dchange(trans,de2,[Theta(xi,zeta),xi,zeta]);
>  newde:=simplify(%*R/q0);
>  debook:=(1-xi^2)*D[2](Theta)(xi,zeta)-diff(xi*diff(Theta(xi,zeta),xi),xi)/xi; The difference between the LHS & the RHS of eq. 10.8-19
>  newde-debook;
>  simplify(%); Our equation agrees with the text.
>  BC2:=k*diff(T(r,z),r)=q0; eq.10.8-14 at r=R
>  newBC2:=dchange(trans,BC2,[Theta(xi,zeta),xi,zeta])/q0; eq. 10.8-21 for xi=1
>  Theta:=(xi,zeta)->C0*zeta+psi(xi); Assuming we want the solution for large zeta
>  newde;
>  newBC2;
>  s:=dsolve({newde,psi(0)=finite,D(psi)(1)=1},psi(xi)); This does not give a reasonable answer.
>  s:=dsolve({newde,D(psi)(1)=1},psi(xi)); Leaving out the BC at xi = 0.
>  assign(s); psi:=unapply(psi(xi),xi);
>  C0:=4; To keep the solution finite at xi = 0 ( ie r=0) by setting the coef. of log(xi)=0
>  psi(xi);
>  Theta(xi,zeta);
>  trans;
>  T:=(r,z)->T1+Theta(r/R,z*k/(rho*Cp*vzmax*R^2))*q0*R/k;
>  simplify(T(r,z),assume=positive);
>  eq:=rho*Cp*int((T(r,L)-T1)*vz(r)*2*Pi*r,r=0...R)=2*Pi*R*L*q0; using 10.8-24
>  simplify(solve(eq,_C1));
>  eq:=zeta=int(Theta(xi,zeta)*(1-xi^2)*xi,xi=0...1); eq. 10.8-25
>  _C1:=solve(eq,_C1);
>  Theta(xi,zeta); This agrees with eq. 10.8-31

>