Sectionn 10.1: The Equations of Energy: looking at the derivation of terms in eq. 10.1.8

>  restart;
with(linalg): Using the linear algebra package
>  Ax:=dy*dz; Ay:=dx*dz; Az:=dx*dy; dV:=dx*dy*dz; 
Areas perpendicular to each of the axes and the volume of the element 
v:=(x,y,z,t)->vector([vx(x,y,z,t),vy(x,y,z,t),vz(x,y,z,t)]); 
velocity vector
>  q:=(x,y,z,t)->vector([qx(x,y,z,t),qy(x,y,z,t),qz(x,y,z,t)]); 
heat conduction vector
>  g:=(x,y,z,t)->vector([gx(x,y,z,t),gy(x,y,z,t),gz(x,y,z,t)]); 
gravitational vector
>  IKE:=(x,y,z,t)->rho(x,y,z,t)*(U(x,y,z,t)+Vsq(x,y,z,t)/2); 
internal plus kinetic energy per volume of fluid
>  Accum:=dV*(IKE(x,y,z,t+dt)-IKE(x,y,z,t)); accumulation in the 
volume element
>  LHS:=limit(Accum/(dt*dV),dt=0); Getting the left side of eq. 10.1-8. 
Note these are the same terms, but the derivative of rho times U and 1/2v^2 
is expanded by Maple.


>  simplify(LHS-D[4](rho*(U+Vsq/2))(x,y,z,t)); Subtracting off the 
LHS of eq. 10.1-8

>  CCE:=(x,y,z,t)->matadd(q(x,y,z,t),v(x,y,z,t),1,IKE(x,y,z,t)); 
Conduction plus convection fluxes
>  CCEx:=(x,y,z,t)->CCE(x,y,z,t)[1]*Ax;  Flow of convection plus 
conduction across an element perpendicular to the x axis
>  CCEy:=(x,y,z,t)->CCE(x,y,z,t)[2]*Ay;  Flow of convection plus 
conduction across an element perpendicular to the y axis
>  CCEz:=(x,y,z,t)->CCE(x,y,z,t)[3]*Az;  Flow of convection plus 
conduction across an element perpendicular to the z axis
>  CCEx(x,y,z,t);
>  dCCEx:=limit((CCEx(x,y,z,t)-CCEx(x+dx,y,z,t))/dV,dx=0);
>  simplify(dCCEx+D[1](vx*rho*(U+Vsq/2)+qx)(x,y,z,t));

>  GWork:=(x,y,z,t)->rho(x,y,z,t)*multiply(v(x,y,z,t),g(x,y,z,t))*dV; 
work associated with gravity.

>  GWork(x,y,z,t)/dV; This checks with the gravity term in eq. 10.1-8

>  tau:=(x,y,z,t)->matrix([[txx(x,y,z,t),txy(x,y,z,t),txz(x,y,z,t)],
[tyx(x,y,z,t),tyy(x,y,z,t),tyz(x,y,z,t)],[tzx(x,y,z,t),tzy(x,y,z,t),
tzz(x,y,z,t)]]); The shear or momentum transport matrix or tensor
>  TaupV:=(x,y,z,t)->matadd(multiply(tau(x,y,z,t),v(x,y,z,t)),
v(x,y,z,t),1,p(x,y,z,t)); The three components of the rate of work per 
area on an element perpendicular to each face.
This includes the effect of pressure in addition to the shear forces.
>  TaupV(x,y,z,t)[1]; The first component of the work vector.

>  limit((TaupV(x+dx,y,z,t)[1]-TaupV(x,y,z,t)[1])*Ax/dV,dx=0); 
Getting the terms in eq. 10.1-8 that arise from work term differences in 
the x direction.
>  simplify(%-D[1]((p+txx)*vx+txy*vy+txz*vz)(x,y,z,t)); Checking with 
the terms in 10.1-8.