Solving the Maxwell-Stefan equations using the orthogonal collocation 
and the shooting methods 
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

In the bottom of a Stefan tube there is a quiescent liquid with a mixture of acetone (component 1) and methanol (component 2). Vapor that evaporates from this pool of liquid diffuses to the top of the tube. Air (component 3) keeps the mole fractions of components 1 and 2 equal to zero. 

(Subscript[F, ij]) = ((Subscript[cD, ij])/δ)where c=(P/RT) , T=328.5 K , P=99.4 kPa , δ=0.24 m.

The Maxwell-Stefan diffusion coefficients of the binary pairs are : 
(Subscript[D, 12])= 8.48 (mm^2)/s,  (Subscript[D, 13])= 13.72 (mm^2)/s and  (Subscript[D, 23])= 19.91 (mm^2)/s.

c=99.4/(328.5 8.314) 10^3

δ = 0.24

(Subscript[D, 12]) = 8.48

(Subscript[F, 12]) = c/δ (Subscript[D, 12]) 10^-6

(Subscript[D, 13]) = 13.72

(Subscript[F, 12]) = c/δ (Subscript[D, 13]) 10^-6

(Subscript[D, 23]) = 19.91

(Subscript[F, 23]) = c/δ (Subscript[D, 23]) 10^-6

(Subscript[F, 12]) = 1.286 10^-3; (Subscript[F, 13]) = 2.081 10^-3; 
    (Subscript[F, 23]) = 3.019 10^-3;

Since air does not dissolve in the liquid, we have N3=0

N3=0;

At x=δ (top of the Stefan tube) , Air (component 3) keeps the mole fractions of components 1 and 2 equal to zero. At x=0 (the vapor liquid interface), the mole fractions of components 1 and 2 are equal to their equilibrium values at T=328.5 K and P=99.4 kPa (Carty and Schrodt, 1975).

Y[1,1]=0.319;Y[2,1]=0.528;Y[1,5]=0;Y[2,5]=0;

Orthogonal Collocation Method

Three internal collocation points (N=3)

Solve[JacobiP[3,1,1,-1+2 u]==0.,u]

x[1] = 1.0 10^-20; x[2] = 0.172673; x[3] = 0.5; x[4] = 0.827327; 
    x[5] = 1. ;

Matrix Q and (Q^(-1))

Q = 
      Chop[Table[Table[x[j]^(i - 1), {i, 1, 5}], {j, 1, 5}]]

IQ = Inverse[Q]

Matrix C

Ci = 
      Chop[Table[
          Table[(i - 1) x[j]^(i - 2), {i, 1, 5}], {j, 1, 
            5}]]

Matrix A

A=Ci.IQ

Ci . IQ // MatrixForm

Matrix D

Di = 
      Chop[Table[
          Table[(i - 1) (i - 2) x[j]^(i - 3), {i, 1, 5}], {j, 
            1, 5}]]

Matrix B

B=Di.IQ

Di . IQ // MatrixForm

Maxwell-Stefan relations written using the orthogonal collocation method gives a set of simultaneous algebraic equations that are easily solved using FindRoot.

For[i = 1, i < 6, 
      eq[1, i] = (∑_(j = 1)^5) A[[i, j]] Y[1, j] == 
          1/(Subscript[F, 12]) (Y[1, i] N2 - Y[2, i] N1) + 
            1/(Subscript[F, 13]) (Y[1, 
                      i] N3 - (1 - Y[1, i] - 
                        Y[2, i]) N1); i++]

For[i = 1, i < 6, 
      eq[2, i] = (∑_(j = 1)^5) A[[i, j]] Y[2, j] == 
          1/(Subscript[F, 12]) (Y[2, i] N1 - Y[1, i] N2) + 
            1/(Subscript[F, 23]) (Y[2, 
                      i] N3 - (1 - Y[1, i] - 
                        Y[2, i]) N2); i++]

sol=FindRoot[Flatten[Join[Table[eq[1,i],{i,2,5}],Table[eq[2,i],{i,2,5}]]],{Y[1,2],0.3},{Y[1,3],0.2},{Y[1,4],0},{Y[2,2],0.4},{Y[2,3],0.25},{Y[2,4],0},{N1,1.286 10^-3},{N2,1.286 10^-3}]

