Example 11.4-1: Steady-State Forced-Convection Heat Transfer in Laminar Flow in a Circular Tube

Problem:

Show how to set up the equations for the problem considered in section 10.8 - namely, that of finding the fluid temperature profiles for the fully developed laminar flow in a tube.

We will start with the general equations of continuity, motion, and energy in cylindrical coordinates, make the necessary assumptions, and then show that our result is that obtained in section 10.8.

> restart;

> s:=vector([r,Theta,z]); "S" is the position vector in cylindrical coordinates

s := vector([r, Theta, z])

> rho:=(s,t)->rho; v:=vector([vr(r,Theta,z,t), vTheta(r,Theta,z,t), vz(r,Theta,z,t)]); mu:=(s,t)->mu; The assumption is made that the fluid has constant properties, i.e., constant density and constant viscosity; Furthermore, velocity has three components, each of which can depend on space and time.

rho := proc (s, t) options operator, arrow; rho end...

v := vector([vr(r,Theta,z,t), vTheta(r,Theta,z,t), ...

mu := proc (s, t) options operator, arrow; mu end p...

> d[r]:=f->diff(f,r); d[Theta]:=f->diff(f,Theta); d[z]:=f->diff(f,z); d[t]:=f->diff(f,t); These functions represent the differential operators in space and time.

d[r] := proc (f) options operator, arrow; diff(f,r)...

d[Theta] := proc (f) options operator, arrow; diff(...

d[z] := proc (f) options operator, arrow; diff(f,z)...

d[t] := proc (f) options operator, arrow; diff(f,t)...

> d1:=d[r](v[1]): d2:=d[Theta](v[2]): d3:=d[z](v[3]): d1,d2, and d3 represent the derivatives of each velocity component with respect to its direction, e.g., the derivative with respect to r of vr.

> continuity:=d[t](rho(s,t))+rho(s,t)*((1/r)*d[r](r*v[1])+(1/r)*d2+d3); General equation of continuity for cylindrical coordinates, expanded, with the assumption of constant physical properties, as given in Appendix B, eq. B.4-2 in B,S,&L.

continuity := rho*(1/r*(vr(r,Theta,z,t)+r*diff(vr(r...

> motion:=rho(s,t)*(d[t](v[3])+v[1]*d[r](v[3])+(v[2]/r)*d[Theta](v[3])+v[3]*d3)+d[z](P(r,Theta,z,t))-mu(s,t)*((1/r)*d[r](r*d[r](v[3]))+(1/r^2)*d[Theta](d[Theta](v[3]))+d[z](d3)); General equation of motion in cylindrical coordinates, as given in Appendix B, eq. B.6-6, assuming constant physical properties. Gravitational term rho*gz(z) is combined with the pressure term d(p(z))/dz by the term P(z), which is P =p(z)-rho*g*z, a convenient way to sum the gravitational and pressure terms in balances.

motion := rho*(diff(vz(r,Theta,z,t),t)+vr(r,Theta,z...
motion := rho*(diff(vz(r,Theta,z,t),t)+vr(r,Theta,z...
motion := rho*(diff(vz(r,Theta,z,t),t)+vr(r,Theta,z...

> energy:=rho(s,t)*Cp*(d[t](T(r,Theta,z,t))+v[1]*d[r](T(r,Theta,z,t))+v[2]/r*d[Theta](T(r,Theta,z,t))+v[3]*d[z](T(r,Theta,z,t)))-k*((1/r)*d[r](r*d[r](T(r,Theta,z,t)))+(1/r^2)*d[Theta](d[Theta](T(r,Theta,z,t)))+d[z](d[z](T(r,Theta,z,t))))-mu(s,t)*Phi(r,Theta,z,t); General energy equation, expanded, assuming constant physical properties, as given in Appendix B, eq. B.9-2

energy := rho*Cp*(diff(T(r,Theta,z,t),t)+vr(r,Theta...
energy := rho*Cp*(diff(T(r,Theta,z,t),t)+vr(r,Theta...
energy := rho*Cp*(diff(T(r,Theta,z,t),t)+vr(r,Theta...

> Phi:=(r,Theta,z,t)->2*(d1)^2+2*(d2/r+v[1]/r)^2+2*(d3^2)+(r*d[r](v[2]/r)+(1/r)*d[Theta](v[1]))^2 +((1/r)*d[Theta](v[3])+d[z](v[2]))^2+(d[z](v[1])+d[r](v[3]))^2-(2/3)*((1/r)*d[r](r*v[1])+(1/r)*d2+d3)^2:

> Phi(r,Theta,z,t); Phi is the equation for the general dissipation function in cylindrical coordinates, as given in Appendix B by eq. B.8-2.

2*diff(vr(r,Theta,z,t),r)^2+2*(1/r*diff(vTheta(r,Th...
2*diff(vr(r,Theta,z,t),r)^2+2*(1/r*diff(vTheta(r,Th...
2*diff(vr(r,Theta,z,t),r)^2+2*(1/r*diff(vTheta(r,Th...
2*diff(vr(r,Theta,z,t),r)^2+2*(1/r*diff(vTheta(r,Th...

> rho:='rho': vr:='vr': vTheta:='vTheta': vz:='vz': p:='p': T:='T': mu:='mu': We need to free up these variable names so we can assign them again as functions to simplify our equations of continuity,motion, and energy.

We now make some assumptions to reduce these general equations. We want to assume that the only velocity component is in the z-direction, and that it varies only with radial position. We also assume that pressure varies only with height. Furthermore, we assume that the temperature T varies only with height and radial position:

> vr(r,Theta,z,t):=0: vTheta(r,Theta,z,t):=0: vz(r,Theta,z,t):=vz(r): T(r,Theta,z,t):=T(r,z): P(r,Theta,z,t):=P(z): mu(r,Theta,z,t):=mu:

> continuity; Eq. 11.4-1

0

> motion; Eq. 11.4-2, with the derivative of r*dvz/dr expanded.

diff(P(z),z)-mu/r*(diff(vz(r),r)+r*diff(vz(r),`$`(r...

> energy; Eq. 11.4-3, with the term which contains a second derivative of T expanded.

rho*Cp*vz(r)*diff(T(r,z),z)-k*(1/r*(diff(T(r,z),r)+...

Next, we need to solve for the velocity distribution from the equation of motion.

> sol:=dsolve(motion,vz(r));

sol := vz(r) = 1/4*diff(P(z),z)/mu*r^2+_C1*ln(r)+_C...

In order that the velocity be finite at r=0, the constant _C1 must be 0 to prevent the natural log term from blowing up:

> _C1:=0; sol;

_C1 := 0

vz(r) = 1/4*diff(P(z),z)/mu*r^2+_C2

> assign(sol); vzr:=unapply(vz(r),r);

vzr := proc (r) options operator, arrow; 1/4*diff(P...

In order to solve for the last constant, we use the Boundary Condition that at the pipe edge, the velocity is zero:

> bc2:=vzr(R)=0;

bc2 := 1/4*diff(P(z),z)/mu*R^2+_C2 = 0

> _C2:=solve(bc2,_C2); Solve for the constant.

_C2 := -1/4*diff(P(z),z)/mu*R^2

> simplify(vzr(r));

-1/4*diff(P(z),z)*(-r^2+R^2)/mu

Pressure varies with height in a hydrostatic manner, and the gravitational force varies monotonically with height: i.e., the derivative of P with respect to z the for the length of the pipe is simply -(P(0)-P(L))/L

> eqP:=diff(P(z),z)=-(P0-PL)/L:

> sol:=dsolve(eqP,P(z)):

> assign(sol): P:=unapply(P(z),z): vzr(r);

1/4*(-P0+PL)/L/mu*r^2-1/4*(-P0+PL)/L/mu*R^2

This is the parabolic velocity distribution expected, as used in section 10.8 as Eq. 10.8-1.

We now want to substitute this result into the energy equation.

But first, we need to make some assumptions about the equation of energy: i) heat conduction in the z-direction is much smaller than heat convection, so that we may neglect the second order differential of T, and ii) the viscous heating effects are minimal, so that we may ignore the mu*(d^2(vz)/dr^2) term:

> newenergy:=simplify(rho*Cp*vzr(r)*d[z](T(r,z))=k*(1/r)*d[r](r*d[r](T(r,z))));

newenergy := -1/4*rho*Cp*(-P0+PL)*(-r^2+R^2)/L/mu*d...

we will define vzmax as in eq. 2.3-19:

> mu:=(P0-PL)*R^2/4/L/vzmax: simplify(newenergy);

rho*Cp*(-r^2+R^2)/R^2*vzmax*diff(T(r,z),z) = k/r*(d...

Comparing this result with the equation 10.8-12:

> eq10812:=rho*Cp*vzmax*(1-(r^2/R^2))*d[z](T(r,z))=k*((1/r)*d[r](r*d[r](T(r,z))));

eq10812 := rho*Cp*vzmax*(1-r^2/R^2)*diff(T(r,z),z) ...

> difference:=simplify(eq10812-newenergy);

difference := 0 = 0

Thus, using the equations of motion to solve for the velocity distribution in the fully developed laminar region region, we get the same results as with the shell balance method.