Example 19.1-1 Unsteady-State Evaporation in a Binary System

>  restart;
>  eq1:=D[1](cA)(t,z)=-diff(NAz(t,z),z); eq. 19.1-1. 
(using diff instead of D or the RHS seems to be required in 
order to get eq7 a little further along.
>  eq2:=D[1](cB)(t,z)=-diff(NBz(t,z),z); eq. 19.1-2
>  eqc:=eq1+eq2; The total molar concentration does not 
change so the left side of this should be zer0.
>  cB:=(t,z)->c-cA(t,z); We can tell Maple that c is a constant 
by:
>  eqc; Then:
>  NBz:=(t,z)->Ntotz(t,z)-NAz(t,z); Now replace NBz with the 
total flux.
>  eqc; Eq. c becomes simple enough that:
>  s:=pdsolve(eqc,Ntotz); pdsolve can find a solution.
>  Ntotz:=(t,z)->NA0(t); We will call _F1(t), NA0(t)
>  NA0:=t->-c*DAB*D[2](xA)(t,0)/(1-xA0); Eq. 19.1-4
>  NAz:=(t,z)->-c*DAB*D[2](xA)(t,z)+xA(t,z)*Ntotz(t,z); 
eq. 17.0-1
>  cA:=(t,z)->c*xA(t,z); For an ideal gas.
>  eq7:=simplify(eq1/c); eq. 19.1-7 (except for a few negative 
signs in the top and bottom.)
>  eqX:=D(D(X))(Z)+2*(Z-phi)*D(X)(Z); Using dchange to switch 
variables in eq7 has not been sucessful so far. Thus we will see 
if the solution to this dimensionless equation does satisfy eq7. 
This is eq. 19.1-12
>  s:=dsolve({eqX,X(0)=1,X(infinity)=0},X(Z)); Solving 19.1-12 
together with the BC: 19.1-14 and 15.
>  assign(s);X:=unapply(X(Z),Z);
>  xA:=(t,z)->xA0*X(z/sqrt(4*DAB*t)); Using the definition of 
X in 19.1-11 to recover xA.
>  eq7; Unfortunately, the relation for xA does not get inserted 
into this.
>  xA(t,z);
>  eq7a:=diff(xA(t,z),t)-DAB*diff(xA(t,z),z,z)-DAB*diff(xA(t,z),
z)*Dif/(1-xA0); We need to rewrite eq7 in terms of diff instead 
of D. At the same time we will put all terms on the same side of 
the equation.
>  Dif:=unapply(diff(xA(t,z),z),z);
>  Dif:=Dif(0); This is the long way to get D[2](xA)(t,0).
>  eq7a:=simplify(eq7a,assume=positive); This is not zero since 
we have not used the realtion between phi and xA0 given in 
eq. 19.1-13.
>  -(1/2)*xA0*D(X)(0)/(1-xA0); This is the RHS of eq. 19.1-13. 
It is a simple function of xA0, but a rather complicated on in phi.
>  eqphi:=phi=%; Eq. 19.1-13
>  xA0:=solve(eqphi,xA0); There is no hope of solving for phi, 
but this works.
>  simplify(eq7a,assume=positive); Eq. 7a is satisfied!
>  xA0f:=unapply(xA0,phi); Here is xA0 as a function of phi. 
It should be identical to eq. 19.1-17. We will show that it gives 
the same results shown in Table 19.1-1
>  evalf(xA0f(.1562)); This checks.
>  evalf(xA0f(.3578)); So does this.
>  evalf(xA0f(.6618)); And this.
>  NA0(t); This has the same problem with D as before.
>  Dif:=unapply(diff(xA(t,z),z),z);
>  Dif(0);
>  NA0:=t->-c*DAB*Dif(0)/(1-xA0);
>  NA0(t); This only depends on phi and t.
>  VA:=simplify(int(NA0(t)/c,t=0...t1),assume=positive); 
eq. 19.1-19 omitting the area S.
>  VAFick:=xA0*sqrt(4*DAB*t1/Pi); eq. 19.1-20
>  psi2:=simplify(VA/VAFick,assume=positive); The binary 
correction coefficient in 19.1-21.
>  simplify(%-phi*sqrt(Pi)/xA0,assume=positive); It agrees 
with the formula in Table 19.1-1

>