Plotting the results

Tbl1=Table[{x[i],Y[1,i]}/.sol,{i,1,5}]

Tbl2=Table[{x[i],Y[2,i]}/.sol,{i,1,5}]

Tbl3=Table[{x[i],1-Y[1,i]-Y[2,i]}/.sol,{i,1,5}]

plt1=ListPlot[Tbl1,PlotStyle->RGBColor[1,0,0],PlotJoined->True,DisplayFunction->Identity]

plt2=ListPlot[Tbl2,PlotStyle->RGBColor[0,0,1],PlotJoined->True,DisplayFunction->Identity]

plt3=ListPlot[Tbl3,PlotStyle->RGBColor[0,1,0],PlotJoined->True,DisplayFunction->Identity]

Show[plt1,plt2,plt3,DisplayFunction->$DisplayFunction]

plt4=ListPlot[Tbl1,PlotStyle->{PointSize[0.02],RGBColor[1,0,0]},DisplayFunction->Identity]

plt5=ListPlot[Tbl2,PlotStyle->{PointSize[0.02],RGBColor[0,0,1]},DisplayFunction->Identity]

plt6=ListPlot[Tbl3,PlotStyle->{PointSize[0.02],RGBColor[0,1,0]},DisplayFunction->Identity]

Show[plt4,plt5,plt6,DisplayFunction->$DisplayFunction]

Shooting Method

Acknowledgement : 
Since in this problem we have split boundary conditions, its resolution is based on the Shooting Method developped in Mathematica by Professor Brian Higgins (Chemical Engineering Department, University of California, Davis, http://www.higgins.ucdavis.edu/chemmath.php  ).

Solving the Maxwell-Stefan equations using the shooting method

system[N1_, 
        N2_] = {D[Y1[x], {x, 1}] == 
          1/(Subscript[F, 12]) (Y1[x] N2 - Y2[x] N1) + 
            1/(Subscript[F, 13]) (Y1[
                      x] N3 - (1 - Y1[x] - 
                        Y2[x]) N1), 
D[
            Y2[x], {x, 1}] == 
          1/(Subscript[F, 12]) (Y2[x] N1 - Y1[x] N2) + 
            1/(Subscript[F, 23]) (Y2[
                      x] N3 - (1 - Y1[x] - 
                        Y2[x]) N2), 
Y1[0] == 
          0.319, Y2[0] == 0.528}

myODEsoln[N1_,N2_]:=NDSolve[system[N1,N2],{Y1[x],Y2[x]},{x,0,1}]

yend1[N1_?NumericQ,N2_?NumericQ]:=Flatten[(Y1[x]/.myODEsoln[N1,N2])/.x->1]
yend2[N1_?NumericQ,N2_?NumericQ]:=Flatten[(Y2[x]/.myODEsoln[N1,N2])/.x->1]

bc=FindRoot[{yend1[N1,N2]==0,yend2[N1,N2]==0},{N1,0.05},{N2,0.057}];
c={N1,N2}/.bc;

c

Values of the fluxes N1 and N2

c[[1]]

c[[2]]

Plotting the results

plt7=Plot[Evaluate[Y1[x]/.myODEsoln[c[[1]],c[[2]]]],{x,0,1},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity]

plt8=Plot[Evaluate[Y2[x]/.myODEsoln[c[[1]],c[[2]]]],{x,0,1},PlotStyle->RGBColor[0,0,1],DisplayFunction->Identity]

plt9=Plot[Evaluate[(1-Y1[x]-Y2[x])/.myODEsoln[c[[1]],c[[2]]]],{x,0,1},PlotStyle->RGBColor[0,1,0],DisplayFunction->Identity]

Show[plt7,plt8,plt9,DisplayFunction->$DisplayFunction]

Comparing results from the shooting method and the orthogonal collocation method. We get perfect agreement between both approaches.

Show[plt4,plt5,plt6,plt7,plt8,plt9,DisplayFunction->$DisplayFunction]

Reference: Example 1.15 of Principles and Modern Applications of Mass Transfer Operations by Jaime Benitez, Wiley-Interscience, 2002.