By rolling your own Method Of Lines (MOL), you get more flexibility on what you can do with the derivatives. If you look at "vstar", then you see it was only evaluated at the interface
$u_0$ (note that part 1 of delconc is also at the interface).
vstar = -1/(1 - Subscript[u, 0][t]) delconc[[1]];
So, vstar is evaluated at the interface only where the flux of B is zero and it is only a function of
$t$ because its divergence is zero. Therefore, it should satisfy your conditions of only being evaluated at the interface and not throughout the domain.
Regarding the figure 20.1-1, if the characteristic length scale is,
$L$, and the characteristic time scale is
$\frac{L^2}{D_{AB}}$, then Z using dimensionless units becomes
$$Z=\frac{z}{\sqrt{4t}}$$
Given that figure 20.1-1 has a Z value that goes from 0 to 2 and we know that z is between 0 and 1, we can solve for t in terms of z.
Solve[z/(2*Sqrt[t]) == 2, t]
(* {{t -> z^2/16}} *)
So, I chose
$z=\frac{1}{4}<1$ and that meant the time when
$Z=2$ was
$t=\frac{1}{256}$. Then, I solved for
$z$ in terms of
$Z$ using.
Solve[bigZ == z/(2*Sqrt[t]), z]
(* {{z -> 2*bigZ*Sqrt[t]}} *)
Taken together, that means the following is the numerically derived
$X$ for
$x_{A0}=0.5$ @
$z=\frac{1}{4}$ (note that z in the Mathematica implementation is actually
$Z$ to avoid beginning with capitals).
f50[1/256, 2 Sqrt[1/256] z]/0.5
If we want to look at the case where
$z=1$, then
$t=\frac{1}{16}$. Now, instead of following the infinite tube, we expect the numerical result to satisfy the
$x_{A0}=0$ boundary condition.
Plot[{bigX[z, 0.75],
Callout[f75[1/16, 2 Sqrt[1/16] z]/0.75, "3/4", 1,
CalloutMarker -> "CirclePoint"], bigX[z, 0.5],
Callout[f50[1/16, 2 Sqrt[1/16] z]/0.5, "1/2", 0.9,
CalloutMarker -> "CirclePoint"], bigX[z, 0.25],
Callout[f25[1/16, 2 Sqrt[1/16] z]/0.25, "1/4", 0.8,
CalloutMarker -> "CirclePoint"], bigX[z, 0]}, {z, 0, 2},
PlotRange -> {0, 1}, GridLines -> Automatic,
Filling -> {{1 -> {2}}, {3 -> {4}}, {5 -> {6}}},
PlotStyle -> {Dashed, Thick, Dashed, Thick, Dashed, Thick, Dashed},
ImageSize -> Large, Frame -> True,
FrameLabel -> {Style["Z", 18], Style["X", 18],
Style[StringForm["Dimensionless Concentration Profiles"], 18]},
PlotLegends -> Placed[legend, {Right, Center}]]
It appears that the BC is satisfied.
Now, we can test your assertion that my model does not allow for mass conservation by simply evaluating
$N_{Bz}$ at
$z=1$. You can use the relations stated previously to derive an expression for the molar flux of B. Start with the dimensionless molar velocity.
$$\mathbf{v}^*=N_A+N_B$$
Use the expression for
$N_{B z}$ in terms of
$x_A$.
$$N_{B z} = \left (1-x_A \right ) \left(N_{A z}+N_{B z}\right)+\frac{\partial x_A}{\partial z} $$
Substitute
$\mathbf{v}^*$ for
$N_A+N_B$
$$N_{B z} = \left (1-x_A \right ) \mathbf{v}^*+\frac{\partial x_A}{\partial z}$$
Evaluate at the exit (
$z=1$).
$${N_{Bz}} = {\left. {\left( {1 - {x_A}} \right){{\mathbf{v}}^*} + \frac{{\partial {x_A}}}{{\partial z}}} \right|_{z = 1}}$$
Evaluate
$\mathbf{v}^*$ at the interface
$${N_{Bz}} = {\left. {\left( {\left( {1 - {x_A}} \right)\left( {{{\left. {\frac{1}{{1 - {x_{A,sat}}}}\frac{{\partial {x_A}}}{{\partial z}}} \right|}_{z = 0}}} \right) + \frac{{\partial {x_A}}}{{\partial z}}} \right)} \right|_{z = 1}}$$
We can implement and plot the last equation in Mathematica code for the
$x_{A,sat}=0.5$ case.
Plot[((1 - f50[t, 1]) (-D[f50[t, x], x]/(1 - 0.5)) /. {t -> tp,
x -> 0}) + ((D[f50[t, x], x]) /. {t -> tp, x -> 1}), {tp, 0.0002,
0.1}, PlotRange -> {-1, 40}, ImageSize -> Large, Frame -> True,
FrameLabel -> {Style["Dimensionless Time", 18],
Style["\!\(\*SubscriptBox[\(N\), \(Bz\)]\)", 18],
Style[StringForm["Flux of B at z=1"], 18]}]
Now, we did not achieve infinite flux because this is a numerical model, but the flux of B at the exit starts out very high and asymptotically approaches 0 at large t. So, the assertion of the model having zero B flux along the axis is false.
