Example 10.5-5 One Dimensional Compressible Flow

>  restart;
>  econt:=D(rhovx)(x)=0; Continuity
>  emot:=rhovx(x)*D(vx)(x)=-D(p)(x)+(4/3)*mu*D(D(vx))(x); 
Motion
>  een:=rhovx(x)*Cp*D(T)(x)=k*D(D(T))(x)+vx(x)*D(p)(x)+
(4/3)*mu*(D(vx)(x))^2; Energy
>  s:=dsolve({econt,rhovx(0)=v1*rho1},rhovx(x)); rho*vx is 
constant
>  assign(s);rhovx:=unapply(rhovx(x),x);
>  een2:=lhs(emot)*vx(x)+lhs(een)=rhs(emot)*vx(x)+rhs(een); 
Getting rid of the pressure term in the energy equation by using 
the equation of motion.
>  een2:=simplify(%/(v1*rho1));
>  int(een2,x); We can not integrate it.
>  int(lhs(een2),x); We can't even integrate one side of the 
equation.
>  int(D(vx)(x)*vx(x)+Cp*D(T)(x),x);
>  int(rhs(een2),x);
>  diff(Cp*T(x)+(1/2)*vx(x)^2,x)-lhs(een2); Going the other 
direction to check on the validity of BS&L's answer in 10.5-57
>  simplify(%); The left side checks.
>  dif:=k/(Cp*rho1*v1)*diff((4/3)*Pr*D((vx^2)/2)(x)+
Cp*D(T)(x),x)-rhs(een2);
>  simplify(%);
>  Cp:=Pr*k/mu; Replacing Cp with the Prandtl Number
>  simplify(dif); The rest of 10.5-57 checks.

>  Pr:=3/4; Suggested as the only value for which we can get a solution.
>  een2; Here is our equation.
>  T:=x->(C1+C2*exp(rho1*v1*Cp*x/k)-vx(x)^2/2)/Cp; Here is 
the solution in 10.5-58
>  een2; Is this OK?
>  simplify(%); Yes. The solution given as 10.5-58 in the text 
satisfies the energy equation.