Example 11.4-7 One Dimensional Compressible Flow
a second version
>  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(rhs(een2),x);
>  simplify(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

>  Cp:=Pr*k/mu;  Replacing Cp with the Prandtl Number
>  simplify(k/(Cp*rho1*v1)*diff((4/3)*Pr*D((vx^2)/2)(x)+Cp*D(T)(x),x)-rhs(een2)); The right hand side also checks out.

>  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.

>  T(x);
>  C2:=0;eq:=T(x1)=T1;
>  C1:=solve(eq,C1);
>  C1:=subs(vx(x1)=v1,C1); Remember that Pr=3/4=Cp*mu/k, so 3/4(k/mu)=Cp.
>  simplify(een2);

>  ien:=simplify((1/2)*vx(x)^2+Cp*T(x)=(1/2)*v1^2+Cp*T1); The integrated energy equation is now satisfied.

>  emot;

>  iemot:=x->v1*rho1*vx(x)+p(x)-(4/3)*mu*D(vx)(x)-C3;
>  D(iemot)(x); This can be seen to be identical to the motion equation.

>  C3:=p1+rho1*v1^2;
>  iemot(x);

>