Section 9.8 Forced Convection

pp. 291-295

>  restart;
>  vz:=r->vzmax*(1-(r/R)^2); eq. 9.8-1
>  C:=r->2*Pi*r; The cicumference 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)-T0)*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. 
We now get D[2] which implies differentiation wrt the second 
variable: z.
>  de2:=simplify(de2=0);
>  qr:=(r,z)->-k*D[1](T)(r,z); qz:=0; eq. 9.8-10 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. 9.8-12
>  trans:={r=R*xi,z=zeta*rho*Cp*vzmax*R^2/k,T=T0+Theta*q1*R/k}; All three transformations in eqs. 9.8-16,17,18
with(PDEtools,dchange);
>  newde:=dchange(trans,de2,[Theta(xi,zeta),xi,zeta]);
>  newde:=simplify(%*R/q1);
>  nops(newde); This is one way to get the LHS of newde
>  op(1,newde);
>  debook:=-(1-xi^2)*D[2](Theta)(xi,zeta)+
diff(xi*diff(Theta(xi,zeta),xi),xi)/xi; The difference 
between the RHS & the LHS of eq. 9.8-19
>  op(1,newde)-debook;
>  simplify(%); Our equation agrees with the text.
>  Theta:=(xi,zeta)->C0*zeta+psi(xi); Assuming we want the 
solution for large zeta
>  newde;
>  s:=dsolve({newde,psi(0)=finite,D(psi)(1)=-1},psi(xi)); 
This does not work
>  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.
>  psi(xi);
>  Theta(xi,zeta);
>  eq:=-zeta=int(Theta(xi,zeta)*(1-xi^2)*xi,xi=0...1); eq. 9.8-25
>  _C1:=solve(eq,_C1);
>  Theta(xi,zeta); This agrees with eq. 9.8-31

>  Theta:='Theta'; Let's see if we can transform BC#2 from dimensional 
to dimensionless form using the same tranformation as the one used on the DE.
>  trans;
>  BC2:=-k*D[1](T)(r,z)=q1; eq. 9.8-14
>  newBC2:=simplify(dchange(trans,BC2,[Theta(xi,zeta),xi,zeta])/q1); 
Finally we get eq. 9.8-21

>